Integrating Kernel Methods and Deep Neural Networks for Solving PDEs

Venue: ICLR_cc_2024_Workshop_AI4DiffEqtnsInSci
Published: 2024-02-07 05:29:02

Abstract: Physics-informed machine learning (PIML) has emerged as a promising alternative to conventional numerical methods for solving partial differential equations (PDEs). 
PIML models are increasingly built via deep neural networks (NNs) whose performance is very sensitive to the NN's architecture, training settings, and loss function. 
Motivated by this limitation, we introduce kernel-weighted Corrective Residuals (CoRes) to integrate the strengths of kernel methods and deep NNs for solving nonlinear PDE systems. To achieve this integration, we design a modular and robust framework which consistently outperforms competing methods in a broad range of benchmark problems. This performance improvement has a theoretical justification and is particularly attractive since we simplify the training process while negligibly increasing the inference costs. 
Our studies also indicate that the proposed approach considerably decreases the sensitivity of NNs to factors such as random initialization, architecture type, and choice of optimizer.

Main paper: 
 introduction: 
 Partial differential equations (PDEs) elegantly explain the behavior of many engineered and natural systems such as power grids (Brunton et al., 2016) , advanced materials (Mozaffar et al., 2019) , and biological agents (Brunton et al., 2016) . Since most PDEs cannot be analytically solved, numerical approaches such as the finite element method are frequently used to solve them. Recently, a new class of methods known as physics-informed machine learning (PIML) has been developed and successfully used in many applications.While PIML models have fueled a renaissance in modeling complex systems, their performance heavily depends on optimizing the model's training mechanism, loss function, and architecture (Pestourie et al., 2023) . To reduce the time and energy footprint of developing PIML models while improving their accuracy, we re-envision solving PDEs via machine learning (ML) and introduce Corrective Residuals (CoRes) that integrate the strengths of kernel methods and deep neural networks (NNs). 
related works: 
 We can broadly classify PIML models into two categories. The first group of methods relies on variants of neural networks (NNs) and can be traced back to Lagaris et al. (1998) . Physics-informed neural networks (PINNs) (Raissi et al., 2019; Sirignano & Spiliopoulos, 2018) and their extensions are the most widely adopted member of these methods and their basic idea is to parameterize the PDE solution via a deep NN. The parameters of this NN are optimized by minimizing a multicomponent loss function which encourages the NN to satisfy the PDE as well as the initial and/or boundary conditions (IC and BCs), see Figure A2 . This minimization is known to be very sensitive to the optimizer, loss function formulation, and NN's architecture. To decrease this sensitivity, recent works have developed adaptive loss functions (Chen et al., 2018) and tailored architectures that improve gradient flows or automatically satisfy BCs (Wang et al., 2021a; Lagaris et al., 1998) . These advancements, however, fail to generalize to a diverse set of PDEs and substantially increase the cost and complexity of training.The second group of PIML models leverage kernel methods which have long been used in ML but their application in solving PDEs is largely unexplored. The few existing works (see Zhang et al. (2022) ; Iwata & Ghahramani (2017) ; Meng & Yang (2023) ) exclusively employ zero-mean Gaussian processes (GPs) . With this choice, solving the PDE amounts to designing the GP's kernel whose parameters are obtained via either maximum likelihood estimation (MLE) or a regularized MLE where the penalty term quantifies the GP's error in satisfying the PDE system. In a recent work by Chen et al. (2021) , solving PDEs via a zero-mean GP is cast as an optimal recovery problem that aims to estimate the solution at a finite number of interior nodes in the domain. Once these values are estimated, the PDE solution is approximated anywhere in the domain via kernel regression. 
neural networks with corrective residuals: 
 Kernel methods such as GPs have less extrapolation and scalability powers compared to deep NNs. They also struggle to approximate PDE solutions that have large gradients or involve coupled dependent variables. However, GPs locally generalize better than NNs and are interpretable and easy to train. Grounded on these properties, we introduce deep architectures with kernel-weighted CoRes that integrate the attractive features of NNs and GPs for solving PDEs.We argue that the sole reliance on the kernel serves as a double-edged sword when solving PDEs. To demonstrate, we consider the task of emulating the function u(x) given the n samples X = {x 1 , • • • , x n } with corresponding outputs u = {u 1 , • • • , u n } where u i = u(x i ). If we endow u(x) with a GP prior with the mean function m(x; θ) and kernel c(x, x ′ ; ϕ), the conditional process is also a GP whose expected value at x * is:EQUATIONw := w(x * , X; ϕ) = c -1 (X, X; ϕ)c(X, x * ; ϕ) (1b) r := r(X, u; θ) = um(X; θ).Here, θ and ϕ are the model parameters (typically estimated via MLE), c(X, x * ; ϕ) = [c(x 1 , x * ; ϕ), • • • , c(x n , x * ; ϕ)] T , r denotes the residuals on the training data, w are the kernelinduced weights, and C = c(X, X; ϕ) is the covariance matrix with ij th entry c(x i , x j ; ϕ). The covariance function can be a deep NN as in Wilson et al. (2016) or the simple Gaussian kernel: c(x, x ′ ; ϕ) = exp -(xx ′ ) T diag(ϕ)(xx ′ ) .(2)Since many kernels can approximate an arbitrary continuous function (Rasmussen, 2006) , zeromean GPs are used in many regression problems as eliminating m(x; θ) reduces the number of trainable parameters while increasing numerical stability (see Appendix A 1 for more detailed discussions). Unlike regression, PDE systems cannot be accurately solved via zero-mean GPs without any in-domain samples since the posterior process in Equation (1) predicts zero for any point that is sufficiently far from the boundaries. This reversion to the mean behavior is due to the exponential decay of the correlations as the distance between two points increases, see Equation (2).Following the above discussions, we make two important observations on Equation (1): it heavily relies on m(x; θ) in data scarce regions and it regresses u regardless of the values of m(X; θ). These observations suggest that a GP with an NN mean function provides an attractive prior for solving PDEs since functions formulated as in Equation (1a) satisfy the BCs/IC and their smoothness can be controlled through the mean and covariance functions. This approach, however, presents two major challenges. First, the posterior distribution in this case will not be Gaussian in general since most practical PDEs are nonlinear. Second, jointly optimizing ϕ and θ is a computationally expensive and unstable process due to the repeated need for constructing and inverting C.We address the above challenges via modularization and formulating the training process based on maximum a posteriori (MAP) instead of MLE. Our framework consists of two sequential modules that aim to solve PDE systems with deep NNs that substantially benefit from kernel-weighted CoRes. These modules seamlessly integrate the best of two worlds: (1) the local generalization power of kernels close to the domain boundaries, and (2) the substantial capacity of deep NNs in learning multiple levels of distributed representations in the interior regions where there are no labeled data.In the first module, we endow the PDE solution with a GP prior whose mean and covariance functions are a deep NN and the Gaussian kernel in Equation ( 2), respectively. Conditioned on the data, u, the posterior distribution of the solution is again a GP and follows Equation (1a) where r and w denote the residuals and kernel-induced weights, respectively. Importantly, in this module we fix θ to some random values and choose ϕ such that the GP can faithfully reproduce u.In the second module, we obtain the final model by conditioning the GP on the (nonlinear) constraints that require the in-domain predictions to satisfy the PDE system. We achieve this conditioning by fixing ϕ from module 1 and optimizing θ to ensure that the model in Equation (1a) satisfies the PDE at n P DE randomly selected collocation points (CPs) in the domain, see Figure 1 .Model Characteristics: As extensively studied in Appendix A, our approach provides four unique features. First, the training cost in our approach mainly depends on the second module since selecting ϕ does not rely on MLE and is an inexpensive process. Additionally, our experiments consistently indicate that the performance of the final model is quite robust to ϕ as long as the BCs/IC are sufficiently sampled. This robustness is independent of the random values assigned to θ in module 1. Based on these two observations, in all of our experiments we simply assign 10 2 to all the kernel parameters and sample 40 points at each boundary.Second, the computational cost of coupling GPs and deep NNs in our framework is negligible during both training and testing since C does not change in the second module and its size only depends on u. When solving PDEs such as the Navier-Stokes equations that have multiple dependent variables, the size of u can grow rapidly since it will store boundary and initial data on multiple outputs.Hence, for such PDEs we decouple the kernel-weighted CoRes of the outputs to keep C and u small. As shown in Figure A5 , we formulate this decoupling by endowing the dependent variables with a collection of GP priors which share the same mean function but have independent kernels.The third feature of our model is its ability to exactly satisfy the BCs/IC as the number of sampled boundary points increases. This behavior is independent of the domain geometry and the potential noise corrupting the data; see Appendix A for the proof and examples. Due to this feature, the loss function in Figure 1 only minimizes the error in satisfying the PDE and excludes data loss terms.The above three features imply that training a PINN or its extensions costs similarly to the case where the same network is used as m(x; θ) in our framework. This behavior is very attractive since our approach consistently and substantially improves the performance of existing NN-based methods while also simplifying the training process.Lastly, our framework allows to perform data fusion and system identification with fast convergence rates by incorporating the additional measurements in the kernel structure in exactly the same way that BC/IC data are handled, see Appendix A 5.5 for details and results. 
results: 
 We compare the performance of our approach against four methods on the PDE systems given by Equations ( 11) to ( 14) in Appendix A 3. We study each problem under two settings to understand the effect of PDE complexity on the results. We obtain the reference solution u(x) for each PDE system as detailed in Appendix A 4.5 and use it to quantify the accuracy of the PIML models based on the Euclidean error norm, L 2 e , at n t = 10 4 randomly chosen points. Throughout, [x, y], t, and u denote the spatial coordinates, time, and the PDE solution, respectively. In the case of Navier-Stokes equation, u(x) = [u(x), v(x), p(x)] T denotes the solution.ℒ 𝑃𝐷𝐸 𝜽 = 1 𝑛 𝑃𝐷𝐸 σ 𝑖=1 𝑛 𝑃𝐷𝐸 𝜂 𝑡 + 𝜂𝜂 𝑥 -As detailed in Appendix A 4, we use a simple feed-forward NN in our framework and design its input and output dimensionality based on the PDE system. We denote our model via NN-CoRes and compare it against (1) GP OR of Chen et al. (2021) , (2) PINNs whose architectures are exactly the same as the mean function of NN-CoRes, (3) PINN DW which leverages dynamic weights for loss terms, and (4) PINN HC which is a PINN whose output is designed to strictly satisfy the BCs/IC.The results of our studies are summarized in Table 1 and indicate that our approach consistently outperforms other methods by relatively large margins. Interestingly, in most cases even the small NN-CoRes achieve lower errors than the high capacity version of the competing methods; indicating NN-CoRes more effectively use their networks' capacity to learn the PDE solution. To visually compare the capacity utilization across different NN-based models, in Figure A12 we provide the histogram of the PDE loss gradients with respect to θ at the end of training. We observe that NN-CoRes achieve the most near-zero gradients while satisfying the BCs/IC. In contrast, PINN HC , which is also designed to automatically satisfy the BCs/IC, struggles to minimize the PDE loss.We observe in Table 1 that the performance of all the methods drops as either the problem complexity increases (e.g., Burger's vs. LDC) or PDE parameters are changed to introduce nonlinearity (e.g., A = 3 vs A = 5 in LDC). This trend is expected as the architecture and training settings across our experiments are fixed. That is, we can increase the accuracy of all methods, especially NN-CoRes, by increasing their capacity or improving the training, see Appendix A 5.3 for multiple experiments. 
conclusions and future works: 
 We introduce kernel-weighted CoRes that integrate the strengths of kernel methods and deep NNs for solving nonlinear PDEs. We design a modular framework to achieve this integration and show that it improves the accuracy without complicating the training process. As extensively studied in Appendix A, our findings not only are very robust to the choice of optimizer and initial parameter values, but also applicable to various neural architectures.The current major limitation of our approach is that the contributions of the kernel-weighted CoRes decrease in the absence of boundary data. We believe devising periodic kernels is a promising direction for addressing this limitation which will be particularly useful in multiscale simulations where PDEs with periodic BCs frequently arise in the fine-scale analyses. 
 
 Appendix: We appreciate the support from the Office of the Naval Research, NASA's Space Technology Research Grants Program, and National Science Foundation. As demonstrated in Figure A1 our experiments involve sampling from a sinusoidal function where we study the effects of frequency, noise, data distribution, function differentiability, adopted prior mean function, and hyperparameter optimization on the behavior of GPs. For all of these studies we endow the GP with the following parametric kernel:EQUATIONwhere λ = [σ, ϕ, δ] T are the kernel parameters. In this equation, σ 2 is the process variance which, looking at Equation 1 in the main text, does not affect the posterior mean and hence we simply set it to 1 in our framework (this feature of our framework is in sharp contrast to other methods such as Chen et al. (2021) whose performance is quite sensitive to the selected kernel parameters). The rest of the parameters in Equation ( 3) are defined as follows. ϕ = 10 ω where ω is the length-scale or roughness parameter that controls the correlation strength along the x-axis, 1{•} returns 1/0 if the enclosed statement is true/false, and δ is the so-called nugget or jitter parameter that is added to the kernel for modeling noise and/or improving the numerical stability of the covariance matrix. We quantify the numerical stability of the covariance matrix via its condition number or κ.Given some training data, λ can be quickly estimated via maximum likelihood estimation (MLE).We denote parameter estimates obtained via this process by appending the subscript MLE to them, i.e., λ M LE . Alternatively, we can manually assign specific values to λ.We first study the effect noise by training two GPs where both GPs aim to emulate the same underlying function but one has access to noise-free responses while the other is trained on noisy data, see Figure A1 (a) and Figure A1 (b), respectively. We observe in Figure A1 (a) that the estimated value for δ M LE is very small since the data is noise-free (the small value is added to reduce κ) while in Figure A1 (b) the estimated nugget parameter is much larger and close to the noise variance (2.48e -3 vs. 2.50e -3). Additionally, comparing Figure A1 (a) and Figure A1 (c) we observe a direct relation between the frequency of the underlying function and the estimated kernel parameters. In particular, the magnitude of ω M LE increases as u(x) becomes rougher since the correlation between two points on it quickly dies out as the distance between those points increases (for this reason, ω is also sometimes called the roughness parameter). Further increasing the frequency of u(x) to the extend that it resembles a noise signal directly increases ω M LE . These points indicate that the kernel parameters of a GP are interpretable.We next study the reversion to the mean behavior and numerical instabilities of GPs in The reversion to the mean behavior is also seen in Figure A1 (e) but this time it is not undesirable since the functional form of the chosen parametric mean function is similar to u(x) (note that a large neural network can also reproduce the training data but such a network cannot match u(x) in interior regions where there are no labeled data). This similarity forces the kernel to regress residuals that are mostly zero (i.e., the kernel must regress a constant value in the entire domain). Since any two points on a constant function have maximum correlation, regressing such residuals requires ϕ → 0 which, in turn, renders the covariance matrix ill-conditioned to the extend that κ → +∞. Based on these observations, in our framework we do not estimate the kernel parameters jointly with the weights and biases of the deep neural network (NN).Lastly, in Figure A1 (e) we demonstrate that GPs can interpolate non-differentiable functions as long as they are provided with sufficient training data. The power and efficiency of GPs in learning from data is quite robust to the hyperparameters. As shown in Figure A1 (g) through Figure A1 (i) GPs with manually selected ω can accurately surrogate u(x) regardless of its frequency (the nugget value in these three cases is chosen such that κ does not exceed a pre-determined value). This attractive behavior forms the basis of our choice to manually fix ϕ in the first module of our framework. It is highlighted that the manual parameter selection results in sub-optimal prediction intervals but this issue does not affect our framework since we do not leverage these intervals.A 2 NEURAL NETWORKS WITH KERNEL-WEIGHTED CORRECTIVE RESIDUALS REPRODUCE THE DATAWe prove that the error of our model in reproducing the boundary data converges to zero as we increase the number of sampled boundary data. For the sake of completeness, we begin by a definition and invoking two theorems and then proceed with our proof.Reproducing kernel Hilbert space (RKHS): Let H be a Hilbert space of real functions u defined on an index set X . Then, H is called an RKHS with the inner product ⟨•, •⟩ H if the function c : X × X → R with the following properties exists:• For any x, c(x, x ′ ) as a function of x ′ is in H, • c has the reproducing property, that is ⟨u(x ′ ), c(x ′ , x)⟩ H = u(x). -1.0 -0.5 0.0 0.5 1.0 x -1.0 -0.5 0.0 0.5 1.0 (a) u(x) = sin (2πx) m(x) = 0 ωMLE = 0.551, δMLE = 4.182e-07, κ = 2.505e+07 -1.0 -0.5 0.0 0.5 1.0 x -1.0 -0.5 0.0 0.5 1.0 (b) u(x) = sin (2πx), ∼ N (0, 0.05 2 ) m(x) = 0 ωMLE = 0.703, δMLE = 2.482e-03, κ = 1.946e+07 -1.0 -0.5 0.0 0.5 1.0 x -1.0 -0.5 0.0 0.5 1.0 (c) u(x) = sin (4πx) m(x) = 0 ωMLE = 1.207, δMLE = 1.622e-07, κ = 1.176e+07 -1.0 -0.5 0.0 0.5 1.0 x -2 -1 0 1 2 (d) u(x) = sin (2πx) m(x) = 0 ωMLE = 1.182, δMLE = 5.244e-04, κ = 1.282e+04 -1.0 -0.5 0.0 0.5 1.0 x -1.0 -0.5 0.0 0.5 1.0 (e) u(x) = sin (2πx) m(x; θ) = θ sin (2πx), θMLE = 1.000 ωMLE = -7.829, δMLE = 5.197e-08, κ = ∞ -1.0 -0.5 0.0 0.5 1.0 x 0.0 0.2 0.4 0.6 0.8 1.0 (f) u(x) = | sin (2πx)| m(x) = 0 ωMLE = 1.330, δMLE = 1.567e-06, κ = 1.030e+07 -1.0 -0.5 0.0 0.5 1.0 x -1.0 -0.5 0.0 0.5 1.0 (g) u(x) = sin (πx) m(x) = 0 ω = 2.000, δMLE = 1.517e-05, κ = 3.323e+05 -1.0 -0.5 0.0 0.5 1.0 x -1.0 -0.5 0.0 0.5 1.0 (h) u(x) = sin (2πx) m(x) = 0 ω = 2.000, δMLE = 1.546e-05, κ = 3.323e+05 -1.0 -0.5 0.0 0.5 1.0 x -1.0 -0.5 0.0 0.5 1.0 (i) u(x) = sin (4πx) m(x) = 0 ω = 2.000, δMLE = 1.516e-05, κ = 3.323e+05 u(x) η(x) 95% Prediction Interval Training dataFigure A1: Properties of Gaussian processes: We demonstrate that GPs have interpretable hyperparameters and can regress a wide range of functions. The dependency of regression quality on the hyperparameters rapidly decreases as the size of the training data increases.Note that the norm of u is∥ u ∥ H = ⟨u, u⟩ H and that ⟨c(x, •), c(x ′ , •)⟩ H = c(x, x ′ ) since both c(x ′ , •) and c(x, •) are in H.Mercer's Theorem: The eigenfunctions of the real positive semidefinite kernel c(x, x ′ ) whose eignenfunction expansion with respect to measureπ is c(x, x ′ ) = N i=1 α i ψ i (x)ψ i (x ′ ), are or- thonormal. That is: ψ i (x)ψ j (x)dπ = δ ij (4)where δ ij denotes the Kronecker delta function. Following this theorem, we note that for a Hilbert space defined by the linear combinations of the eigenfunctions, that is uEQUATION)whereh = [h 1 , • • • , h n ]T is the observation vector, u(x) denotes the function that we aim to fit to h, Kimeldorf & Wahba (1971) .u = [u(x 1 ), • • • , u(x n )] T = [u 1 , • • • , u n ] T areIn our case, to prove that our model can reproduce the boundary data we first assume that the initial and boundary conditions are sufficiently smooth functions and that the neural network (i.e., the mean function of the GP) produces finite values on the boundaries. These assumptions simplify the proof by allowing us to work with the difference of these two terms.We now consider a specific form of Equation ( 5):EQUATION)where u(x) is the zero-mean GP predictor and ∥ u(x) ∥ H is the RKHS norm with kernel c(•, •). The second term on the right hand side corresponds to the negative log-likelihood of a Gaussian noise model with precision λ 2 and hence the minimizer of Equation ( 6) is the posterior mean of the GP Szeliski (1987) . Hence, we now need to show that as n → ∞ the minimizer of Equation ( 6), which is our GP, can reproduce the data h. We denote the ground truth function that we aim to discover and the variance around it by, respectively, h(x) and τ 2 (x) = (hh(x)) 2 dπ(h|x) where π(x, h) is the probability measure that generates the data (x i , h i ).We rewrite the second term on the right hand side of Equation ( 6) as:EQUATION)where the zero on the last line is due to the definition of h(x), i.e., h(x) = E[h|x]. Since τ 2 (x) is independent of u(x), we can use Equation ( 7) to rewrite Equation (6) as:EQUATIONWe now invoke Mercer's theorem to write u(x) = ∞ i=1 u i ψ i (x) and h(x) = ∞ i=1 h i ψ i (x)where ψ i are the eigenfunctions of the nondegenerage kernel of the GP. Since {ψ i } form an orthonormal basis, we can write:EQUATION)We take the derivative of Equation ( 9) with respect to u i and set it to zero to obtain:EQUATION)Since 1/nλ 2 → 0 as n → ∞, in the limit u i → h i , i.e., our zero-mean GP predictor corrects for the error that m(x, θ) has on reproducing the initial and boundary conditions. Note that the convergence in Equation ( 10) does not depend on τ 2 (x) and hence holds for the case where the observation vector h is noisy.A 3 DETAILS ON THE BENCHMARK PROBLEMS Below, we provide the governing equations as well as the initial and boundary conditions for the four PDE systems studied in this paper.We consider a viscous system subject to IC and Dirichlet BC in one space dimension:EQUATIONwhere x = [x, t] and ν is the kinematic viscosity. Equation ( 11) frequently arises in fluid mechanics and nonlinear acoustics. In our studies, we investigate the performance of different PIML models in solving Equation ( 11) for ν = 0.01 π , 0.02 π which controls the solution smoothness at x = 0 where a shock wave forms as ν approches zero.Nonlinear Elliptic PDE To assess the ability of our approach in learning high-frequency solutions, we study the boundary value problem developed in Chen et al. (2021) :EQUATIONwhere x = [x, y] and α = {20, 30} is a constant that controls the nonlinearity degree. f (x, y) is designed such that the solution is u(x, y) = sin (πx) sin (πy) + 2 sin (4πx) sin (4πy).Eikonal Equation We consider the two-dimensional regularized Eikonal equation Chen et al.(2021) which is typically encountered in the context of wave propagation:EQUATIONwhere x = [x, y] and ϵ = {0.01, 0.05} is a constant that controls the smoothing effect of the regularization term.Lid-Driven Cavity (LDC) The two-dimensional steady state LDC problem has become a gold standard for evaluating the ability of PIML models in solving coupled PDEs. This problem is governed by the incompressible Navier-Stokes equations:u x + v y = 0, ∀x ∈ (0, 1) 2 uu x + vu y = - 1 ρ p x + ν(u xx + u yy ), ∀x ∈ (0, 1) 2 uv x + vv y = - 1 ρ p y + ν(v xx + v yy ), ∀x ∈ (0, 1) 2 v(x, 0) = v(x, 1) = v(0, y) = v(1, y) = 0, ∀x, y ∈ [0, 1] u(x, 0) = u(0, y) = u(1, y) = 0, ∀x, y ∈ [0, 1] u(x, 1) = A sin (πx), ∀x ∈ [0, 1] p(0, 0) = 0 (14)where x = [x, y], ν = 0.01 is the kinematic viscosity, ρ = 1.0 denotes the density, and A = {3, 5} is a scaling constant. The Reynolds number for this LDC problem can be computed via Re = ρūL ν where ū = 1 0 A sin (πx)dx is the characteristic speed of the flow and L = 1 is the characteristic length. For the two cases A = {3, 5}, we obtain Re = {191, 318}.We first briefly introduce the four PIML models that we have used in our comparative studies and then provide some details on how the reference solution for each PDE system is obtained. These solutions are used to quantify the accuracy of the PIML models.To be able to directly compare the implementation of the four PIML models, we use Burgers' equation in the following descriptions. The PDE system is:EQUATIONwhere x = [x, t] are the independent variables, u is the PDE solution, and ν is a constant that denotes the kinematic viscosity. Also, we denote the output of the NN models via m(x; θ) throughout this section. Note that we also employ m(x; θ) for denoting the NN in the mean function of NN-CoRes. 15). To calculate L(θ), we first obtain the network's output at n BC points on the x = -1 and x = 1 boundaries, n IC points on the t = 0 boundary which marks the initial condition, and n P DE collocation points (CPs) inside the domain, see Figure A2b . For the n BC + n IC points on the boundaries, we can directly compare the network's outputs to the specified boundary and initial conditions in Equations ( 15b) and (15c). For each of the n P DE CPs, we evaluate the partial derivatives of the output and calculate the residual in Equation (15a). Once these three terms are calculated, we obtain L(θ) by summing them up as follows:EQUATIONThe loss function in Equation ( 16) is typically minimized via either the Adam Kingma & Ba (2014) or L-BFGS Liu & Nocedal (1989) methods which are both gradient-based optimization algorithms.With either Adam or L-BFGS, the parameters of the network are first initialized and then iteratively updated to minimize L(θ). These updates rely on partial derivaties of L(θ) with respect to θ which can be efficiently obtained via automatic differentiation Baydin et al. (2018) .While Adam and L-BFGS are both gradient-based optimization techniques, they have some major differences Sun et al. (2019) . Adam is a first-order method while L-BFGS is not since it is a quasi-Newton optimization algorithm. Compared to Adam, L-BFGS is more memory-intensive and has a higher per-epoch computational cost since it uses an approximation of the Hessian matrix during the optimization. Moreover, Adam scales to large datasets better than L-BFGS which does not accommodate mini-batch training. However, L-BFGS typically provides lower loss values and requires fewer number of epochs for convergence compared to Adam.One of the challenges associated with minimizing the loss function in Equation ( 16) is that the three terms on the right-hand side disproportionately contribute to L(θ). To mitigate this issue, a popular approach is to scale each loss component independently before summing them up, that is:EQUATION)Since the scale of the three loss terms can change dramatically during the optimization process, these weights must be dynamic, i.e., their magnitude must be adjusted during the training. In our experiments, we follow the process described in Wang et al. (2021b) for dynamic loss balancing and highlight that this approach is only applicable to cases where Adam is used.An alternative approach to dynamic weight balancing is to eliminate L BC (θ) and L IC (θ) from Equation ( 16) by requiring the model's output to satisfy the boundary and initial conditions by construction Berg & Nyström (2018) . To this end, we now denote the output of the network by m(x; θ) and then formulate the final output of the model as:EQUATIONwhere a(x) and b(x) are analytic functions that ensure m(x; θ) satisfies Equations ( 15b) and ( 15c) regardless of what m(x; θ) produces at x. A common strategy is to choose a(x) to be the signed distance function that vanishes on the boundaries and produces finite values inside the domain. The construction of b(x) is application-specific since one has to formulate a function that satisfies the applied boundary and initial conditions while generating finite values inside the domain. For the PDE system in Equation ( 15), one option is b(x) = -2 sin πx 1+e -t .A 4.4 OPTIMAL RECOVERY This recent approach leverages zero-mean GPs for solving nonlinear PDEs Chen et al. (2021) . Specifically, let us denote the kernel of a zero-mean GP via c(•, •). We associate c(•, •) with the reproducing kernel Hilbert space (RKHS) U where the RKHS norm is defined as ∥u∥. Following these definitions, we can approximate u(x) by finding the minimizer of the following optimal recovery problem:EQUATIONwhere n P DE , n BC , n IC are the number of nodes inside the domain, on the x = -1 and x = 1 lines where the boundary conditions are specified, and on the t = 0 line where the initial condition is specified, respectively. We denote the collection of these n P DE + n BC + n IC points via X.The optimization problem in Equation ( 19) is infinite-dimensional and hence Chen et al. variable z = z (1) , z (2) , z (3) , z (4) :EQUATIONwhere N = 4(n P DE + n BC + n IC ) + 3n P DE and Θ is the covariance matrix (see Section 3.4.1 of Chen et al. ( 2021) for details on Θ). Equation ( 20) can be reduced to an unconstrained optimization problem by eliminating the equality constraints following the process described in Subsection 3.3.1 of of Chen et al. (2021) . Once z is estimated, the PDE solution can be estimated at the arbitrary point x in the domain via GP regression.We note that the process of defining the slack variables and obtaining the equivalent finitedimensional optimization problem needs to be repeated for different PDE systems (e.g., in a PDE system one may have to define some of the slack variables as the Laplacian of the solution rather than the solution itself). Also, per the recommendations in Chen et al. ( 2021), c(•, •) is set to an anisotropic kernel and its parameters are chosen manually (i.e., they do not need to be jointly estimated with z) but, unlike our approach, this choice must be done carefully since it affects the results. In our studies, we use the values reported in Chen et al. ( 2021) for the kernel parameters.A 4.5 IMPLEMENTATION DETAILS IN OUR COMPARATIVE STUDIES Below, we describe the training procedure of the PIML models used throughout out paper and also comment on how the reference solutions are obtained for each PDE system. All of our codes, data, and models will be made publicly available upon publication.Training The NN-based approaches (i.e., NN-CoRes, PINN, PINN DW , and PINN HC ) are all implemented in PyTorch Paszke et al. (2019) and use hyperbolic tangent activation functions in all their layers except the output one where a linear activation function is used. The number and size of the hidden layers (see Table 1 in the main text) are exactly the same across these methods to enable a fair and straightforward comparison.To optimize NN-CoRes, PINNs, and PINN HC we leverage L-BFGS with a learning rate of 10 -2 while PINN DW is optimized using Adam with a learning rate of 10 -3 (note that the performance of L-BFGS deteriorates if dynamic weights are used in the loss function). To ensure these NN-based methods produce optimum models, we use a very large number of epochs during training. Specifically, we employ 1, 000 and 2, 000 epochs for single-and multi-output problems, respectively. Since Adam typically requires more epochs for convergence, we train PINN DW for 40, 000 epochs across all problems. To evaluate the loss function, we use 10, 000 collocation points within the domain in all cases. For PINN and PINN DW we uniformly sample boundary and/or initial conditions at 1, 000 locations while we only sample 40 points for NN-CoRes. This significant difference is due to the fact that we observed that NN-CoRes with just 40 boundary points can outperform other methods.Leveraging more boundary data improves the performance of NN-CoRes in solving PDE systems especially in satisfying the boundary and initial conditions.We fit GP OR based on the code and specifications provided by Chen et al. (2021) which leverages a variant of the Gauss-Newton algorithm for optimization. The performance of GP OR depends on the kernel parameters and the number of interior nodes n P DE where z needs to be estimated. For the former, we use the recommended values in Chen et al. (2021) and for the latter we choose two values (1, 000 and 2, 000) in our experiments.NN-CoRes, PINN, PINN DW , and PINN HC are trained on an NVIDIA GeForce RTX 3060 with 64 GB of RAM whereas GP OR is trained on a CPU equipped with a 11th Gen Intel-Core i7-11700K running at a base clock speed of 3.6 GHz.Reference Solutions We obtain the reference solutions for the PDE systems as follows:Table A1 : Summary of comparative studies for the LDC problem: We report L 2 e × 10 3 of different methods as a function of model capacity and A. The symbol ⊗ indicates the network architecture (e.g., 4⊗10 is an NN which has four 10neuron hidden layers). Unlike NN-based methods, GP OR 's accuracy relies on the number of interior nodes which we set to 1, 000 or 2, 000. GP OR is not applied to LDC as it relies on manual derivation of the equivalent variational problem which, unlike the first three PDEs, is not done by the developers Chen et al. (2021) . • Elliptic PDE: The analytical solution for this problem is u(x, y) = sin (πx) sin (πy) + 2 sin (4πx) sin (πy).• Eikonal Equation: We leverage the solution method provided by Chen et al. ( 2021) which applies the transformation u(x, y) = -ϵ log g(x, y) leading to the linear PDE g-ϵ 2 ∆g = 0 that can be solved via the finite difference method.• Lid-Driven Cavity: we use the finite element method implemented in the commercial software package COMSOL Multiphysics (1998).A 5 ADDITIONAL EXPERIMENTSIn the following subsections, we summarize the findings of some additional experiments that emphasize the attractive properties of NN-CoRes. We first provide some details on prediction errors that complement the results reported in the main text. Afterwards, we elaborate on the extended version of our framework that accommodates PDE systems such as the Navier-Stokes equations that have multi-variate solutions. Then, we conduct rigorous sensitivity analyses to characterize the effect of factors such as random initialization, noise (on data obtained from the initial and/or boundary conditions), optimization settings, and architecture on our results. These analyses demonstrate that our framework is substantially less sensitive to such factors compared to competing methods. We conclude this section by investigating the behavior of NN-CoRes' loss function during training and extending our approach for solving inverse problems.A 5.1 DETAILED ERROR ANALYSISThe solution of the LDC problem consists of three dependent variables which are the pressure p(x) and the two velocity components in the x and y directions, u(x) and v(x), respectively. In Table 1 we report the mean of the Euclidean norm of the error on the three outputs. In Table A1 we provide the errors for the individual outputs of this benchmark problem and observe the same trend where NN-CoRes consistently outperforms other methods. We also notice that all the models predict pressure with less accuracy compared to the velocity components. This trend is due to the facts that not only the scale of p(x, y) is smaller than the velocity components, but also p(x, y) is known at a single point on the boundaries whereas u(x, y) and v(x, y) are known everywhere on the boundaries.To gain more insight into the performance of each method, we visualize the error maps in Figure A3 . We observe that GP OR is least accurate either in regions with sharp solution gradients or inside the domain where boundary information is not effectively propagated inward by the zero-mean GP. For PINN DW , the errors are predominantly close to either the boundaries or where solution discontinuities are expected to appear. PINNs' errors in reproducing the BCs/IC is eliminated in PINN HC but at the expense of significant loss of accuracy elsewhere in the domain. These issues are largely addressed by NN-CoRes which reproduce BCs/IC and approximate high gradient solutions quite well.Figure A3 : Reference solutions and absolute error maps: Error maps of NN-CoRes are consistently smaller than the other three methods.Similar to Figure A3 , we provide the reference solution and error maps of different approaches for the Eikonal problem in Figure A4a where we observe similar patterns. Specifically, GP OR fails to properly propagate the boundary information inwards as it relies on a zero-mean GP. PINN is quite accurate inside the domain but cannot faithfully satisfy the boundary conditions. PINN DW addresses the inaccuracy of PINN close to the boundaries but incurs significant errors inside the domain as the reformulation in Equation ( 18) complicates the training dynamics. These issues are effectively addressed by NN-CoRes which achieve small errors inside the domain and on the boundaries.In Figure A4b we solve a canonical PDE system known as Helmholtz Wang et al. (2021b) which is defined as:EQUATIONIn Equation ( 21), q(x, y) is constructed such that the analytic solution is u(x, y) = sin(a 1 πx) sin(a 2 πy) where a 1 and a 2 are two constants that control the frequency along the x and y directions, respectively. The Helmholtz equation is a well-studied benchmark problem since PINNs fail to accurately solve it. To address this shortcoming, recent works have introduced quite complex architectures which typically leverage adaptive loss functions. We test our framework on this benchmark problem by setting a 1 = 1 and a 2 = 4 while using the same architecture and training procedure that are used in our comparative studies. As shown in Figure A4b our predictions accurately capture both the high-and low-frequency features of the solution. We note that the solution in Figure A4b is 5 times more accurate than the one reported in Wang et al. (2021b) which employs a considerably larger architecture (4 ⊗ 50) and leverages the adaptive loss function described in Equation (17).A 5.2 EXTENSION TO COUPLED SYSTEMS As schematically illustrated in Figure A5 we slightly modify our framework to solve coupled PDE systems such as the Navier-Stokes equations which have multiple dependent variables that interact with one another. The essential idea behind this modification is to endow each dependent variable with a GP prior. These GPs have independent kernels but a shared mean function that is parameterized via a deep neural network. While a single kernel can help in learning the inter-variable relations, we avoid this formulation for two main reasons. Firstly, it increases the size and condition number of the covariance matrix especially if the boundary conditions on these variables are significantly different. For instance, on the top edge (y = 1) in the LDC benchmark problem, pressure is unknown while the vertical and horizontal velocity components are equal to, respectively, zero Figure A5 : Solving the 2D incompressible Navier-Stokes equations for the lid-driven cavity problem: With minor architectural changes, our framework can also solve coupled PDE systems. Specifically, we endow each dependent variable with a GP prior. These GPs have independent kernels but a shared mean function that is parameterized via a deep neural network.and A sin (πx). Secondly, our empirical findings indicate that the shared mean function is able to adequately learn the hidden interactions between these dependent variables.A 5.3 SENSITIVITY ANALYSESIn this section, we conduct a wide range of sensitivity studies to assess the impact of factors such as random initialization, noise, network architecture, and optimization settings on the summary results reported in the main text.We first analyze the effect of the roughness parameter, ω, on the results. We use the simple Gaussian kernel in Equation (3) with σ 2 = 1 and ω = 2 in all of our studies. The nugget or jitter parameter of the kernel is chosen such that the covariance matrix is numerically stable. We ensure this stability by imposing an upper bound of approximately κ max ≈ 10 6 on the condition number of the covariance matrix, i.e., κ < κ max . This constraints typically results in a nugget value of around 10 -5 or 10 -4 . We have not optimized the performance of NN-CoRes with respect to κ max as we have found our current results to be sufficiently accurate.As stated in the main text, the performance of NN-CoRes is quite robust to the values chosen for ϕ = 10 ω as long as they lie within a certain range. To obtain this range, we conduct the following inexpensive experiment using the Burgers' equation and the extension of the kernel in Equation (3) to two-dimensional inputs, i.e., c x, x ′ ; ϕ, δ,σ 2 = 1 = exp -ϕ(x -x ′ ) 2 -ϕ(t -t ′ ) 2 + 1{x == x ′ }δ.We first sample n train equally spaced boundary samples using the provided analytic initial and boundary conditions. To quantify the effect of data size on the results, we consider 5 scenarios where n train ∈ {10, 20, 40, 80, 160}. For each of these five cases, we build 200 independent GPs whose only difference is the value that we assign to ω. Specifically, we consider 200 equally spaced values in the [-2, 6] range for ω and use each of these values in one of the GPs which all have a non-zero mean function (we use a deep NN whose parameters are randomly initialized and frozen as the mean function). Once these GPs are built, we use them to predict on n test = 10 4 boundary points (see Equation 1 in the main text for the prediction formula). The results of this study are shown in the left and middle plots in Figure A6a and indicate that as more training data are sampled on the boundaries a wider range of values for ω result in small test errors. We highlight that this study is computationally very fast since none of the GPs are optimized; rather their parameters are either chosen by us (i.e., ω), or fixed (i.e., δ, σ 2 , and parameters of the NN mean).Following the above study, we have decided to use 40 boundary points in NN-CoRes. Based on the left and middle plots in Figure A6a , ω = 2 seems to be a good choice (but not the optimum one) for minimizing the error in reproducing the initial and boundary conditions. To see the effect of this choice on the performance of a trained NN-CoRes, we again vary ω (50 equally spaced values in the [-2, 6] range) but this time we train an NN-CoRes model for each value of ω. We evaluate the performance of these models in solving the Burgers' equation by reporting the Euclidean norm of the error L e 2 at n test = 10 4 points randomly located in the domain. The results are shown in the right plot in Figure A6a and indicate that although ω = 2 is not the optimum choice, it yields a model whose performance is close to optimal (the optimum model is achieved via an ω close to 3).We now conduct a few extensive experiments to study the effect of network size and optimization settings on the performance of various NN-based models. First, we fix everything and increase the number of neurons in each hidden layer from 10 to 50 (at increments of 10) and solve the Burgers' and Elliptic PDEs via both NN-CoRes and PINNs. We then repeat this experiment but this time we fix the architecture to 4 ⊗ 20 and incrementally increase n P DE from 10 3 to 10 4 . The results of these two experiments are summarized in Figure A6b and indicate that NN-CoRes is much less sensitive to the problem than PINNs which perform quite well on Burgers' but fail at accurately solving the Elliptic PDE that has direction-dependent frequency. We also observe that NN-CoRes provide lower errors than PINNs in most simulations.In our next experiment, we study the effects of optimizer (L-BFGS vs Adam), random initialization, and architecture type on the performance of various models. To this end, we again consider the Burgers' and Elliptic PDE systems and solve them with six NN-based methods and GP OR . For each case we repeat the training process of each model 10 times to quantify the effect of random initialization on the models' solution accuracy. For these experiments, we also consider a new network architecture that we denote by M3 which is introduced in Wang et al. (2021b) and aims to improve gradient flows by designing feed-forward networks with connections that resemble transformers Vaswani et al. (2017) . In our framework, we replace the architecture that is used in all of our studies (which is a feed-forward neural network or an FFNN) with M3 and train the model with Adam (the resulting model is denoted by M3-CoRes). We also train another NN-based model denoted by M4 Wang et al. (2021b) whose architecture is the same as M3 but leverages dynamic weights in its loss function. We highlight that the simulations that leverage M3 as their architecture have more parameters (and hence learning capacity) than cases where FFNNs are used so we expect M3-based simulations to provide lower errors.The results of these simulations are summarized in Figure A6c and indicate that (1) NN-CoRes and GP OR are less sensitive to random initializations compared to PINNs and their variations, (2) unlike other models, NN-CoRes performs well in both PDE systems, i.e., our framework provides a more (c) Effect of optimizer, random initialization, and architecture type on errors for the Burgers' and Elliptic problems.2 × 10 -1 3 × 10 -1 4 × 10 -1 6 × 10 -1 L 2 e -u(x) 2 × 10 -1 3 × 10 -1 4 × 10 -1 6 × 10 -1L 2 e -v(x) transferable method for solving PDEs via machine learning, and (3) architectures besides simple FFNNs (such as M3) can also be used in our framework to achieve higher accuracy.10 0 3 × 10 -1 4 × 10 -1 6 × 10 -1 L 2 e -p(The above experiments are based on the Burgers' and Elliptic PDE problems but our studies indicate that similar trends appear in other problems. To demonstrate this, we solve the LDC problem via four NN-based models that either use L-BFGS or Adam as their optimizer. We repeat the training Comparing Figure A6c and Figure A6d we observe that the errors reported for the LDC problem are consistently larger than those reported for the Burgers' and Elliptic problems. The reason behind this trend is that not only LDC is a more complex problem where the PDE solution consists of three inter-dependent variables (compared to only one variable in the case of Burgers' or Elliptic PDEs), but also pressure is only known at a single point on the boundary (rather than everywhere on the boundary). To test the first assertion, we increase the networks sizes for both NN-CoRes and PINN DW and observe in Figure A7 that both models provide higher accuracy compared to the reported numbers in Table 1 . We highlight that the percentage of improvement is noticeably higher in the case of NN-CoRes and these errors keep reducing as the network size increases, as demonstrated in Figure A8 .Finally, we investigate the effect of noisy boundary data on our results. Specifically, we corrupt the solution values that we sample from the initial and/or boundary conditions before using them in our approach. We use a zero-mean normal distribution to model the noise and set the standard deviation to either 0.5% or 1% of the solution range. As shown in Figure A9 , the solution accuracy decreases as the noise variance increases (this trend is expected) but in all cases NN-CoRes are able to quite effectively eliminate the noise and solve the Burgers' and Elliptic PDE systems.To gain more insights into the training dynamics of our approach, we visualize the loss and accuracy during the training process in Figure A10 and Figure A11 where in the latter figure we track the errors individually for each output. We provide these plots for both PINNs and NN-CoRes where the loss function of the former is based on Equation ( 16) while NN-CoRes only use L P DE (θ) in their loss function. The solution accuracy is measured based on L 2 e and L 2 e 2 for points inside the domain and on its boundaries. Note that we square L 2 e on the boundaries to be able to directly see its contribution to PINNs' loss, see L BC (θ) in Equation ( 16). In the case of NN-CoRes, we also report the accuracy of its NN part on predicting the PDE solution to quantify the contributions of kernel-weighted CoRes towards the model's predictions.As it can be observed in Figures A10 and A11 , NN-CoRes typically converge faster than PINNs, see the plots whose y-axis title is L 2 e -Domain. We attribute this trend to the fact that, unlike in PINNs, the initial and boundary conditions are automatically satisfied in our models thanks to the kernel-weighted CoRes which are smooth functions. This features enables NN-CoRes to focus on satisfying the PDE system in module two of our framework while PINNs have to struggle with both the differential equations as well as the initial and boundary conditions. A10 and A11 is that the errors of NN-CoRes are consistently lower than their NN components both in the domain and on the boundaries. That is, the kernel-weighted CoRes positively contribute to the model's predictions both on the boundaries and inside the domain. This behavior is in sharp contrast to most approaches such as PINN HC that satisfy the boundary conditions at the expense of complicating the training process.Another interesting trend that we observe in Figures A10 and A11 is that PINNs achieve lower loss values than NN-CoRes in the case of Eikonal and LDC problems. While lower loss values are desirable, in these cases the observed trends are misleading. To explain this behavior, we note that the . loss function of NN-CoRes is simply L P DE (θ) as the boundary and initial conditions are automatically satisfied. However, the loss function of PINNs minimizes both L P DE (θ) and L BC (θ). That is, since PINNs do not strictly satisfy the boundary conditions, they are less regularized and hence can minimize L P DE (θ) (which dominates the overall loss) in a more flexible manner. However, this behavior provides less accuracy since the boundary conditions are not learnt sufficiently well.Finally, to visually compare the capacity utilization across different NN-based models, in Figure A12 we provide the histogram of the PDE loss gradients with respect to θ at the end of training. We observe that NN-CoRes achieve the most near-zero gradients while satisfying the BCs/IC. In contrast, PINN HC , which is also designed to automatically satisfy the BCs/IC, struggles to minimize the PDE loss. So far we have only used the differential equations along with the initial and boundary conditions in building NN-CoRes. In this section, we introduce an extension of our framework for solving inverse problems where (1) there are some (possibly noisy) labeled data available inside the domain (we refer to these samples as observations to distinguish them from data obtained from the initial and/or boundary conditions), and (2) one or more parameters in the differential equations are unknown.Our goal in such applications is to solve the PDE system while estimating the unknown parameters. As shown in Figure A13a , we modify our framework in two ways to solve the PDE system in Equation (15) assuming ν is unknown but u(x) is known at n obs random points in the domain. Specifically, we (1) use the n obs observations in the kernel of NN-CoRes in exactly the same way that the n BC + n IC boundary data are handled by the kernel, and (2) treat ν as one additional parameter that must be optimized along with the weights and biases of the NN.To evaluate the performance of our approach in solving inverse problems, we consider the Burgers', Elliptic, and Eikonal PDE systems. We solve each problem in two scenarios where there are either n obs = 100 or n obs = 200 observations available in the domain. As shown in Figure A13b , in all cases NN-CoRes can estimate the unknown PDE parameter quite accurately. The convergence rate in all cases is quite fast and insignificantly reduces as n obs is halved from 200 to 100.-0.001 0.000 0.001  


Decision: The reviewer mentions that this paper proposes an interesting approach to combine classical kernel methods with DNNs. In particular, the proposed method improves the performance of GPs and PINNs. Another advantage over classical PINNs that the authors are able to satisfy the initial and boundary conditions exactly rather than as soft constraints in the loss function. I think the comments from Reviewers 8ReM and p7T5 should be addressed in the camera-ready version and I vote for acceptance.

Reviews: review: title: title: An okay paper with critical experiments missing
main_review: main_review: This work talks about using kernels in physics-informed ML setting. The problem setup is interesting as it incorporates IC and BC together with the PDE. However, the reviewer has the following issues:
1. Lack of proper evaluation. It seems like only 0 Dirichlet BC is studied in the paper. The other cases such as non-zero Dirichlet, periodic, Neumann etc. should also be addressed, or at least mentioned.
2. It is not clear why robustness to optimizer is valid here. An intuitive reasoning would strengthen the work. Apart from Fig-A6c, there is no other mention that reviewer could find.
title: title: Nice paper and contribution
main_review: main_review: ## Summary
This paper proposes to merge kernel methods with DNNs with applications to solving PDEs. In particular, they propose Corrective Residuals (CoRes). 

## Strengths
- Nice theoretical analysis and guarantees of the proposed method
- Clear overview of SciML methods in the intro and references
- Nice experimental results on broad set of canonical PDEs including the hyperbolic Burgers and Eikonal equation.
- The paper is well-written.
- Very important that the proposed methods satisfies the BCs/IC as the number of sampled
boundary points increases and is independent of the domain geometry and the potential noise corrupting the data. Typically this is ignored in many ML approaches, so it is nice to see the authors include it with theoretical proofs.
- Very detailed appendix

## Weaknesses
- Please provide reference to the first sentence of sentence 3 that kernel methods such as GPs have less extrapolation and scalability powers compared to deep NNs.
- Baselines: I think kernel based Neural Operator methods, e.g., FNO (Li et al., 2020) and Multi-wavelet Neural Operator (Gupta et al., 2021) should be compared to as well since these methods also use a kernel approach within DNNs rather than just comparing to GPs and PINNs