RoPINN: Region Optimized Physics-Informed Neural Networks

Authors: haixu  wu
huakun  luo
yuezhou  ma
jianmin  wang
mingsheng  long
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-02 08:11:26

Abstract: Physics-informed neural networks (PINNs) have been widely applied to solve partial differential equations (PDEs) by enforcing outputs and gradients of deep models to satisfy target equations. Due to the limitation of numerical computation, PINNs are conventionally optimized on finite selected points. However, since PDEs are usually defined on continuous domains, solely optimizing models on scattered points may be insufficient to obtain an accurate solution for the whole domain. To mitigate this inherent deficiency of the default scatter-point optimization, this paper proposes and theoretically studies a new training paradigm as region optimization. Concretely, we propose to extend the optimization process of PINNs from isolated points to their continuous neighborhood regions, which can theoretically decrease the generalization error, especially for hidden high-order constraints of PDEs. A practical training algorithm, Region Optimized PINN (RoPINN), is seamlessly derived from this new paradigm, which is implemented by a straightforward but effective Monte Carlo sampling method. By calibrating the sampling process into trust regions, RoPINN finely balances optimization and generalization error. Experimentally, RoPINN consistently boosts the performance of diverse PINNs on a wide range of PDEs without extra backpropagation or gradient calculation. Code is available at this repository: https://github.com/thuml/RoPINN.

Main paper: 
 introduction: 
 Solving partial differential equations (PDEs) is the key problem in extensive areas, covering both engineering and scientific research [38, 40, 49] . Due to the inherent complexity of PDEs, they usually cannot be solved analytically [10] . Thus, a series of numerical methods have been widely explored, such as spectral methods [22, 40] or finite element methods [5, 7] . However, these numerical methods usually suffer from huge computational costs and can only obtain an approximate solution on discretized meshes [26, 43] . Given the impressive nonlinear modeling capability of deep models [4, 14] , they have also been applied to solve PDEs, where physics-informed neural networks (PINNs) are proposed and have emerged as a promising and effective surrogate tool for numerical methods [44, 36, 35] . By formalizing PDE constraints (i.e. equations, initial and boundary conditions) as objective functions, the outputs and gradients of PINNs will be optimized to satisfy a certain PDE during training [36] , which successfully instantiates the PDE solution as a deep model.Although deep models have been proven to enjoy the universal approximation capability, the actual optimization process of PINNs still faces thorny challenges [8, 24, 35] . As a basic topic of PINNs, the optimization problem has been widely explored from various aspects [17, 44] . Previous methods attempt to mitigate this problem by using novel architectures to enhance model capacity [58, 3, 50] , reweighting multiple loss functions for more balanced convergence [47] , resampling data to improve important areas [51] or developing new optimizers to tackle the rough loss landscape [55, 37] , etc. Orthogonal to the above-mentioned methods, this paper focuses on a foundational problem, which is the objective function of PINNs. We notice that, due to the limitation of numerical calculation, it is almost impossible to optimize the loss function in the complete continuous domain. Thus, the conventional PINN loss is only defined on a series of selected points [36] (Figure 1 ). However, the scatter-point loss function obviously mismatches the PDE-solving objective, which is approximating the solution on a continuous domain. This mismatch may fundamentally limit the performance of PINNs. Several prior works also try to improve the canonical PINN loss function, which can be roughly categorized into the following two paradigms. One paradigm enhances the optimization by adding high-order derivatives of PDEs as a regularization term to the loss function [55] . However, calculating high-order gradients is numerically unstable and time-consuming, even with automatic differentiation in well-established deep learning frameworks [2, 34] . The other paradigm attempts to bypass the high-order derivative calculation in the PINN loss function with variational formulations [18, 19, 20] . Nevertheless, these variational methods still face difficulties in calculating the integral of deep models and will bring extra computations, thereby mainly limited to very shallow models or relying on massive sampled quadrature points and elaborative test functions [11, 57] .This paper proposes and studies a new training paradigm for PINNs as region optimization. As shown in Figure 1 , we extend the optimization process from selected scatter points into their neighborhood regions, which can theoretically decrease the generalization error on the whole domain, especially for hidden high-order constraints of PDEs. In practice, we seamlessly transform this paradigm into a practical training algorithm, named Region Optimized PINN (RoPINN), which is implemented through simple but effective Monte Carlo sampling. In addition, to control the estimation error, we adaptively adjust the sampling region size according to the gradient variance among successive training iterations, which can constrain the sampling-based optimization into a neighborhood with low-variance loss gradients, namely trust region. In experiments, RoPINN demonstrates consistent and sharp improvement for diverse PINN backbones on extensive PDEs (19 different tasks) without any extra gradient calculation. Our contributions are summarized as follows:• To mitigate the inherent deficiency of conventional PINN optimization, we propose the region optimization paradigm, which extends the scatter-point optimization to neighborhood regions that theoretically benefits both generalization and high-order constraints satisfaction.• We present RoPINN for PINN training based on Monte Carlo sampling, which can effectively accomplish the region optimization. A trust region calibration strategy is proposed to reduce the gradient estimation error caused by sampling for more trustworthy optimization.• RoPINN can consistently improve the performance of various PINN backbones (i.e. canonical and Transformer-based) on a wide range of PDEs without extra gradient calculation. 
preliminaries: 
 A PDE with equation constraints, initial (ICs) and boundary conditions (BCs) can be formalized asEQUATIONwhere F, I, B denote the PDE equations, ICs and BCs respectively [6] . u : R d+1 → R m is the target PDE solution. x ∈ Ω ⊆ R d+1 represents the input coordinate, which is usually a composition of spatial and temporal positions, namely x = (x 1 , • • • , x d , t). Ω 0 corresponds to the t = 0 situation.Correspondingly, the PINN loss function (point optimization) is typically defined as follows [17, 36] :EQUATIONwhere u θ represents the neural network parameterized by θ. N Ω , N Ω0 , N ∂Ω are the numbers of sampled points in Ω, Ω 0 , ∂Ω respectively. λ * is the corresponding loss weight. Note that there is an additional data loss term in Eq. ( 2) when we can access the ground truth of some points [36] . Since we mainly focus on PDE constraints throughout this paper, we omit the data loss term in the above formalization, which is still maintained in our experiments. In this paper, we try to improve PINN solving by defining a new surrogate loss in place of the canonical definition of PINN loss in Eq. (2) .In contrast, the relevant literature mainly improves the objective function in two different directions as follows. Appendix F provides a more comprehensive discussion on other relative topics.High-order regularization The first direction is to add the high-order constraints of PDEs as regularization terms to the loss function [55] . Specifically, since PDEs are sets of identical relations, suppose that the solution u is a K-order differential function, Eq. ( 1) can naturally derive a branch of high-order equations, where the k-th derivative for the j-th dimension is∂ k ∂x k j F(u)(x) = 0,x ∈ Ω, 1 ≤ j ≤ (d + 1), 1 ≤ k ≤ K, corresponding to the following regularization:EQUATIONwhere N k,j denotes the number of sampled points with weight λ k,j . Although this design can explicitly enhance the model performance in satisfying high-order constraints, the calculation of highorder derivatives can be extremely time-consuming and unstable [39] . Thus, in practice, the previous methods [55, 31] only consider a small value of K. In the next sections, we will prove that RoPINN can naturally incorporate high-order constraints. Besides, as presented in Eq. 3, this paradigm still optimizes PINNs on scattered points, while this paper extends optimization to neighborhood regions.Variational formulation As a classical tool in traditional PDE solvers, the variational formulation is widely used to reduce the smoothness requirements of the approximate solution [42] . Concretely, the target PDEs are multiplied with a set of predefined test functions {v 1 , • • • , v M } and then the PDE equation term of the loss function is transformed as follows [18, 19, 20] :EQUATIONwhere F (xj ) defines the antiderivative of F on the j-th dimension. Using integrals by parts, the derivative operation in F is transferred to test functions {v k } M k=1 , thereby able to bypass high-order derivatives. However, the integral on Ω is still hard to compute, which requires massive quadrature points for approximation [18] . Besides, test function selection requires extra manual effort and will bring M times computation costs [57] . In contrast, RoPINN does not require test functions and will not bring extra gradient calculations. Also, RoPINN employs a trust region calibration strategy to limit the optimization in low-variance regions, which can control the estimation error of sampling.Question: Does the paper discuss the limitations of the work performed by the authors? Answer: [Yes] Justification: We have discussed our limitations in Appendix G. Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
method: 
 As aforementioned, we propose the region optimization paradigm to extend the optimization from scatter points to a series of corresponding neighborhood regions. This section will first present the region optimization and its theoretical benefits in both reducing generalization error and satisfying high-order PDE constraints. Then, we implement RoPINN in a simple but effective sampling-based way, along with a trust region calibration strategy to control the sampling estimation error.For clarity, we record the point optimization loss defined in Eq. ( 2) at x as L(u θ , x), where x ∈ Ω ∪ Ω 0 ∪∂Ω denotes the point selected from inner domain, initial state or boundaries. We adopt S to denote the finite set of selected points. Then Eq. ( 2) can be simplified asL(u θ , S) = 1 |S| x∈S L(u θ , x).Correspondingly, we define the objective function of our region optimization innovatively asEQUATIONwhere Ω r = [0, r] (d+1) represents the extended neighborhood region with hyperparameter r. Although this definition seems to require more sampling points than point optimization, we can develop an efficient algorithm to implement it without adding sampling points (see next section). Besides, this formalization also provides us with a convenient theoretical analysis framework. Next, we will discuss the theoretical properties of the two optimization paradigms. All proofs are in Appendix A.Generalization bound Here we discuss the generalization error in expectation [13] , which is independent of the point selection, thereby quantifying the error of PINN optimization more rigorously. Definition 3.1. The generalization error in expectation of model trained on dataset S is defined asEQUATIONwhere A denotes the training algorithm and A(S) represents the optimized model parameters. Assumption 3.2. The loss function L is L-Lipschitz and β-smooth with respect to model parameters, which means that ∀x ∈ Ω the following inequalities hold:EQUATION)Theorem 3.3 (Point optimization). Suppose that the loss function L is L-Lipschitz-β-smooth for θ.If we run stochastic gradient descent with step size α t at the t-th step for T iterations, we have that: [13, 52] ). ( 2) If L is bounded by a constant C for all θ, x and is non-convex for θ with monotonically nonincreasing step sizes α t ≤ 1 βt , then(1) If L is convex for θ and α t ≤ 2 β , then E gen ≤ 2L 2 |S| T t=1 α t (proved byE gen ≤ C |S| + 2L 2 (T -1)β(|S|-1) (tighter bound than [13, 52] ). Lemma 3.4. If L is bounded for all θ, x and is convex, L-Lipschitz-β-smooth with respect to model parameters θ, then L region r is also bounded for all θ, x and convex, L-Lipschitz-β-smooth for θ. Theorem 3.5 (Region optimization). Suppose that the point optimization loss function L is L-Lipschitz and β-smooth for θ. If we run stochastic gradient descent with step size α t for T iterations based on region optimization loss L region r in Eq. (5), the generalization error in expectation satisfies:EQUATION) If L is bounded by a constant C for all θ, x and is non-convex for θ with monotonically nonincreasing step sizes α t ≤ 1 βt , thenE gen ≤ C |S| + 2L 2 (T -1) β(|S|-1) -JL( |Ωr| |Ω| ) 2, where J is a finite number that depends on the training property at the several beginning iterations.Proof. Based on the Lipschitz assumption, E gen can be bounded by a term relating to the expectation of distance between parameter θ optimized from different training sets. The region optimization paradigm brings a more "consistent" gradient optimization direction than point optimization at each iteration, thereby benefiting the generalization property. See Appendix A.3 for complete proof.From Theorems 3.3 and 3.5, we can observe that region optimization can reduce the generalization error E gen . Furthermore, the region optimization theorem also provides a more general theoretical framework. For example, the conventional point optimization is equivalent to the case of Ω r = 0, where only one single point is selected for each region. For another extreme case, enlarging the region size to the whole domain (i.e. Ω r = Ω), Eq. ( 5) is equivalent to directly optimizing the loss defined on Ω, where the generalization error will be reduced to zero. Unfortunately, this ideal situation cannot be satisfied in practice, since Eq. ( 5) requires precise calculation of the integral over the whole domain. More discussions of practical implementation are deferred to the next section.High-order PDE constraints In our proposed region optimization (Eq. ( 5)), the integral operation on the input domain can also relax the smoothness requirements of the loss function L. For example, without any additional assumption of the smoothness of L(u θ , x) on x, we can directly derive the generalization error for the first-order loss ∂ ∂xj L region r (u θ , x) on the j-th dimension as follows.Input: number of iterations T , number of past iterations T 0 retained to estimate the trust region, default region size r, trust region calibration value σ 0 = 1, and initial PINN parameters θ 0 .Output: optimized PINN parameters θ T . Initialize an empty buffer to record gradients as g. for t = 0 to T do // Region Optimization with Monte Carlo Approximation Sample points from neighborhood regions:S ′ = {x i + ξ i } |S| i=1 , x i ∈ S, ξ i ∼ U [0, r σt ] (d+1) Calculate loss function L t = L (u θt , S ′ )Update θ t to θ t+1 with optimizer (Adam [21] , L-BFGS [27] , etc) to minimize loss function L t // Trust Region Calibration Record the gradient of parameters g t throughout optimization Update gradient buffer g by adding g t and keeping the latest T 0 elements Trust region calibration with σ t+1 = ∥σ(g)∥ end for Corollary 3.6 (Region optimization for first-order constraints). Suppose that L is bounded by C for all θ, x and is L-Lipschitz and β-smooth for θ. If we run stochastic gradient method based on first-order j-th dimension loss function ∂ ∂xj L region r for T iterations, the generalization error in Theorem 3.5(2) still holds when we adopt the monotonically non-increasing step size α t ≤ 1 2βt .Corollary 3.6 implies that the integral on the input domain in region optimization can help training PINNs with high-order constraints, which is valuable for high-order PDEs, such as wave equations. In contrast, this valuable property cannot be achieved by the classic point optimization. See Example 3.7. Example 3.7 (Point optimization fails in optimizing with first-order constraints). Under the same assumption with Corollary 3.6, we cannot obtain the Lipschitz and smoothness property of∂ ∂xj L(u θ , x). For example, suppose that L(u θ , x) = |θ T √ x|, x ∈ [0, 1] (d+1) , which is 1-Lipschitz- 1-smooth. However, ∇ θ ∂ ∂xj L(u θ , x) is unbounded when x → 0, thereby not Lipschitz constant.Derived from our theoretical insights of region optimization, we implement RoPINN as a practical training algorithm. As elaborated in Algorithm 1, RoPINN involves the following two iterative steps: Monte Carlo approximation and trust region calibration, where the former can efficiently approximate the optimization objective and the latter can effectively control the estimation error. Next, we will discuss the details and convergence properties of RoPINN. All proofs can be found in Appendix B.Monte Carlo approximation Note that the region integral in Eq. ( 5) cannot be directly calculated, so we adopt a straightforward implementation based on the Monte Carlo approximation. Concretely, to approximate the gradient descent on the region loss L region r , we uniformly sample one point within the region Ω r for the gradient descent at each iteration, whose expectation is equal to the gradient descent of the original region optimization in Eq. ( 5):EQUATIONIn addition to efficiently approximating region optimization without adding sampling points, our proposed sampling-based strategy is also equivalent to a high-order loss function, especially for the first-order term, which is essential in practice [55] . Concretely, with Taylor expansion, we have that:EQUATIONwhere Ω r = [0, r] d+1 , and L 1 represents the first order of loss function, namely ∂ ∂x L(u θ , x). Theorem 3.8 (Convergence rate). Suppose that there exists a constant H, s.t. ∀v and ∀x ∈ Ω, v T ∇ θ L region r (u θ , x)v ≤ H∥v∥ 2 . If the step size α t = 1 √ t+1 decreases over time for T iterations, the region optimization based on Monte Carlo approximation will converge at the speed ofEQUATION)Theorem 3.9 (Gradient estimation error). The estimation error of gradient descent between Monte Carlo approximation and the original region optimization satisfies:EQUATIONwhere σ represents the standard deviation of gradients in region Ω r .Trust region calibration Although the expectation of the Monte Carlo sampling is equal to region optimization as shown in Eq. ( 8), this design will also bring estimation error in practice (Theorem 3.9). A large estimation error will cause unstable training and further affect convergence. To ensure a reliable gradient descent, we propose to control the sampling region size r towards a trustworthy value, namely trust region calibration. Unlike the notion in optimization [56] , here trust region is used to define the area of input domain where the variance of loss gradients for different points is relatively small. Formally, we adjust the region size in inverse proportion to gradient variance:EQUATION)In practice, we initialize the trust region size as a default value r and calculate the gradient estimation error during the training process for calibration (Algorithm 1). However, the calculation of the standard deviation of gradients usually requires multiple samples, which will bring times of computation overload. In pursuit of a practical algorithm, we propose to adopt the gradient variance among several successive iterations as an approximation. Similar ideas are widely used in deep learning optimizers, such as Adam [21] and AdaGrad [48] , which adopt multi-iteration statistics as the momentum of gradient descent. The approximation process is guaranteed by the following theoretical results. Lemma 3.10 (Trust region one-iteration approximation). Suppose that loss function L is L-Lipschitz and β-smooth for θ and the t-th step parameter is θ t . Two gradient difference sequences between successive iterations, ∥∇ θ L(u θt , z 1 )-∇ θ L(u θt-1 , z 2 )∥ and ∥∇ θ L(u θt , z 1 )-∇ θ L(u θt , z 2 )∥, share the same limit, as the difference of the two sequences is dominated by the following inequality:EQUATIONwhere α t-1 represents the step size at the (t -1)-th iteration, which approaches 0 as t tends to ∞. Theorem 3.11 (Trust region multi-iteration approximation). Suppose that loss function L is L-Lipschitz and β-smooth for θ and the learning rate α t ≤ 1 βL converges to zero over time t, then the estimation error can be approximated by the variance of optimization gradients in multiple successive iterations. Given hyperparameter T 0 , our multi-iteration approximation is guaranteed byEQUATIONIt is worth noting that, as presented in Algorithm 1, since the gradient of each iteration has already been on the shelf, our design will not bring any extra gradient or backpropagation calculation in comparison with point optimization. Besides, our algorithm is not limited to a certain optimizer, and in general, we can effectively obtain the gradients of parameters by retrieving the computation graph.Balance between generalization and optimization Recall that in Theorem 3.5, we observe that a larger region size will benefit the generalization error, while Theorem 3.9 demonstrates that too large region size will also cause unstable training because it will result in excessive gradient estimation error of Monte Carlo sampling in our implementation. The above analysis reveals the underlying trade-off between generalization and optimization of PINN models, which is formally stated below. Theorem 3.12 (Region Optimization with gradient estimation error). Based on the same assumptions in Theorem 3.5 but optimizing the PINN model with the approximate region optimization loss Based on Theorem 3.12, we have pinpointed that classical point optimization (Ω r = 0) and sampling points globally (Ω r = Ω) discussed in Theorem 3.5 correspond to two extreme cases. The former makes E r,grad = 0 but cannot reduce the generalization error, while the latter holds a large gradient estimation error. Thus, neither case can yield perfect generalization for PINNs . In contrast, the design for calibrating trust regions in RoPINN provides an adaptive strategy to better balance generalization and optimization, which can adjust the region size according to multi-iteration training stability.L approx r (u θ , x) = ∇ θ L(u θ , x + ξ), ξ ∼ U (Ω r ) for T iterations,Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?• The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
experiments: 
 To verify the effectiveness and generalizability of our proposed RoPINN, we experiment with a wide range of PDEs, covering diverse physics processes and a series of advanced PINN models. Benchmarks For a comprehensive evaluation, we experiment with four benchmarks: 1D-Reaction, 1D-Wave, Convection and PINNacle [12] . The first three benchmarks are widely acknowledged in investigating the optimization property of PINNs [47, 37] . Especially, 1D-Reaction and Convection are highly challenging and have been used to demonstrate "PINNs failure modes" [24, 33] . As for PINNacle [12] , it is a comprehensive family of 20 tasks, including diverse PDEs, e.g. Burgers, Poisson, Heat, Navier-Stokes, Wave and Gray-Scott equations in 1D to 5D space and on complex geometries. In this paper, to avoid meaningless comparisons, we remove the tasks that all the methods fail and leave 16 tasks.Base models To verify the generalizability of RoPINN among different PINN models, we experiment with five base models, including canonical PINN [36] , activation function enhanced models: QRes [3] and FLS [50] , Transformer-based model PINNsFormer [58] and advanced physics-informed backbone KAN [28] . PINNsFormer [58] and KAN [28] are the most advanced PINN models.Baselines As stated before, this paper mainly focuses on the objective function of PINNs. Thus, we only include the gradient-enhanced method gPINN [55] and variational-based method vPINN [18] as baselines. Notably, there are diverse training strategies for PINNs focusing on other aspects than objective function, such as sampling-based RAR [51] or neural tangent kernel (NTK) approaches [47] .We also experimented with them and demonstrated that they contribute orthogonally to RoPINN.Implementations In RoPINN (Algorithm 1), we select the multi-iteration hyperparameter T 0 from {5, 10} and set the initial region size r = 10 -4 for all datasets, where the trust region size will be adaptively adjusted to fit the PDE property during training. For 1D-Reaction, 1D-Wave and Convection, we follow [58] and train the model with L-BFGS optimizer [27] for 1,000 iterations. As for PINNacle, we strictly follow their official configuration [12] and train the model with Adam [21] for 20,000 iterations. Besides, for simplicity and fair comparison, we set the weights of PINN loss as equal, that is λ * = 1 in Eq. ( 2). Canonical loss formalized in Eq. ( 2), relative L1 error (rMAE) and relative L2 error (rMSE) are recorded. All experiments are implemented in PyTorch [34] and trained on a single NVIDIA A100 GPU. See Appendix C for more implementation details.As shown in Table 2 , we investigate the effectiveness of RoPINN on diverse tasks and base models and compare it with two well-acknowledged PINN objectives. Here are two key observations.RoPINN can consistently boost performance on all benchmarks, justifying its generality on PDEs and base models. Notably, since the PDEs under evaluation are quite diverse, especially for PINNacle (Table 1 ), it is extremely challenging to obtain such a consistent improvement. We can find that the previous high-order regularization and variational-based methods could yield negative effects in many cases. For example, gPINN [55] performs badly on 1D-Wave, which may be due to second-order derivatives in the wave equation. Besides, vPINN [18] also fails in 1D-Reaction and QRes. As we stated in Table 1 , 1D-Reaction and Convection are hard to optimize, so-called "PINNs failure modes" [24, 33] . In contrast, empowered by RoPINN, PINNs can mitigate this thorny challenge to some extent. Specifically, with RoPINN, canonical PINN [36] , QRes [3] and FLS [50] achieve more than 90% improvements in 1D-Reaction. Besides, RoPINN can further enhance the performance of PINNsFormer [58] and KAN [28] , which have already performed well in 1D-Recation or Convection, further verifying its effectiveness in helping PINN optimization. Combining with other strategies Since RoPINN mainly focuses on the objective function design, it can be integrated seamlessly and directly with other strategies. As shown in Table 3 , we experiment with the widely-used loss-reweighting method NTK [47] and data-sampling strategy RAR [51] .Although NTK can consistently improve the performance, it will take extra computation costs due to the calculation of neural tangent kernels [15] . Based on NTK, our RoPINN can obtain better results with slightly more time cost. As for RAR, it performs unstable in different tasks, while RoPINN can also boost it. These results verify the orthogonal contribution and favorable efficiency of RoPINN w.r.t. other methods.To provide an intuitive understanding of RoPINN, we plot the curves of training statistics in Figure 2 , including temporally adjusted region size log( r σt ), train loss, and test performance. From Figure 2 (a), we can find that even though we initialize the region size as distinct values, RoPINN will progressively adjust the trust region size to similar values during training. This indicates that our algorithm can capture a potential "balance point" between training stability and generalization error, where the fluctuation of trust region size reveals the balancing process. Further, as shown in Figure 2(b-c ), if r is initialized as a value closer to the balance point (e.g. 1e-4 and 1e-5 in this case), then the training process will converge faster. And too large a region size (e.g. 1e-3) will decrease the convergence speed due to the optimization noise (Theorem 3.9). Number of sampling points in Eq. (8) For efficiency, RoPINN only samples one point within the trust region to approximate the region gradient descent. However, it is worth noticing that sampling more points will make the approximation in Eq. ( 8) more accurate, leading to a lower gradient estimation error. Further, since RoPINN employs an adaptive strategy to adjust region Ω r , a lower gradient estimation error will also make the optimization process adapt to a larger region size r. Therefore, we observe in Figure 3(a-b ) that sampling more points will also increase the finally learned region size r and speed up the convergence. In addition, Figure 3 (c) shows that adding sampled points can also improve the final performance, which has also been theoretically justified in Theorem 3.12 that the upper bound of generalization error is inversely proportion to gradient estimation error. Efficiency analysis As we discussed above, sampling more points can benefit the final performance, while we choose only sample one point as the default setting of RoPINN in the spirit of boosting PINNs without extra backpropagation or gradient calculation, which has already achieved significant promotion w.r.t. original PINNs (Table 2 ). To provide a more comprehensive understanding of algorithm property, we plot the efficiency-performance curve in Figure 4 , where we can obtain the following observations. Firstly, computation costs will grow linearly when adding points. Secondly, more points will bring better performance but will saturate around 10 points, where the performance fluctuations of 9, 13, and 30 points are within three times the standard deviations (Appendix D.3).Ablations To verify the effectiveness of our design in RoPINN, we present ablations in Figure 5 . It is observed that although we only sample one point, even fixed-size region optimization can also boost the performance of PINNs in most cases, demonstrating the effectiveness of introducing "region" to PINN optimization. However, as illustrated in Theorem 3.9, the sampling process may also cause gradient estimation error, so the relative promotion is inconsistent and unstable among different PDEs and base models. With our proposed trust region calibration, we can obtain a more significant and consistent improvement, indicating that achieving a better balance between optimization and generalization (formally stated in Theorem 3.12) performs an essential role in training PINN models.Loss landscape Previous research [24] has studied why PINN cannot solve the Convection equation and found that it is not caused by the limited model capacity but by the hard-to-optimize loss landscape.Here we also provide a loss landscape visualization in Figure 6 , which is obtained by perturbing the trained model along the directions of the first two dominant Hessian eigenvectors [24, 25, 54] .We can find that vanilla PINN optimized by PINN loss in Eq. ( 2) presents sharp cones. In contrast, empowered by RoPINN, the loss landscape is significantly smoothed. This visualization intuitively interprets why RoPINN can mitigate "PINN failure modes". See Appendix D for more results.  
conclusion: 
 This paper presents and analyzes a new PINN optimization paradigm: region optimization. Going beyond previous scatter-point optimization, we extend the optimization from selected points to their neighborhood regions. Based on this idea, RoPINN is implemented as a simple but effective training algorithm, where an efficient Monte Carlo approximation process is used along with a trust region calibration strategy to control the gradient estimation error caused by sampling, theoretically manifesting a better balance of generalization and optimization. In addition to theoretical advantages, RoPINN can consistently boost the performance of various PINN models without extra backpropagation or gradient calculation, demonstrating favorable efficiency, training stability and general capability.A Generalization Analysis in Section 3.1This section will present the proofs for the theorems in Section 3.1. 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA]Justification: Not applicable.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes] Justification: All the data and code are used with proper citation and license.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. 
 
 Appendix: The proof for the convex case is derived from previous papers [13, 52] under the Assumption 3.2. We derive a more compact upper bound for generalization error in expectation for the non-convex setting. Lemma A.1. Given two finite sets of selected pointsS = (x 1 , • • • , x N ) and S ′ = (x ′ 1 , • • • , x ′ N ), let S (i) = (x 1 , • • • , x i-1 , x ′ i , x i+1 , • • • , x N )be the set that is identical to S except the i-th element, the generalization error in expectation is equal to the expectation of the error difference between these two sets, which can be formalized as follows:EQUATIONProof. Directly deriving from the in-domain loss, we have:EQUATIONThen, according to Definition 3.1, we have:EQUATION)Lemma A.2 (Convex case). Given the stochastic gradient method with an update rule asG α,x (θ) = θ -α∇ θ L(θ, x) and L is convex in θ, then for α ≤ 2 β , we have ∥G α,x (θ 1 ) -G α,x (θ 2 )∥ ≤ ∥θ 1 -θ 2 ∥.Proof. For clarity, we denoteg = ∥∇ θ L(θ 1 , x) -∇ θ L(θ 2 , x)∥.Then we have:EQUATIONConvex setting Next, we will give the proof for the convex case of Theorem 3.3 (1) .Proof. According to Lemma A.1, we attempt to bound the generalization error in expectation E gen by analyzing the error difference between two selected sample sets. Denote S and S ′ as two identical sample sets of size |S| except for one sample. Suppose that with the stochastic gradient method on these two sets, we can obtain two optimization trajectories {θ t } T t=1 and {θ ′ t } T t=1 respectively. According to Assumption 3.2, we have the following inequality:EQUATIONWe assume two optimization trajectories, both obtained under the same random update rule and random permutation rule. Note that at the t-th step, with probability (1 -1 |S| ), the example selected by the stochastic gradient method is the same in both S and S ′ . As for the other 1 |S| probability, we have to deal with different selected samples. Thus, according to Lemma A.2, we have:EQUATIONDue to the L-Lipschitz assumption of L, the gradient ∇ θ L(u θ , x) is uniformly smaller than L, then:EQUATIONIn summary, since both optimization trajectories start from the same initialization, namely θ 0 = θ ′ 0 , the following inequality holds:EQUATIONFrom Lemma A.1, we haveE gen ≤ 2L 2 |S| T t=1 α t .Lemma A.3 (Non-convex case). Given the stochastic gradient method with an update rule asG α,x (θ) = θ -α∇ θ L(θ, x), then we have ∥G α,x (θ 1 ) -G α,x (θ 2 )∥ ≤ (1 + αβ)∥θ 1 -θ 2 ∥.Proof. This inequality can be easily obtained from the following:EQUATIONNon-convex setting Finally, we will give the proof for the non-convex case in Theorem 3.3(2).Proof. We also consider the optimization trajectory {θ t } T t=1 and {θ ′ t } T t=1 from S and S ′ , which are identical except for one element. We assume two optimization trajectories, both obtained under the same random update rule and random permutation rule. Let δ t = ∥θ t -θ ′ t ∥ and t 0 ∈ {1, • • • , |S|} be a considered iteration. Here, t ≤ |S| because for t > |S|, we must have δ t0 ̸ = 0. Then we have:EQUATIONSimilar to the convex case, we analyze the expectation of parameter difference in the (t + 1)-th iteration as follows. Since α ≤ 1 βt , then we have:EQUATIONEQUATIONAccumulating the above in equations recursively, we have:EQUATIONOrganizing the above inequalities, we have:EQUATIONAccording to Eq. ( 24), for arbitrary T ≥ 1, we just choose t 0 = 1, thenEQUATION)This boundary is tighter than [13, 52] , where the latter omits the 2L β(|S|-1) term in Eq. (27) .A.2 Proof for Properties of Region Optimization (Lemma 3.4)Since L region r is defined as a region integral of L, Lemma 3.4 can be easily obtained by:EQUATIONA.3 Proof for Region Optimization Generalization Error (Theorem 3.5)Similar to the proof in Appendix A.1, we will discuss the generalization error on region optimization.Convex setting Firstly, we would like to prove the convex case.Proof. According to Lemma 3.4, the region loss also holds the convexity, Lipschitz and smoothness properties, which ensures that Lemma A.2 and A.3 still work for L region r . We also focus on the selected sample sets S and S ′ , which are identical except for one element. Thus, at t-step, the following equation is satisfied:EQUATIONwhereG region α,x (θ) = θ -α∇ θ L region r (θ, x).As for the second item on the right part, we also consider the upper bound of ∇ θ L region r (u θ , x). However, different from point-wise optimization, x + Ω r could overlap with x ′ + Ω r . For clarity, we define the overlapped area as Ω in , whose size is larger than zero when x ′ fall into the area centered at x with size 2 (d+1) |Ω r |. Actually, due to the boundary of Ω, we cannot always ensure (x ′ + Ω r ) ⊂ Ω. Thus, for simplification, we assume that the domain Ω can be projected to a torus, where the out-of-domain samples will be re-included to Ω.Further, E x,x ′ ∈Ω, s.t. |Ωin|=0 is simplified as E Ωin=0 and E x,x ′ ∈Ω, s.t. |Ωin|>0 is shorted as E Ωin>0 .And the operator I is defined as I(x) = max (0, min(1, x)) Thus, we can obtain the estimation for the difference between the updated model parameters through the following derivations:EQUATIONIn the above inequalities, the third inequality is based on the following derivations. Firstly, denote g =∥ 1 |Ωr| Ωin ∇ θ L(u θt , x + ξ)dξ -1 |Ωr| Ωin ∇ θ L(u θ ′ t , x + ξ)dξ ∥. According to Assumption 3.2, we have g ≤ |Ωin| |Ωr| β∥θ t -θ ′ t ∥ ≤ β∥θ t -θ ′ t ∥.Thus, the following inequality holds:EQUATIONThus, recursively accumulating the residual at the t-th step, we have:EQUATIONFor the more general case, we no longer assume that the domain Ω can be projected to a torus, resulting in a non-symmetric scenario when |Ω in | > 0. This asymmetry arises due to the presence of boundaries, since points that could potentially intersect with the set x + Ω r may be truncated by the boundary. Specifically, we consider Ω = [0, l] (d+1) and Ω r = [0, r] (d+1) . The concrete probability P(|Ω in | > 0) and the expectation E |Ωin|>0 (|Ω in |) can be calculated as follows:EQUATIONThus, the general case of Eq. 31 can be reformulated using the identities above. We assume r l < 0.5, as when r l ≥ 0.5, it follows that P(|Ω in | > 0) = 1. Specifically, Eq. ( 31) can be rewrite as:EQUATIONThus, recursively accumulating the residual at the t-th step, we have:EQUATIONAlthough the specific forms differ and the general case is much more complex, these two inequalities both share the same intuitive meaning: region optimization benefits from the overlap between x + Ω r and x ′ + Ω r . Moreover, within a certain range, the benefit increases as the value of r becomes larger.Thus, to keep the bound simple and easy to understand, the main text theorems are under the assumption that Ω can be projected to a torus. Otherwise, |Ωr| |Ω| should be replaced by ( r 2 (2l-3r)(3l-4r)3(l-r) 4) (d+1) .Non-convex setting Next, we will prove the non-convex setting. Similarly, we assume Ω can be projected to a torus. Otherwise, the |Ωr| |Ω| in the final bound should be replaced by ( r 2 (2l-3r)(3l-4r)3(l-r) 4) (d+1) .Proof. For clarity, let M = |Ωr| |Ω| . Then, we can rewrite the Eq. ( 25) as follows.If E(δ t ) ≤ 2L β , we have:EQUATIONwhere the second inequality is based on the following derivations:EQUATIONNotably, E ∥G region αt,x (θ t ) -G region αt,x ′ (θ ′ t )∥ has an obvious upper bound, i.e. E [∥θ t -θ ′ t ∥] + 2L βt . And only when E(δ t ) ≤ 2L β , the bound derived by Eq. ( 38) is tighter. Furthermore, the condition that E(δ t ) ≤ 2L β can be easily satisfied at the beginning several iterations since E(δ 0 ) = 0. Otherwise, we still take the following equation:EQUATIONwhere we do not consider the benefits brought by the overlap area of x + Ω r and x ′ + Ω r .Suppose that at the first K steps E[δ t0+K ] ≤ 2L β , Accumulating the above in equations recursively, we have the generalization error bound accumulated to the first K steps as follows:EQUATION)where J is a finite value that depends on the training property of beginning iterations, namely K and t 0 . The last inequality is from( t0+K t ) M |S| ≤ (1 + M ), when |S| is sufficient enough.Then, considering the all T steps, we haveEQUATION) Thus, following a similar proof process as Theorem 3.3(2), we can obtain:EQUATION)With t 0 = 1, we have the generalization error under the non-convex case satisfies:EQUATIONA.4 Proof for High-order Constraint Optimization (Corollary 3.6)First, we would like to prove the following Lemma. Lemma A.4. Suppose that L is bounded by C for all θ, x and is L-Lipschitz and β-smooth for θ, then the first-order j-th dimension loss function ∂ ∂xj L region r is also bounded by C for all θ, x and is 2L-Lipschitz and 2β-smooth for θ.Proof. For the bounded property, with the non-negative property of loss function, we have:EQUATIONwhereξ r = (• • • , r, • • • ) ∈ Ω t \x j and ξ 0 = (• • • , 0, • • • ) ∈ Ω t \x j .As for the Lipschitz and smoothness, we can obtain the following inequalities:EQUATIONThus, ∂ ∂xj L region r is also bounded by C for all θ, x and is 2L-Lipschitz and 2β-smooth for θ.Next, we will give the proof for Corollary 3.6.Proof. According to Lemma A. EQUATIONLet M = |Ωr| |Ω| , since α t ≤ 1 2βt , we can rewrite the Eq. ( 37) as follows:EQUATIONwhere the second term is derived from Eq. ( 38) by substituting L to 2L and β to 2β, which is:EQUATIONThus, following the same derivation as Theorem 3.5, we have Corollary 3.6 holds.This contains the proof for the theoretical analysis of our proposed algorithm in Section 3.2.The crux of proof is to take expectation for Monte Carlo sampling.Proof. From Taylor expansion, there exist x ′ such that:EQUATION)Taking expectations to ξ on both sides, since E[∇ θ L(u θt , x + ξ)] = ∇ θ L region r (u θt , x + ξ), we have:EQUATIONRearranging the terms and accumulating over T iterations, we have the following sum:EQUATIONwhere u * represents the global optimum. Here we run the gradient descent for a random number of iterations τ . For τ = t iterations with probability:EQUATIONThus, with α t = 1 √ t+1 , we have the gradient norm is bounded by:EQUATIONB.2 Proof for Estimation of RoPINN (Theorem 3.9)As presented in Eq. ( 8) and ( 50), we approximate the region optimization with the Monte Carlo sampling method. For better efficiency, we propose only to sample one point at each iteration. However, this will cause an estimation error formalized in Theorem 3.9, which can be directly derived by the definition of standard deviation as follows:Proof. According to the definition of L region r in Eq. ( 5), we get First, we give the proof for Lemma 3.10.EQUATIONProof. According to Assumption 3.2, there exist x ′ , such that the following equation holds:EQUATIONThus, the following inequality holds:EQUATIONNext, we will prove Theorem 3.11.Proof. This theorem can be proved by demonstrating that: for all i, j ∈ {1, • • • , T 0 }:EQUATIONFor the first equation, since α t → 0, given ∀ϵ, there exists a constant M such that any t > M , α t ≤ ϵ T0Lβ . Thus, for any t > M and any i, j ∈ {1, • • • , T 0 }, the following equation is satisfied:EQUATION)Thus,lim t→∞ ∇ θ L(u θt-i+1 , z i ) = ∇ θ L(u θt , z i ). Therefore, given ∀ϵ ′ , there exist a constant M ′ , ∀t > M ′ , ∥∇ θ L(u θt-i+1 , z i ) -∇ θ L(u θt , z i )∥ ≤ ϵ 2 .As for the second equation, for any t > M ′ , the following equation is satisfied:EQUATION)Thus, Theorem 3.11 can be proved by replacing the gradient of past iterations with their limitations.Optimization with Gradient Estimation Error (Theorem 3.12)Convex setting Firstly, we would like to prove the convex case as follows.Proof. Similar to the proof of region optimization in Appendix A.3, at the t-th step, we can obtain the following equation:EQUATIONwhereG approx α,x (θ) = θ -α∇ θ L approx r (θ, x) = θ -α∇ θ L(θ, x + ξ), ξ ∼ U (Ω r ). Suppose that we have sampled ξ, ξ ′ ∈ Ω r , the second term on the right part can be bounded as follows:EQUATIONEQUATIONThus, recursively accumulating the residual at the t-th step, we have:EQUATION)Non-convex setting Similarly, we can prove the non-convex setting as follows.Proof. It is easy to prove that L approx r (θ, x) is still L-Lipchitz-β-smoothness for θ. For clarity, we define that M = |Ωr| |Ω| . Thus, based on Eq. ( 38), we have the following derivations.If E(δ t ) ≤ 2L β -2 βM E r,grad , we have:EQUATIONOtherwise, we still consider the following inequality:EQUATIONSuppose that at the first K ′ steps E[δ t0+K ′ ] ≤ 2L β -2 βM E r,grad . Accumulating the above in equations recursively, we have the generalization error bound accumulated to the first K ′ steps as follows:EQUATIONwhere J ′ is a finite value that depends on the training property of beginning iterations, namely K ′ and t 0 . The last inequality is from( t0+K ′ t ) M |S| ≤ (1 + M ), when |S| is sufficient enough.Then, considering the all T steps, we haveE [∥θ T -θ ′ T ∥|δ t0 = 0] ≤ ∆ + T t=t0+K+1 Π T k=t+1 exp 1 t - 1 t|S| 2L βt|S| ≤ T t=t0+1 exp (1 - 1 |S| ) log T t 2L βt|S| -J ′ M 2 + J ′ E r,grad (1 + M ) ( T k=t+1 1 k ≤ log T t ) = 2L β|S| T 1-1 |S| T t=t0+1 t -(1-1 |S| )-1 -J ′ M 2 + J ′ E r,grad (1 + M ) ≤ 2L β|S| T 1-1 |S| 1 1 -1 |S| t -(1-1 |S| ) 0 -T -(1-1 |S| ) -J ′ M 2 + J ′ E r,grad (1 + M ). (Integral approximation) (66)Next, following the proof in Appendix A.3, we can obtain the generalization bound as follows:EQUATION)Note that K ′ does not exist when 2L β < 2 βM E r,grad , then J ′ = 0, which corresponds to the situation that the region size is too large and brings serious gradient estimation error. Introducing "region" cannot bring a better generalization bound in this case.This section provides experiment details, including benchmarks, metrics and implementations.To comprehensively test our algorithm, we include the following four benchmarks. The first three benchmarks cover three typical PDEs (plotted in Figure 7 ), which are widely used in exploring the PINN optimization [24, 37] . The last one is an advanced comprehensive benchmark with 20 different PDEs. Here are the details. This problem is a one-dimensional non-linear ODE, which describes the chemical reactions. The concrete equation that we studied here can be formalized as follows:EQUATIONThe analytical solution to this problem is u(x, t) = h(x)e ρt h(x)e ρt +1-h(x) and h(x) = exp -(x-π) 2(π/4) 2 . In our experiments, we set ρ = 5. This problem is previously studied as "PINN failure mode" [24] , which is because of the non-linear term of the equation [29] . Besides, as shown in Figure 7 (a), it contains sharp boundaries for the center high-value area, which is also hard to learn for deep models.Following experiments in PINNsFormer [58] , we uniformly sampled 101 points for initial state Ω 0 and boundary ∂Ω and a uniform grid of 101×101 mesh points for the residual domain Ω. For evaluation, we employed a 101×101 mesh within the residual domain Ω. This strategy is also adopted for 1D-Wave and Convection experiments.1D-Wave This problem presents a hyperbolic PDE that is widely studied in acoustics, electromagnetism, and fluid dynamics [1] . Concretely, the PDE can be formalized as follows:EQUATIONThe analytic solution for this PDE is u(x, t) = sin(πx) cos(2πt) + 1 2 sin(βπx) cos(2βπt). We set β = 3 for our experiments. As presented in Figure 7 (b), the solution is smoother than the other two datasets, thereby easier for deep models to solve in some aspects. However, the equation contains second-order derivative terms, which also brings challenges in automatic differentiation. That is why gPINN [55] fails in this task (Table 2 ). Table 4 : Details of datasets in PINNacle [12] (16 different PDEs included in our experiments), including the dimension of inputs, highest order of PDEs, number of train/test points and concrete equations. Here we only present the simplified PDE formalizations for intuitive understanding. More detailed descriptions of PDE type and coefficient meanings can be found in their paper [12] . -∆u + k 2 u = f (x, y) 3d-CG 3D 2 49152 40960 -µi∆u + k 2 i u = f (x, y, z), i = 1, 2 2d-MS 2D 2 12288 10329 -∇(a(x)∇u) = f (x, y) Heat 2d-VC 2D+Time 2 65536 49189 ∂u ∂t -∇(a(x)∇u) = f (x, t) 2d-MS 2D+Time 2 65536 49189 ∂u ∂t - 1 (500π) 2 uxx -1 π 2 uyy = 0 2d-CG 2D+Time 2 65536 49152 ∂u ∂t -∆u = 0 NS 2d-C 2D 2 14337 12378 u • ∇u + ∇p -1 Re ∆u = 0, ∇ • u = 0 2d-CG 2D 2 14055 12007 Wave 1d-C 1D+Time 2 12288 10329 utt -4uxx = 0 2d-CG 2D+Time 2 49170 42194 ∇ 2 -1 c(x) ∂ 2 ∂t 2 u(x, t) = 0 Chaotic GS 2D+Time 2 65536 61780 ut = ε1∆u + b(1 -u) -uv 2 vt = ε2∆v -dv + uv 2 High- PNd 5D 2 49152 67241 -∆u = π 2 4 n i=1 sin 2 xi dim HNd 5D+Time 2 65537 49152 ∂u ∂t = k∆u + f (x, t)Convection This problem is also a hyperbolic PDE that can be used to model fluid, atmosphere, heat transfer and biological processes [41] . The concrete PDE that we studied in this paper is:EQUATIONThe analytic solution for this PDE is u(x, t) = sin(x -βt), where β is set as 50 in our experiments. Note that although the final solution seems to be quite simple, it is difficult for PINNs in practice due to the highly complex and high-frequency patterns. And the previous research [24] has shown that the loss landscape of the Convection equation contains many hard-to-optimize sharp cones.PINNacle This benchmark [12] is built upon the DeepXDE [30] , consisting of a wide range of PDEs and baselines. In their paper, the authors included 20 different PDE-solving tasks, covering diverse phenomena in fluid dynamics, heat conduction, etc and including PDEs with high dimensions, complex geometrics, nonlinearity and multiscale interactions. To ensure a comprehensive evaluation, we also benchmark RoPINN with PINNacle.During our experiments, we found that there are several subtasks that none of the previous methods can solve, such as the 2D Heat equation with long time (Heat 2d-LT), 2D Navier-Stokes equation with long time (NS 2d-LT), 2D Wave equation with long time (Wave 2d-MS) and Kuramoto-Sivashinsky equation (KS). In addition to the challenges of high dimensionality and complex geometry mentioned by PINNacle, we discover unique challenges in these tasks caused by long periods and high-order derivatives of governed PDEs, making them extremely challenging for current PINNs. To solve these problems, we might need more powerful PINN backbones. Since we mainly focus on the PINN training paradigm, we omit the abovementioned 4 tasks to avoid the meaningless comparison and experiment with the left 16 tasks. Our datasets are summarized in Table 4 .In our experiments, we adopt the following three metrics. Training loss, rMAE and rMSE. And the training loss has been defined in Eq. ( 2). Here are the calculations for rMSE and rMAE: rMAE:EQUATIONwhere u * denotes the ground truth solution. Note that the model output and ground truth can be negative and positive, respectively. Thus, these two metrics could be larger than 1.For classical base models PINN [36] , QRes [3] and FLS [50] , we adopt the conventional configuration from previous papers [58] . As for the latest model PINNsFormer [58] and KAN [28] , we use their official code. Next, we will detail the implementations of optimization algorithms.RoPINN As we described in the main text, we set the initial region size r = 10 -4 , past iteration number T 0 ∈ {5, 10} and only sample 1 point for each region at each iteration for all datasets. The corresponding analyses have been included in Figure 2 for r, Figure 3 for sampling points and Appendix D.1 for T 0 to demonstrate the algorithm property under different hyperparameter settings.In addition, our formalization for region optimization in Eq. ( 8) only involves the equation, initial and boundary conditions, where we can still calculate their loss values after random sampling in the extended region. This definition perfectly matches the setting of 1D-Reaction, 1D-Wave and Convection. However, in PINNacle [12] , some tasks also involve the data loss term, such as the inverse problem (Appendix D.2), which means we can only obtain the correct values for several observed or pre-calculated points. Since these points are pre-selected, we cannot obtain their new values after sampling. Thus, we do not apply region sampling to these points in our experiments. Actually, the data loss term only involves the forward process of deep models, which is a pure data-driven paradigm and is distinct from the other PDE-derived terms in PINNs. Therefore, the previous methods, gPINN and vPINN, also do not consider the data loss term in their algorithms.gPINN For the first three benchmarks, we add the first-order derivatives for spatial and temporal dimensions as the regularization term. We also search the weights of regularization terms in {1, 0.1, 0.01} and report the best results. As for the PINNacle, we report the results of canonical PINN following their paper [12] and experiment with other base models by only replacing the model.vPINN We follow the code base in PINNacle, and implement it to the first three benchmarks. The test functions are set as Legendre polynomials and the test function number is set as 5. The number of points used to compute the integral within domain Ω is set as 10, and the number of grids is set differently for each subtask, with values of {4, 8, 16, 32} for PINNacle and the same to other baselines for the first three benchmarks.Other baselines In Table 3 of the main text, we also experiment with the loss-reweighting method NTK [47] and data-resampling method RAR [51] . For NTK, we follow their official code and recalculate the neural tangent kernel to update loss weights every 10 iterations. And the kernel size is set as 300. As for RAR, we use the residual-based adaptive refinement with distribution algorithm.In this section, we provide more results as a supplement to the main text, including additional hyperparameter analysis, new experiments and more showcases.As we stated in Algorithm 1, we adopt the gradient variance of past T 0 iterations to approximate the sampling error defined in Theorem 3.9. In our experiments, we choose T 0 from {5, 10}, which can [36] and PINNsFormer [58] on different benchmarks. We change T 0 in {1, 5, 10, 15, 20, 25, 30} and record the rMSE.achieve consistently good and stable performance among different benchmarks and PDEs. To analyze the effect of this hyperparameter, we further add experiments with different choices in Figure 8 .As shown in Figure 8 , we can find that under all the choices in {1, 5, 10, 15, 20, 25, 30}, RoPINN performs better than the vanilla PINN. Specifically, in both 1D-Reaction and 1D-Wave (Figure 8(a-b )), the model performs quite stable under different choices of T 0 . As for Convection in Figure 8(c-d ), the influence of T 0 is relatively significant in PINN. This may caused by the deficiency of PINN in solving Convection, where all the PINN-based experiments fail to generate an accurate solution for Convection (rMSE>0.5, Table 2 ). If we adopt a more powerful base model, such as PINNsFormer [58] , this sensitivity will be alleviated. Also, it is worth noticing that, even though in Convection, RoPINN surpasses the vanilla PINN under all hyperparameter settings of T 0 .Besides, we can observe that the model performance slightly decreases when we set T 0 with a relatively large value. This may come from the difference between parameters θ t and θ t+29 , which will make the gradient variance approximation less reliable (Eq. ( 58) in the Theorem 3.11 proof). As we stated in the implementations (Appendix C.3), RoPINN can also be applied to tasks with data loss. Here we also include an inverse problem in PINNacle to testify to the performance of RoPINN in this case, which requires the model to reconstruct the diffusion coefficients of the Poisson equation from observations on 2500 uniform grids with additional Gaussian noise.As presented in Table 5 , in this task, RoPINN can also boost the performance of PINN with over 10% in the rMSE metric and outperform the other baselines (gPINN and vPINN) that cannot bring improvements. Note that in this experiment, we failed to reproduce the performance of vPINN reported by PINNacle. Thus, we report the results of vPINN by directly running the official code in PINNacle.Considering the limited resources, we repeat all the experiments on the first three typical benchmarks and our method on the PINNacle three times and other experiments one time. The official paper of PINNacle has provided the standard deviations for PINN, gPINN and vPINN on all benchmarks.We summarize the standard deviations of PINN in Table 6 . As for other base models, the standard deviations of FLS, QRes and KAN are within 0.005 on 1D-Wave and Convection, and within 0.001 for 1D-Reaction. PINNsFormer's standard deviations are smaller than 0.001 for all three benchmarks.Table 6 : Standard deviations for canonical PINN on three typical benchmarks. The confidence for RoPINN achieving the best performance is over 99% in all three benchmarks.rMSE±Standard Deviations 1D-Reaction 1D-Wave Convection PINN [36] 0.981±5e-4 0.335±1e-3 0.840±5e-4 +gPINN [55] 0.978±3e-4 0.399±3e-3 0.935±3e-3 +vPINN [18] 0.982±3e-3 0.173±1e-3 0.743±2e-3 +RoPINN (Ours) 0.095±8e-4 0.064±1e-3 0.720±2e-3 Burges 1d-C 1.1e-2 3.3e-2 3.7e-1 5.1e-1 4.0e-2(-272.2%) 3.5e-1(-952.3%) 9.1e-3(15.8%) 1.4e-2(56.6%) 2d-C 4.5e-1 5.2e-1 4.9e-1(-8.7%) 5.4e-1(-3.8%) 6.6e-1(-46.4%) 6.4e-1(-23.0%) 4.3e-1(4.0%) 4.9e-1(5.3%) Poisson 2d-C 7.5e-1 6.8e-1 7.7e-1(-3.0%) 7.0e-1(-4.0%) 4.6e-1(38.6%) 4.9e-1(27.4%) 4.1e-1(44.9%) 6.6e-1(2.5%) 2d-CG 5.4e-1 6.6e-1 7.4e-1(-37.1%) 7.9e-1(-21.0%) 2.4e-1(54.7%) 2.9e-1(56.4%) 4.1e-1(24.1%) 6.0e-1(8.1%) 3d-CG 4.2e-1 5.0e-1 4.3e-1(-4.2%) 5.2e-1(-3.5%) 8.0e-1(-91.4%) 7.4e-1(-46.7%) 4.7e-1(-11.9%) 4.6e-1(8.7%) 2d-MS 7.8e-1 6.4e-1 6.7e-1(13.2%) 6.2e-1(2.3%) 9.6e-1(-23.6%) 9.7e-1(-52.1%) 7.7e-1(0.4%) 6.4e-1(0.3%) QRes [3] Burges 1d-C 5.8e-3 2.0e-2 3.6e-1 5.1e-1 2.0e-2(-241.0%) 8.1e-2(-309.3%) 5.7e-3(1.4%) 1.8e-2(6.8%) 2d-C 3.2e-1 4.8e-1 4.9e-1(-53.1%) 5.4e-1(-11.1%) 1.1e+0(-256.9%)1.3e+0(-161.2%) 3.2e-1(0.3%) 4.7e-1(1.8%) Poisson 2d-C 2.9e-1 6.8e-1 7.6e-1(-159.4%) 7.2e-1(-5.9%) 3.2e-1(-7.6%) 3.0e-1(55.6%) 2.9e-1(0.6%) 7.0e-1(-3.7%) 2d-CG 3.1e-1 7.4e-1 7.3e-1(-136.9%) 7.8e-1(-6.0%) 6.4e-1(-108.9%) 7.3e-1(0.4%) 2.9e-1(4.0%) 7.2e-1(2.5%) 3d-CG 9.0e-2 5.8e-1 5.9e-1(-558.8%) 6.4e-1(-10.7%) 8.0e-1(-790.4%) 7.4e-1(-28.4%) 8.9e-2(1.2%) 5.6e-1(3.5%) 2d-MS 1.7e+0 7.9e-1 5.5e-1(67.2%) 5.0e-1(36.5%) 9.7e-1(41.8%) 9.8e-1(-24.4%) 1.9e+0(-13.1%) 9.1e-1(-15.0%) Heat 2d-VC 2.0e-1 1.3e+0 7.2e+0 5.9e+0(-363.0%) 3.3e-1(-63.9%) Burges 1d-C 9.0e-3 1.4e-2 3.3e-1 4.8e-1 6.5e-2(-620.9%) 3.0e-1 9.0e-3(0.2%) 1.3e-2(8.6%) 2d-C 4.4e-1 4.9e-1 4.9e-1(-13.0%) 5.4e-1(-10.1%) 1.3e+0(-195.0%)1.3e+0(-171.4%) 4.3e-1(0.4%) 4.9e-1(0.2%) Poisson 2d-C 6.8e-1 6.7e-1 7.3e-1(-6.6%) 6.7e-1(-1.2%) 9.0e-1(-30.9%) 9.4e-1(-40.5%) 7.2e-1(-5.4%) 6.3e-1(5.2%) 2d-CG 6.3e-1 7.0e-1 7.8e-1(-22.7%) 8.1e-1(-16.2%) 6.0e-1(5.6%) 7.0e-1(-0.6%) 5.9e-1(6.8%) 6.9e-1(1.2%) 3d-CG 4.2e-1 5.0e-1 4.5e-1(-7.8%) 4.6e-1(9.5%) 8.1e-1(-95.5%) 7.5e-1(-48.8%) 4.1e-1(1.3%) 5.1e-1(-0.5%) 2d-MS 8.9e-1 7.5e-1 5.1e-1(43.5%) 4.7e-1(36.4%) 9.0e-1(-0.0%) 9.4e-1(-25.7%) 8.1e-1(9.5%) 6.7e-1(10.0%) Poisson 2d-C 7.3e-1 6.8e-1 7.4e-1(-1.4%) 6.8e-1(-0.1%) 1.6e-1(78.2%) 1.5e-1(78.4%) 6.3e-1(13.4%) 6.8e-1(-0.1%) 2d-CG 7.7e-1 8.1e-1 8.1e-1(-5.5%) 8.6e-1(-6.9%) 6.5e-1(14.9%) 7.3e-1(9.8%) 6.5e-1(15.8%) 6.9e-1(14.0%) 3d-CG 7.5e-1 1.4e+0 2.0e+0(-163.6%) 1.8e+0(-25.5%) NaN NaN 7.2e-1(4.1%) 6.8e-1(51.8%) 2d-MS 9.5e-1 9.8e-1 1.0e+0(-5.0%) 1.0e+0(-2.7%) 9.9e-1(-4.6%) 1.0e+0(-2.3%) 9.2e-1(3.4%) 9.6e-1(2.1%) Heat 2d-VC 1.9e+1 1.5e+1 3.3e+0(82.5%) 2.7e+0(81.8%) 8.9e-1(95.3%) 9.1e-1(94.0%) 6.0e+0(68.7%) 4.7e+0(68.4%) 2d-MS 1.4e+0 1.1e+0 2.0e+0(-39.8%) 1.4e+0(-28.4%) 8.3e-1(42.8%) 7.5e-1(31.4%) 7.4e-1(49.1%) 6.9e-1(36.8%) 2d-CG 5.0e-1 5.3e-1 8.2e-1(-63.1%) 7.3e-1(-37.2%) 8.1e-1(-60.8%) 9.0e-1(-69.7%) 4.9e-1(1.7%) 5.0e-1(4.9%)2d-C 5.0e-1 4.1e-1 8.1e-1(-60.6%) 6.7e-1(-65.8%) This section will discuss some related works as a supplement to Section 2. We will first discuss some PINN research and then we will also clarify some looking similar but completely distinct topics.PINN optimizers As we mentioned in the second paragraph of the introduction, many previous works focus on developing efficient and effective deep-model optimizers for PINNs [55, 37] , which may help the optimization process tackle the ill-conditioned Hessian matrix or naturally balance multiple loss terms [53] . As we formalized in Algorithm 1, RoPINN is not restricted to a certain optimizer. The researchers can easily replace the Adam [21] or L-BFGS [27] with other advanced optimizers. Since we mainly focus on the objective function, these works are orthogonal to us.Numerical differentiation for objective functions In addition to the regularization or variationalbased methods, some researchers attempt to replace the automatic differentiation with numerical approximations [39, 9] , which can tackle the expensive computation cost caused by calculating high-order derivatives. However, this paradigm does not attempt to change the objection function definition, just focuses on the calculation of point optimization PINN loss, which is distinct from our proposed region optimization paradigm.Sampling-based methods Strategies for sampling collocation points perform an important role in training PINN models [46] . Previous sampling-based methods mainly focus on accumulating collocation points to high-residual areas [51, 23] or considering the temporal causality [45] , whose theoretical analyses are usually based on the quadrature theorem [32] . Distinct from these methods, RoPINN is motivated by the optimization deficiency of PINN models and can be seamlessly integrated with sampling-based methods with significant promotion (Table 3 ), indicating that RoPINN works orthogonally to sampling methods. Besides, the theoretical analysis of RoPINN also starts from the optimization perspective, which reveals that one key advancement of RoPINN is a better balance between optimization error and generalization error (Theorem 3.12).In addition, RoPINN is also distinct from data augmentation or adversarial training techniques in the following aspects: (1) Theorem difference: although our proposed practical algorithm is based on Monte Carlo sampling in a region, the underlying theoretical support and insights are a region-based objective function (Theorem 3.5). ( 2) Implementation difference: In our algorithm, not only is the input changed, but the objective is also correspondingly changed. Thus, this paper is foundationally different from augmentation and adversarial training in that the ground truth label is fixed. Our design is tailored to the physics-informed loss function of PINNs, where we can accurately calculate the equation residual at any point within the input domain.This paper presents region optimization as a new PINN training paradigm and provides both theoretical analysis and practical algorithms, supported by extensive experiments. However, there are still several limitations. In the theoretical analysis, we assume that the canonical loss function is L-Lipschitz-β-smooth, which may not be guaranteed in practice. Besides, RoPINN involves several hyperparameters, such as initial region size r, and number of past iterations T 0 . Although we have studied the sensitivity w.r.t. them in Figures 2 and Appendix D.1 and demonstrate that they are easy to tune in most cases, we still need to adjust them for better performance in practice.According to our experiments and theorems, we provide some recipes for hyperparameter tuning in the following, which may be helpful to the usage of RoPINN:• As shown in Figure 2 , region size r will be progressively adjusted by RoPINN. Setting r in [10 -6 , 10 -4 ] can work well. According to Theorem 3.12, the choice of r should balance optimization and generalization, which may be inherently decided by the PDE smoothness. • As analyzed in Figures 3 and 4 , sampling 1-30 points can gain consistent promotion but will linearly increase the computation costs. Following our default setting (sampling one point) can already achieve a competitive performance in a wide range of PDEs. • As presented in Figure 8 , number of past iterations T 0 is easy to tune in [1, 20] . Setting T 0 ∈ {5, 10} can be a good choice, which has been widely verified in our paper. • Some hyperparameter tuning tools, such as Weights and Bias (Wandb1 ), may mitigate this limitation to some extent, which has already been used in previous related work [37] .In this paper, we develop a new region optimization training paradigm for PINNs and provide both theorem analyses and practical algorithms. This new perspective may inspire the subsequent research of PINNs, especially rethinking the canonical objective function. In addition, our proposed RoPINN shows favorable efficiency and generalizes well in different base models and PDEs, which can be used to boost the precision of PINNs and generally benefit the downstream tasks, such as physics phenomenon simulation, biological property analysis, etc. Since we purely focus on the training algorithm of PINNs, there are no potential negative social impacts or ethical risks.NeurIPS Paper Checklist 2, 8 and 9 ) and detailed discussions of related work (Section 2 and Appendix F) are provided to verify the effectiveness of our design and clarify the scope of our paper. Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper.https://github.com/wandb/wandb 


Decision: The submission has received mixed ratings after extensive post-rebuttal discussion between the authors and the referees. One of the referees recommends rejection based mainly on presentation aspects. The other three referees recommend acceptance, and they also rate the presentation as sufficiently clear. The authors are advised to take into account the detailed reviews as well as the clarifications provided in their rebuttal while revising the submission.

Reviews: review: paper_summary: paper_summary: The paper proposes a novel optimization method for training physics-informed neural networks (PINNs): Region optimization, which extends a regular pointwise optimization of PINNs to neighborhood regions, named RoPINNs. The paper provides theoretical analysis explaining the decrease of generalization error with RoPINNs and high-order PDE constraints satisfaction. Then the paper presents a practical algorithm to enable the RoPINNs paradigm by introducing Monte-Carlo approximation of region integral and a region calibration scheme. Lastly, the paper assesses the performance on several well-known benchmark problems and showed the improved performance over the considered baselines.
strength_weakness: strength_weakness: - Although shown to be very effective in several benchmark problems, the paper does not seem to provide general guidelines on how to set some important hyper-parameters such as initial region size and the number of sampling points. (While acknowledging that the authors indicate this as one of the limitations,) it would be great to see some experts’ guidelines.

- If the authors could provide some analysis with regards to computational wall time, that would provide more complete pictures on how the proposed method performs. For example, it would be information to see a figure depicting a scatter plot showing computational wall time versus final rMSE type information, where a point in the plot corresponds to a different hyper-parameter setting (that is, the number of sample points). 

- [minor] there is a typo in the second paragraph of Section 4.2: line 289 Figure 2 => Figure 3.
questions: questions: - Eq (5) seems to suggest that region optimization is applied to the boundary condition as well as L = L_bc + L_ic + L_pde. Is this the correct understanding or is it a typo?

- The proposed optimization method seems to benefit significantly in the case of 1D reaction case while the benefit in 1D Wave or 1D convection cases are not as significant as that of 1D reaction. That is, rMSE for example in Table 2 of 1D reaction achieves an order (or orders) of magnitude improvement over the second best performing methods. Do the authors have some explanation on why?
limitations: limitations: - some more discussions on practical guidelines would be needed for users who want to utilize this method in different applications 

- some additional experiments (regarding computational wall time) would be needed to provide a complete picture of the proposed method.
paper_summary: paper_summary: The authors developed a region optimized PINN to improve the prediction accuracy compared to the scatter-point based PINN.
strength_weakness: strength_weakness: The practical application scope is limited.
questions: questions: (1)  It is suggested to add some descriptions of training difficulty factors for the canonical PINN on 1D-Reaction in Section 4.2. 
(2)  Can the proposed method find a good number of sampling points well balancing the computational cost and convergence speed in Figure 3.
(3)  The motivation of using Monte Carlo approximation should be elaborated. Why don’t the authors choose some other more advanced methods to provide better accuracy.
(4)  The authors should add more details about the possible practical applications with the canonical loss function of L-Lipschitz-β-smooth.
limitations: limitations: Some initial guess methods can be developed to efficiently determine the preferable region size and sample number.
paper_summary: paper_summary: This paper extends the optimization process of PINNs from isolated points to their continuous neighborhood regions, which can theoretically decrease the generalization error, especially for hidden high-order constraints of PDEs. A practical training algorithm, Region Optimized PINN (RoPINN), is seamlessly derived from this new paradigm, which is implemented by a straightforward but effective Monte Carlo sampling method. By calibrating the sampling process into trust regions, RoPINN finely balances sampling efficiency and generalization error. Experimentally, RoPINN consistently boosts the performance of diverse PINNs on a wide range of PDEs without extra backpropagation or gradient calculation.
strength_weakness: strength_weakness: 1. It is better to include the main proof idea of theoretical results in the main text.
2. Although generalization error bound is provided, an intuitive explanation of the reason behind the success of region optimization is desirable. For example, when sampling one point in each region, why is the total loss decreased compared with point optimization?
questions: questions: 1. Which results in section 4 are for comparisons with the losses with high-order terms and variational methods?
2. How tight is the generalization error bound?
limitations: limitations: The case of sampling more than one points in each region is not discussed.
paper_summary: paper_summary: The preprint proposes to replace the collocation based PINN loss by a sum of local continuous integrals over regions around the collocation points. These continuous integrals are then again discretized using Monte Carlo integration with a single quadrature pint. The authors furthermore propose to adapt the region size during training using gradient statistics.
strength_weakness: strength_weakness: The introduction of continuous integrals over regions around the collocation points that subsequentially are discretized by Monte Carlo integration again seems tautological. After all, the loss function in PINNs is already a Monte Carlo discretization of a continuous integral (over the whole computational domain). Furthermore, the analysis that the authors present for the modified loss in equation (5) should not be carried out with the continuous integral over the regions $\Omega_r$ but with its Monte Carlo approximation. Otherwise, the comparison to the discretized PINN loss is unfair.
questions: questions: I struggle to see why the proposed method should work theoretically. I acknowledge the adaptive nature of the regions for the sampling but struggle to see how this might help to accumulate integration points in regions of, e.g., high residual. A visualization that this, or something along these lines that explains why the proposed method works well, happens would be helpful. Furthermore, I am not convinced that the theoretical analysis presented is meaningful, as it analyzes the integrals over the region as continuous objects. Please comment or clarify misconception.
limitations: limitations: See questions

Comments: title: More Discussion about Presentation
Summary of rebuttal to Reviewer mhfV
Thanks for your reply and acknowledging our "good" contribution
Thanks for your response and positive support
Thanks for your response and acknowledge our rebuttal
More discussions about sampling-based methods
Thanks for your response and more explanations about how RoPINN works
comment: Many thanks for your time in reviewing our paper and new opinion about writing.

For a convenient overview, we respectfully list some key reasons why we think our current writing from an "optimization" view is more suitable than sampling:

**(1) All the theoretical analyses and proof of our paper are based on an optimization-based view. Previous sampling papers are under different theoretical frameworks from us.** 

  Here is a comparison of two different perspectives. We believe that the previous theoretical framework of sampling methods is distinct from ours, which has also been clarified in $\underline{\text{Clarify misconception of original rebuttal}}$.

|  Theoretical framework            | Our Paper                            | Previous Sampling Paper (e.g. RAR)                           |
| ------------ | ------------------------------------ | ------------------------------------------------------------ |
| Key Question | Will the PINN model be well-trained? | Will the sampling well approximate integral?                 |
| Theorem      | **Optimization-based theorem**   | **Quadrature-based theorems** [Siddhartha Mishra et al., IMAJNA 2023] |

**(2) The optimization-based view can better present how our method works.** In our current paper, all the analyses are based on an optimization perspective. Since we do not change the collocation points (center of region) during training and only adjust the region size, the visualization of point distribution change is not suitable for us. Here is a comparison.

|   Analysis                        | Our paper                                          | Previous Sampling Paper (e.g. RAR) |
| ------------------------- | -------------------------------------------------- | ---------------------------------- |
| Training Curve            | $\underline{\text{Figures 2,3 of main text}}$      | N/A                                |
| Gradient std Statistics   | $\underline{{\text{Q2-(2) of original rebuttal}}}$ | N/A                                |
| Point distribution change | N/A                                                | The most widely used analysis tool     |

**(3) Our contribution is orthogonal to previous sampling papers with significant promotion**. In $\underline{\text{Table 3 of main text}}$, we present that RoPINN can be combined with RAR to further boost their performance by a great margin, which means RoPINN is from a new idea (may be more foundational). Here are part of the results.

|  rMSE                        | 1D-Reaction | 1D-Wave | Convection |
| ------------------------ | ----------- | ------- | ---------- |
| PINN vanilla             | 0.981       | 0.335   | 0.840      |
| +RAR                     | 0.981       | 0.126   | 0.771      |
| +RAR+RoPINN              | 0.080       | 0.030   | 0.695      |
| **Promotion w.r.t. RAR** | **92%**     | **76%** | **10%**    |

Since the discussion period will end in a few hours, we do hope that the reviewer can reconsider our presentation.

Sincerely thanks for your time and active discussion.
We sincerely thank your dedication in reviewing our paper. 

Since this is the last day of the review-author discussion period, we summarize the key points of our rebuttal as follows for a convenient overview:

- **Clarify the concern about practical application scope:** (1) We highlighted that "Theoretical assumption will not affect the practicability of algorithms". (2) The extensive experiment results (5 base models on 19 PDEs) have well supported RoPINN's practicability. (3) Our theoretical assumption is widely used in other papers.
- **Add new experiments about more sampling points:** We highlight that only simpling 1 point is already a good choice and also add experiments on more points to provide an overview of RoPINN.
- **Add comparison w.r.t. advanced quadrature methods:** We show that the advanced methods need more points to achieve a better performance, which defeats our purpose in "no extra gradient calculation".
- **Give practical guidance to hyperparameters** based on analysis experiments in the main text.

**Sincerely thank you for ranking our paper with "excellent" soundness, "good" presentation and "good" contribution**. We do hope our rebuttal can fully resolve your questions to your satisfaction. If so, with the greatest respect, we also hope you can kindly consider raising the score correspondingly.

Many thanks for your time and looking forward to your response.
We sincerely thank you for your response, appreciating our idea and acknowledging our "good" contribution.

> "I disagree that you propose an optimization method."

**We do NOT attempt to say that RoPINN is an optimization method**. What we want to do is to design a new training paradigm considering the optimization process of PINN models, **which is why we name RoPINN as "region optimized PINN" not "region optimizer" in the title**.

Thus, we respectfully point out that there may exist some misunderstandings in this comment.

> "You propose an adaptive sampling method with different criteria than the residual."

As per your request, we explain RoPINN in a simpling-based method style in the previous reply to provide a better understanding for you.

However, we have to emphasize that **the contribution of region optimization is more than different criteria than residual**. As we stated in the $\underline{\text{Clarify misconception in previous rebuttal}}$, both theorem and algorithm are proposed based on the idea of considering the training process.

**Thus, our key contribution is taking the optimization process into consideration and building a complete theorem and practical algorithm to implement this idea.** We believe that this new theoretical framework can be a good assistant to previous methods, which is also acknowledged by you.

> "I think the manuscript would benefit from a different presentation that puts the sampling viewpoint as the core novelty. I am still against publication in the current form but now mainly because of the presentation."

We do appreciate the reviewer's suggestion in the paper presentation, which gives us another perspective to think about our paper. Many thanks for the detailed descriptions in your question, which help us a lot to understand your concern.

**However, we have to clarify that the sampling-based explanation is just following your request. We do hope that you can rethink our paper from our original "optimization-based" insights (see $\underline{\text{Q1-(1) in original rebuttal}}$).** In this context, current writing about region optimization theorem and practical algorithm is more fluent and smoother than sampling-based ideas. We believe that this disagreement is only about preference in academic perspectives.

Besides, we also want to highlight our experiment results, **consistently improving 5 base models on 19 different PDEs with over 50% gain in most cases**. These significant results well support that our idea may be more foundational than sampling.

## Our promise in revision and respectful request to reconsider our paper presentation

As we stated in the previous response, we will polish our presentation by incorporating the **newly proved theorem, experimental statistics and discussion w.r.t. previous sampling methods** into our paper, which will clearly present our paper from an "optimization-based" view. We will consider your valuable suggestions and do our best to make our paper easy to understand.

With the greatest respect, we also request you to reconsider your score. We believe that **both sample-based and optimization-based viewpoints can give a reasonable explanation to our design**, where the optimization-based writing is also $\underline{\text{accepted by the other two reviewers}}$. Thus, **we kindly say that your concern in the presentation may be caused by two different academic perspectives.**

Greatest thanks for your time in reviewing our paper and discussion.
Thank you for the detailed answer. In fact, I can understand your point more clearly now and like the idea of informing the sampling based on gradient statistics instead of or in combination with residual-based methods.

I still find the explanation in your write-up hard to digest, introducing regions around sampling points which are treated like integrals and then again discretized by one collocation point was confusing to me (and also to another reviewer) and I apologize it took some time before I could understand the underlying idea. 

Furthermore, I disagree that you propose an optimization method. You propose an adaptive sampling method with different criteria than the residual.

To conclude, I will raise my score but I think the manuscript would benefit from a different presentation that puts the sampling viewpoint as the core novelty. I am still against publication in the current form but now mainly because of the presentation.
Thanks for your prompt response and for acknowledging our new experiments in computation-accuracy balancing. Your positive support and valuable suggestions help us a lot in revising our paper. We will include all the new results in the revised paper.
We would like to sincerely thank you for your valuable response and suggestions.

(1) As for "A larger gradient estimation error will lead to an unstable training process", as we stated in the response to Reviewer zNBp, this is just an intuitive understanding. Here are more details.

Please recall the optimization algorithm SGD and GD. SGD generally faces more optimization fluctuations, thereby converging more slowly than GD. Similarly, if the gradient among successive iterations has a low consistency, the optimization direction of PINN models will face fluctuations, which will make it harder to arrive at a convergence point. If the gradient among successive iterations is consistent, the optimization direction is confirmed, thereby converging faster.

(2) As for "considering the smoothness property of loss landscape", we think the loss landscape will affect the gradient consistency, where the gradient consistency is a dependent variable. Suppose that the target PDE is easy to solve, RoPINN can learn to explore larger region size $r$ in pursuing a better balance between generalization bound and optimization (gradient consistency).

Following your suggestion, we will add more discussion about the intuition behind RoPINN (including the newly proved theorem for practical algorithm) and the difference w.r.t. previous sampling methods in the revised paper.

Sincerely thanks for your time in reviewing our paper.
Thank you for sharing the new results. They provide useful insight into the tradeoff between additional computation and accuracy. I remain positive about the paper and will maintain the current score.
Thank the authors for their replies. Aftering reading the authors' responses to me and other reviewers, especially why introducing gradient estimation error is benefit and the comparision between RoPINN and other sampling method,  I got the rough idea of why region optimization works, though I did not check the proof details.  "A larger gradient estimation error (corresponding to a lower consistency of gradients in successive iterations) will lead to an unstable training process", this seems to be related with the loss landscape of PINNs, how does this gradient consistency promote optimization convergence considering the smoothness property of loss landscape? 

I am satisfied with the answers to my other questions. I would like to keep my score, and suggest the authors to explain in detail the intuition behind RoPINN  and the  difference w.r.t. previous sampling methods in the updated submission.
**We believe that accumulating points in high-residual areas is NOT the only principle to design sampling-based methods.**

Actually, although the previous sampling methods (e.g. RAR) can make the collocation point accumulate to high-residual areas, **they may also face optimization difficulty due to too many hard-to-optimize points.** For example, if one point is high-residual but extremely hard to optimize, the model optimization may be over-attracted by this point and be misguided.

RoPINN is based on a distinct idea, which **considers the model optimization process**. Specifically, the calibrated region size $r$ can finely balance optimization difficulty and "exploration" to unseen areas (i.e. larger sampling regions).

|                           | Design                       | Pros                                                         | Cons                                                        |
| - | - | - | - |
| Previous methods e.g. RAR | Sampling high-residual areas | Make the model solve hard areas better       |  optimization difficulty |
| RoPINN                    | Adaptive sampling region size      | Balance optimization difficulty and "exploration" to unseen areas | Without specific optimization to hard areas        |

The Cons of RoPINN can be solved by integrating RoPINN and RAR ($\underline{\text{Table 3 of main text}}$). Note that we are not saying that RoPINN is better than RAR. They contribute orthogonally.

**We do hope that the reviewer can leave the sampling high-residual design behind (which is distinct from our design) and think about our paper based on our proposed "optimization-based" theorems.**

Many thanks for your dedication to our paper.
We sincerely thank the reviewer's response and detailed descriptions of the question.

We kindly remind you that our previous rebuttal has included the requested theorem for practical algorithm in $\underline{\text{Q1-(1)}}$, which shows that introducing "region" is not tautological and ensures a fair theoretical comparison w.r.t. discretized PINN loss. Thus, we think the theoretical value of introducing "region" is clear. 

Next, we will explain how RoPINN helps the PDE-Solving in practice.

> Could the authors explain to me how the proposed method "informs" the sampling procedure?

Firstly, we list the comparison between RoPINN and a sampling method RAR to clarify "what RoPINN informs the sampling" and "why it works". More details are included in the following.

|                       | Information to sampling                                      | Intuitive understanding of benefits                          |
| - | - | - |
| RAR [46 of our paper] | PDE residual                                                 | Make the algorithm aware of which area is hard to solve.     |
| RoPINN                | Gradient estimation error (gradient variance of successive iterations) | Make the algorithm aware if the current sampling range can bring a "good" optimization, which refers to a relatively stable training for convergence and is sufficient to make the model "explore" new areas for generalization. |

**(1) RoPINN "informs" the sampling procedure with "gradient estimation error" (gradient variance of successive iterations).**

As shown in $\underline{\text{Algorithm 1 of main text}}$, RoPINN sets the sampling region size $r$ of each iteration as $r/\sigma_t$, where $\sigma_t$ denotes the gradient variance of successive iterations. In $\underline{\text{Theorem 3.9, 3.11 of main text}}$, we have also proved that $\sigma_t$ can be used to approximate gradient estimation error. Thus, RoPINN's design can introduce the gradient estimation error (gradient variance) into the sampling procedure.

Why introducing gradient estimation error is beneficial? Here are the explanations.

- **Theoretical analysis:** As shown in $\underline{\text{Q1-(1) of previous rebuttal}}$, we have proved that in practice, the generalization bound is also affected by gradient estimation error. The detailed explanations for the effect of $r$ are included in the previous rebuttal. One key point is that a proper sampling region size $r$ can better balance $(1-|\Omega_r|/|\Omega|)L$ and $\mathcal{E}_{r, \mathrm{grad}}$, bringing a better performance. What RoPINN does is change $r$ adaptively.
- **Intuitive understanding:** **The gradient estimation error of MC can be used to represent the “consistency” of optimization direction within a region** ($\underline{\text{Theorem 3.9 of main text}}$). A larger gradient estimation error (corresponding to a lower consistency of gradients in successive iterations) will lead to an unstable training process, which may overwhelm the model or even fail to converge. On the other side, too stable training is also insufficient to make the model "explore" new areas and damage the generalization. RoPINN can achieve a balanced result, which is supported by our new statistics in $\underline{\text{Q2-(2) of previous rebuttal}}$.

The above "intuitive understanding" explanation has been partially discussed in $\underline{\text{Lines 207-210 of main text}}$. In the revised paper, we will rephrase this paragraph for a more detailed explanation.

**(2) RoPINN can be combined with the "sampling techniques" that you mentioned.**

Actually, in $\underline{\text{Q2 of previous rebuttal}}$, we have pointed out that the contribution of RoPINN is orthogonal to previous sampling methods. The key idea of RoPINN is to extend optimization targets from collocation points to their regions. Thus, RoPINN can be combined with your mentioned "sampling techniques". Here are part of the results. Please see $\underline{\text{Table 3 of main text}}$ for full results.

| rMSE            | 1D-Reaction | 1D-Wave   | Convection |
| - | - | - | - |
| PINN            | 0.981       | 0.335     | 0.840      |
| +RAR            | 0.981       | 0.126     | 0.771      |
| +RoPINN         | 0.095       | 0.064     | 0.720      |
| **+RAR+RoPINN** | **0.080**   | **0.030** | **0.695**  |

To better explain how our model works, we will rephrase our paper by:

- Incorporating the requested theorem in $\underline{\text{Q1-(1) of previous rebuttal}}$ in the main text to illustrate the balance between optimization and generalization.
- Incorporating the new statistics in $\underline{\text{Q2-(2) of previous rebuttal}}$ to show the balanced results of RoPINN.
- Explaining that region optimization is a general framework for Fixed or Global sampling point optimization.
- Adding a new section to explain the relation and difference w.r.t. previous sampling methods.

We hope these new results can resolve your concerns and we are happy to answer any further questions.
I thank the authors for their replies. I still struggle to see the value of the proposed region optimization.

Could the authors explain to me how the proposed method "informs" the sampling procedure? All sampling techniques for PINN type loss functions I am aware of reallocate points to regions that are challenging to learn for the network. This is informed by quantities that relate to how well the PDE is currently solved -- typically the PDE residual. I cannot really see that in your proposed method. Does this somehow implicitly happen? 

I acknowledge that Figure 2 in the main text illustrates that the regions change in size during the training process, but this does not answer my question.

I remain critical of the work.