Newton Informed Neural Operator for Computing Multiple Solutions of Nonlinear Partials Differential Equations

Authors: wenrui  hao
xinliang  liu
yahong  yang
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-13 18:20:24

Abstract: Solving nonlinear partial differential equations (PDEs) with multiple solutions is essential in various fields, including physics, biology, and engineering. However, traditional numerical methods, such as finite element and finite difference methods, often face challenges when dealing with nonlinear solvers, particularly in the presence of multiple solutions. These methods can become computationally expensive, especially when relying on solvers like Newton's method, which may struggle with ill-posedness near bifurcation points.
In this paper, we propose a novel approach, the Newton Informed Neural Operator, which learns the Newton solver for nonlinear PDEs. Our method integrates traditional numerical techniques with the Newton nonlinear solver, efficiently learning the nonlinear mapping at each iteration. This approach allows us to compute multiple solutions in a single learning process while requiring fewer supervised data points than existing neural network methods.

Main paper: 
 introduction: 
 Neural networks have been extensively applied to solve partial differential equations (PDEs) in various fields, including biology, physics, and materials science [21, 9] . While much of the existing work focuses on PDEs with a unique solution, nonlinear PDEs with multiple solutions pose a significant challenge [39, 10] but are widely encountered in applications such as [1, 3, 31, 36, 14, 13] . In this paper, we aim to solve the following nonlinear PDEs with multiple solutions:EQUATIONHere, Ω is the domain of equation, f (u) is a nonlinear function in R, u : R d → R and L is a second-order elliptic operator given by Lu = -d i,j=1 a ij (x)u xixj + d i=1 b i (x)u xi + c(x)u, for given coefficient functions a ij , b i , c(i, j = 1, . . . , d) with n i,j=1 a ij (x)ξ i ξ j ≥ λ|ξ| 2 , for a constant λ ≥ 0.Various neural network methods have been developed to tackle partial differential equations (PDEs), including PINN [33] , the Deep Ritz method [45] , DeepONet [29] , FNO [25] , MgNO [16] , and OL-DeepONet [26] . These methods can be broadly categorized into two types: function learning and operator learning approaches. In function learning, the goal is to directly learn the solution. However, these methods often encounter the limitation of only being able to learn one solution in each learning process. Furthermore, the problem becomes ill-posed when there are multiple solutions. On the other hand, operator learning aims to approximate the map between parameter functions in PDEs and the unique solution. This approach cannot address the issue of multiple solutions or find them in a single training session. We will discuss this in more detail in the next section.In this paper, we present a novel neural network approach for solving nonlinear PDEs with multiple solutions. Our method is grounded in operator learning, allowing us to capture multiple solutions within a single training process, thus overcoming the limitations of function learning methods in neural networks. Moreover, we enhance our network architecture by incorporating traditional Newton methods [35, 1] , as discussed in the next section. This integration ensures that the problem of operator learning becomes well-defined, as Newton's methods provide well-defined locally, thereby ensuring a robust operator. This approach addresses the challenges associated with directly applying operator networks to such problems. Additionally, we leverage Newton information during training, enabling our method to perform effectively even with limited supervised data points. We introduce our network as the Newton Informed Neural Operator. To clarify, we do not design a specific neural structure for the neural operator. The Newton information is not incorporated into the architecture of the neural network but rather into the training process. Specifically, the Newton method is incorporated into the loss function, as detailed in Section 3.3.As mentioned earlier, our approach combines the classical Newton method, which translates nonlinear PDEs into an iterative process involving solving linear functions at each iteration. One key advantage of our method is that, once the operator is effectively learned, there is no need to solve the linear equation at every iteration. This significantly reduces computation time, especially in complex systems encountered in fields such as material science, biology, and chemistry. Furthermore, once the Newton Informed Neural Operator is well-trained, it can be applied to compute new solutions with appropriate initial guesses, even those not present in the training data. Details of this capability are demonstrated in the numerical example of the Gray-Scott model. Overall, the Newton Informed Neural Operator efficiently solves nonlinear PDEs with multiple solutions by learning the Newton nonlinear solver. It addresses the time-consuming nature of traditional nonlinear solvers and requires fewer supervised data points compared to existing neural network methods. Additionally, it saves time by eliminating the need for repeatedly solving nonlinear systems, as is required in traditional Newton methods. Once the neural operator is learned, it can also compute new solutions beyond those provided in the supervised data.The following paper is organized as follows: In the next section (Section 2), we will delve into nonlinear PDEs with multiple solutions and discuss related works on solving PDEs using neural network methods. In Section 3, we will review the classical Newton method for solving PDEs and introduce the Newton Informed Neural Operator, which combines neural operators with the Newton method to address nonlinear PDEs with multiple solutions. In this section, we will also analyze the approximation and generalization errors of the Newton Informed Neural Operator. Finally, in Section 4, we present the numerical results of our neural networks for solving nonlinear PDEs. The first example demonstrates that the Newton Informed Neural Operator requires minimal data for training, the second example shows that the speed of the Newton Informed Neural Operator is significantly faster than the traditional Newton method, and the last example highlights that the Newton Informed Neural Operator can discover new solutions not present in the supervised data. 
backgrounds and relative works 2.1 nonlinear pdes with multiple solutions: 
 Significant mathematical models depicting natural phenomena in biology, physics, and materials science are rooted in nonlinear partial differential equations (PDEs) [5] . These models, characterized by their inherent nonlinearity, present complex multi-solution challenges. Illustrative examples include string theory in physics, reaction-diffusion systems in chemistry, and pattern formation in biology [20, 12] . However, experimental techniques like synchrotronic and laser methods can only observe a subset of these multiple solutions. Thus, there is an urgent need to develop computational methods to unravel these nonlinear models, offering deeper insights into the underlying physics and biology [17] . Consequently, efficient numerical techniques for identifying these solutions are pivotal in understanding these intricate systems. Despite recent advancements in numerical methods for solving nonlinear PDEs, significant computational challenges persist for large-scale systems. Specifically, the computational costs of employing Newton and Newton-like approaches are often prohibitive for the large-scale systems encountered in real-world applications. In response to these challenges [15, 19] , we propose an operator learning approach based on Newton's method to efficiently solve nonlinear PDEs.Indeed, there are numerous approaches to solving partial differential equations (PDEs) using neural networks. Broadly speaking, these methods can be categorized into two main types: function learning and operator learning.In function learning, neural networks are used to directly approximate the solutions to PDEs. Function learning approaches aim to directly learn the solution function itself. On the other hand, in operator learning, the focus is on learning the operator that maps input parameters to the solution of the PDE. Instead of directly approximating the solution function, the neural network learns the underlying operator that governs the behavior of the system.In function learning, a commonly employed method for addressing this problem involves the use of Physics-Informed Neural Network (PINN)-based learning approaches, as introduced by Raissi et al. [33] , and Deep Ritz Methods [45] . However, in these methods, the task becomes particularly challenging due to the ill-posed nature of the problem arising from multiple solutions. Despite employing various initial data and training methods, attaining high accuracy in solution learning remains a complex endeavor. Even when a high-accuracy solution is achieved, each learning process typically results in the discovery of only one solution. The specific solution learned by the neural network is heavily influenced by the initial conditions and training methods employed. However, discerning the relationships between these factors and the learned solution remains a daunting task. In [19] , the authors introduce HomPINNs for learning multiple solutions to PDEs, where the number of solutions that can be learned depends on the choice of "start functions." However, if the "start functions" are not appropriately selected, HomPINNs may fail to capture all solutions. In this paper, we present an operator learning approach combined with Newton's method to train the nonlinear solver. While this approach is not specifically designed for computing multiple solutions, it can be employed to compute them if suitable initial guesses are provided.Operator learning methods Various approaches have been developed for operator learning to solve PDEs, including DeepONet [29] , which integrates physical information [7, 26] , as well as techniques like FNO [25] inspired by spectral methods, and MgNO [16] , HANO [27] , and WNO [24] based on multilevel methods, and transformer-based neural operators [4, 8] . These methods focus on approximating the operator between the parameters and the solutions. Firstly, they require the solutions of PDEs to be unique; otherwise, the operator is not well-defined. Secondly, they focus on the relationship between the parameter functions and the solution, rather than the initial data and multiple solutions.3 Newton Informed Neural Operator 
review of newton methods to solve nonlinear partial differential equations: 
 To tackle this problem Eq. ( 1), we employ Newton's method by linearizing the equation.For the Newton method applied to an operator, if we aim to find the solution of F(u) = 0, the iteration can be written as:F ′ (u n ) u n+1 = F ′ (u n ) u n -F (u n ) ⇐⇒ F ′ (u n ) δu = -F (u n ) ,where δu = u n+1 -u n .In this context, F ′ (u)v is the (Fréchet) derivative of the operator, which is a linear operator to v, defined as follows: To find F ′ (u) in X , for any v ∈ X ,lim |v|→0 |F(u + v) -F(u) -F ′ (u)v| |v| = 0,where | • | denotes the norm in X .For solving Eq. ( 1), given any initial guess u 0 (x), for i = 1, 2, . . . , M , in the i-th iteration of Newton's method, we have ũ(x) = u + δu(x) by solvingEQUATIONwhich is based on the fact that the (Fréchet) derivative of 2) has a unique solution, then by solving Eq. ( 2) and repeating the process M times, we will obtain one of the solutions of the nonlinear equation (1) . Denoting the mapping for u and δu, the solution of Eq. ( 2) with parameter u, as G(u) := δu, we know thatL -f (•) at u is L -f ′ (u). If Eq. (lim n→∞ (G + Id) (n) (u 0 ) = u * ,where u * is one of the solutions of Eq. (1). For different initial conditions, this process will converge to different solutions of Eq. (1), making this method suitable for finding multiple solutions. Furthermore, the Newton method is well-posed, meaning that each initial condition u 0 will converge to a single solution of Eq. (1) under appropriate assumptions (see Assumption 1). This approach helps to address the ill-posedness encountered when using PINNs directly to solve Eq. (1). However, repeatedly solving Eq. ( 1) can be computationally expensive, especially in high-dimensional cases or when a large number of discrete points are involved. In this paper, we tackle these challenges by employing neural networks.In this section, we introduce the structure of the neural operator to approximate the operator locally in the Newton methods from Eq.(2), i.e., δu := G(u), where δu is the solution of Eq.( 2), which depends on u. If we can learn the operator G(u) well using the neural operator O(u; θ), then for an initial function u 0 , assume the n-th iteration will approximate one solution, i.e., (G + Id)(n) (u 0 ) ≈ u * . Thus, (O + Id) (n) (u 0 ) ≈ (G + Id) (n) (u 0 ) ≈ u * .For another initial condition, we can evaluate our neural operator and find the solution directly.Then we discuss how to train such an operator. To begin, we define the following shallow neural operators with p neurons for operators from X to Y asEQUATIONwhereW i ∈ L(X , Y), B i ∈ Y, A i ∈ L(Y, Y), and θ denote all the parameters in {W i , A i , B i } p i=1 . Here, L(X , Y) denotes the set of all bounded (continuous) linear operators between X and Y, and σ : →R defines the nonlinear point-wise activation.In this paper, we will use shallow DeepONet [29, 22] to approximate the Newton operator. To provide a more precise description, in the shallow neural network, W i represents an interpolation of operators. With proper and reasonable assumptions, we can present the following theorem to ensure that DeepONet can effectively approximate the Newton method operator. The proof will be provided in the appendix. Furthermore, MgNO is replaced by W as a multigrid operator [38] , and FNO is some kind of kernel operator; our analysis can be generalized to such cases.Before the proof, we need to establish some assumptions regarding the input space X ⊂ H 2 (Ω) of the operator and f (u) in Eq. (1) . The definition of the Sobolev space can be found in Appendix B.1. Assumption 1. (i): For any u ∈ X , we have that the linear equation (L -f ′ (u))δu = -L + f (u) is well-posed for solving δu.(ii): There exists a constant F such that ∥f (x)∥ W 2,∞ (R) ≤ F . , where α k are the coefficients in the expansion of u with respect to the basis {b n }. The canonical projection operator P n is a linear bounded operator on X . According to the properties of Schauder bases, these projections P n are uniformly bounded by some constant C. Furthermore, we assume, for any u ∈ X , ϵ > 0, there exists a n such that ∥u -P n u∥ H 2 (Ω) ≤ ϵ, for all u ∈ X .More discussion about the assumption is shown in the appendix. The sketch of the proof is illustrated in Fig. 1 .Figure 1 : The sketch of proof for Theorem 1.Theorem 1. Suppose X = Y ⊂ H 2 (Ω) and Assumption 1 holds. Then, there exists a neural network O(u; θ) ∈ Ξ p defined asEQUATIONsuch thatEQUATION)where σ is a smooth non-polynomial activation function, n is shown in Assumption 1 and contained in W i , C 1 is a constant independent of m, ϵ, and p, C 2 is a constant depended on p, n and F (see in Assumption 1) is the scale of the P in Assumption 1. And ϵ depends on n. n and ϵ are defined in Assumption 1.The approximation results of DeepONet in Sobolev training can be found in [40] .The Mean Square Error loss function is defined as:E S (θ) := 1 M u • M x Mu j=1 Mx k=1 |G (u j ) (x k ) -O (u j ; θ) (x k )| 2 (6)where u 1 , u 2 , . . . , u Mu ∼ µ are independently and identically distributed (i.i.d) samples in X , and x 1 , x 2 , . . . , x Mx are uniformly i.i.d samples in Ω.However, using only the Mean Squared Error loss function is not sufficient for training to learn the Newton method, especially since in most cases, we do not have enough data {u j , G (u j )} Mu j=1 . Furthermore, there are situations where we do not know how many solutions exist for the nonlinear equation (1) . If the data is sparse around one of the solutions, it becomes impossible to effectively learn the Newton method around that solution.Given that E S (θ) can be viewed as the finite data formula of E Sc (θ), whereE Sc (θ) = lim Mu,Mx→∞ E S (θ).The smallness of E Sc can be inferred from Theorem 1. To understand the gap between E Sc (θ) and E S (θ), we can rely on the following theorem. Before the proof, we need some assumptions about the data in E S (θ): Assumption 2. (i) Boundedness: For any neural network with bounded parameters, characterized by a bound B and dimension d θ , there exists a function Ψ :L 2 (Ω) → [0, ∞) such that |G(u)(x)| ⩽ Ψ(u), sup θ∈[-B,B] d θ |O(u; θ)(x)| ⩽ Ψ(u), sup θ∈[-B,B] d θ |LO(u; θ)(x)| ⩽ Ψ(u)for all u ∈ X , x ∈ Ω, and there exist constants C, κ > 0, such thatEQUATION(ii) Lipschitz continuity: There exists a function Φ :EQUATIONfor all u ∈ X , x ∈ Ω, and Φ(u) ⩽ C 1 + ∥u∥ H 2 (Ω) κ , for the same constants C, κ > 0 as in Eq. ( 7).(iii) Finite measure: There exists α > 0, such thatH 2 (Ω) e α∥u∥ 2 H 2 (Ω) dµ(u) < ∞.Theorem 2. If Assumption 2 holds, then the generalization error is bounded bysup θ∈[-B,B] d θ |E(E S (θ) -E Sc (θ))| ⩽ C 1 √ M u 1 + Cd θ log(CB M u ) 2κ+1/2 + d θ √ log M x √ M x ,where C is a constant independent of B, d θ , M x , and M u . The parameter κ is specified in (7) . Here, B represents the bound of parameters, and d θ is the number of parameters.The proof of Theorem 2 is presented in Appendix B.3. Remark 1. Assumption 2 is easily satisfied if we consider X as the local function set around the solution, which is typically the case in Newton's methods. This aligns with our approach and the working region in the approximation part (see Remark 3). The errorsup θ∈[-B,B] d θ |E(E S (θ) - E Sc (θ))| suggests that the network can perform well based on the loss function E S (θ). The reasoning is as follows: letθ S = arg min θ∈[-B,B] d θ E S (θ) and θ Sc = arg min θ∈[-B,B] d θ E Sc (θ).We aim for EE Sc (θ S ) to be small, which can be written as:EE Sc (θ S ) ≤ E Sc (θ Sc ) + E(E S (θ S ) -E Sc (θ S )) ≤ E Sc (θ Sc ) + sup θ∈[-B,B] d θ |E(E S (θ) -E Sc (θ))|,where E Sc (θ Sc ) is small, as demonstrated by Theorem 1 when B is sufficiently large.As we have mentioned, relying solely on the MSE loss function can require a significant amount of data to achieve the task. However, obtaining enough data can be challenging, especially when the equation is complex and the dimension of the input space is large. Hence, we need to consider another loss function to aid learning, which is the physical information loss function [33, 7, 19, 24] , referred to here as the Network loss function.The Newton loss function is defined as:E N (θ) := 1 N u • N x Nu j=1 Nx k=1 |(L -f ′ (u j (x k )))O (u j ; θ) (x k ) + Lu j (x k ) + f (u j (x k ))| 2 (9)where u 1 , u 2 , . . . , u Nu ∼ ν are independently and identically distributed (i.i.d) samples in X , and x 1 , x 2 , . . . , x Nx are uniformly i.i.d samples in Ω.Given that E N (θ) can be viewed as the finite data formula of E N c (θ), whereE N c (θ) = lim Nu,Nx→∞ E S (θ).To understand the gap between E N c (θ) and E N (θ), we can rely on the following theorem: Corollary 1. If Assumption 2 holds, then the generalization error is bounded bysup θ∈[-B,B] d θ |E(E N (θ) -E N c (θ))| ⩽ C 1 √ N u 1 + Cd θ log(CB N u ) 2κ+1/2 + d θ √ log N x √ N x ,where C is a constant independent of B, d θ , N x , and N u . The parameter κ is specified in (7) . Here, B represents the bound of parameters, and d θ is the number of parameters.The proof of Corollary 1 is similar to that of Theorem 2; therefore, it will be omitted from the paper. Remark 2. If we only utilize E S (θ) as our loss function, as demonstrated in Theorem 2, we require both M u and M x to be large, posing a significant challenge when dealing with complex nonlinear equations. Obtaining sufficient data becomes a critical issue in such cases. In this paper, we integrate Newton's information into the loss function, defining it as follows:EQUATION) where E N (θ) represents the cost associated with the unsupervised learning data. If we lack sufficient data for E S (θ), we can adjust the parameters by selecting a small λ and increasing N x and N u . This strategy enables effective learning even when data for E S (θ) is limited. We refer to this neural operator, which incorporates Newton information, as the Newton Informed Neural Operator.In the following experiment, we will use the neural operator established in Eq. ( 3) and the loss function in Eq. (10) to learn one step of the Newton method locally, i.e., the map between the input u and the solution δu in eq. ( 2). If we have a large dataset, we can choose a large λ in E(θ) (10); if we have a small dataset, we will use a small λ to ensure the generalization of the operator is minimized. After learning one step of the Newton method using the operator neural networks, we can easily and quickly obtain the solution by the initial condition of the nonlinear PDEs (1) and find new solutions not present in the datasets. 
experimental settings: 
 We introduce two distinct training methodologies. The first approach employs exclusively supervised data, leveraging the Mean Squared Error Loss (6) as the primary optimization criterion. The second method combines both supervised and unsupervised learning paradigms, utilizing a hybrid loss function 10 that integrates Mean Squared Error Loss (6) for small proportion of data with ground truth (supervised training dataset) and with Newton's loss (9) for large proportion of data without ground truth (unsupervised training dataset). We call the two methods method 1 and method 2. The approaches are designed to simulate a practical scenario with limited data availability, facilitating a comparison between these training strategies to evaluate their efficacy in small supervised data regimes. We chose the same configuration of the neural operator (DeepONet) which is aligned with our theoretical analysis. One can find the detailed experimental settings and the datasets for each example below in Appendix A.We consider 2D convex problem L(u)-f (u) = 0, where L(u) := -∆u, f (u) : -u 2 +sin 5π(x+y) and u = 0 on ∂Ω. We investigate the training dynamics and testing performance of neural operator (DeepONet) trained with two methods, focusing on Mean Squared Error (MSE) and Newton's loss functions. For method 1, we use 500 supervised data samples (with ground truth), while for method 2, we use 5000 unsupervised data samples (only with the initial state) along with supervised data samples, employing the regularized loss function as defined in Equation 10 with λ = 0.01. Please refer A.1.1 for the dataset generation of convex case. For the detailed experimental settings, refer to Appendix A. We consider a 2D Non-convex problem,EQUATIONwhere Ω = (0, 1) × (0, 1) [3] . In this case, L(u) := -∆(u), f (u) := u 2 -s sin(πx) sin(πy) and it has multiple solutions (see Figure 3 for its solutions).In the experiment, we let one of the multiple ground truth solutions rarely touched in the supervised training dataset such that the neural operator trained via method 1 will saturate in terms of test error because it relies on the ground truth to recover all the patterns for multiple solution cases (as shown by the curves in Figure 3 ). On the other hand, the model trained via method 2 is less affected by the limited supervised data since the utilization of Newton's loss. One can refer to Appendix A for the detailed experiment setting.Efficiency This case study highlights the superior efficiency of our neural operator-based method as a surrogate model for Newton's method. Both methods parallelize operations to solve 500/5000 Newton linear systems simultaneously, each with distinct initial states. The key advantage of the neural operator lies in its ability to batch the computation of these independent systems efficiently. By efficiently batching and sampling a wide variety of initial states, the neural operator improves the likelihood of discovering to multiple solutions, particularly in nonlinear PDEs with complex solution landscapes. Consequently, while Newton's method alone does not inherently guarantee finding multiple solutions, the combination of rapid computation and extensive initial condition sampling enhances the chances of identifying multiple solutions.For a fair comparison, the classical Newton solver was also parallelized using CUDA on a GPU. However, the neural operator naturally handles large batch sizes during inference, allowing it to process all systems in one go. One can find the detailed description of the experiments in A.6. The table demonstrates the significant efficiency gain achieved by batching the computation of independent Newton systems with distinct initial states using the neural operator. For NINO, solving 5000 independent Newton linear systems scales up minimally compared to solving 500 systems, while the classical solver experiences a tenfold increase in computation time. This efficient batching is crucial for improving performance, particularly in complex nonlinear systems like the Gray-Scott model, where solving numerous systems simultaneously is essential for effective pattern discovery.The Gray-Scott model [31, 11] describes the reaction and diffusion of two chemical species, A and S, governed by the following equations:∂A ∂t = D A ∆A -SA 2 + (µ + ρ)A, ∂S ∂t = D S ∆S + SA 2 -ρ(1 -S),where D A and D S are the diffusion coefficients, and µ and ρ are rate constants.Newton's Method for Steady-State Solutions Newton's method is employed to find steady-state solutions ( ∂A ∂t = 0 and ∂S ∂t = 0) by solving the nonlinear system:EQUATION)The Gray-Scott model is highly sensitive to initial conditions, where even minor perturbations can lead to vastly different emergent patterns. Please refer to Figure 5 for some examples of the patterns. This sensitivity reflects the model's complex, non-linear dynamics that can evolve into a multitude of possible steady states based on the initial setup. Consequently, training a neural operator to map initial conditions directly to their respective steady states presents significant challenges. Such a model must learn from a vast functional space, capturing the underlying dynamics that dictate the transition from any given initial state to its final pattern. This complexity and diversity of potential outcomes is the inherent difficulty in training neural operators effectively for systems as complex as the Gray-Scott model. One can refer to A.1.2 for a detailed discussion on the Gray-Scott model.We employ a Neural Operator as a substitute for the Newton solver in the Gray-Scott model, which recurrently maps the initial state to the steady state. In subfigure (a), we use a ring-like pattern as the initial state to test our learned neural operator. This pattern does not appear in the supervised training dataset and lacks corresponding ground truth data. Instead, it is present only in the unsupervised data (Newton's loss), i.e., some data in Newton's loss will converge to this specific pattern. Despite this, our neural operator, trained using Newton's loss, can effectively approximate the mapping of the initial solution to its correct steady state. we further test our neural operator, utilizing it as a surrogate for Newton's method to address nonlinear problems with an initial state drawn from the test dataset. The curve shows the average convergence rate of ∥u -u i ∥ across the test dataset, where u i represents the prediction at the i-th step by the neural operator. In subfigure (c), we compare the Training Loss (Rescaled Newton's Loss) and Absolute L2 Test Error. The magnitudes are not directly comparable as they represent different metrics; however, the trends are consistent, indicating that the inclusion of unsupervised data and training with Newton's loss contributes to improved model performance. 
conclusion: 
 In this paper, we develop neural operators to learn the Newton's solver related to nonlinear PDEs (Eq. ( 1)) with multiple solutions. To speed up the computation of multiple solutions for nonlinear PDEs, we combine neural operator learning with classical Newton methods, resulting in the Newtoninformed neural operator. We provide a theoretical analysis of our neural operator, demonstrating that it can effectively learn the Newton operator, reduce the number of required supervised data, and learn solutions not present in the supervised learning data due to the addition of the Newton loss (9) in the loss function. Our experiments are consistent with our theoretical analysis, showcasing the advantages of our network as mentioned earlier.Reproducibility Statement Code Availability: The code used in our experiments can be accessed via https://github.com/xlliu2017/learn_newton and also the supplementary material. Datasets can be downloaded via URLs in the repository. This encompasses all scripts, functions, and auxiliary files necessary to reproduce our results. Configuration Transparency: All configurations, including hyperparameters, model architectures, and optimization settings, are explicitly provided in the Appendix. Function and Jacobian The function F (u) might typically be defined as:F (u) = -∆u + u 2The Jacobian J(u), for the given function F (u), involves the derivative of F with respect to u, which includes the Laplace operator and the derivative of the nonlinear term:J(u) = -∆ + 2 • u.The dataset are generated by sampling the initial state u 0 ∼ N (0, ∆ -3 ) and then calculate the convergent sequence {u 0 , u 1 , ..., u n } by Newton's method. Each convergent sequence {u 0 , u 1 , ..., u n } is one data entry in the dataset.The analysis of function and Jacobian for the non-convex problem (case 2) is similar to the convex problem except that its Jacobian J(u) = ∆ -2u such that Newton's system is not positive definite. 
experimental setting/details: 
 Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes]Justification: We provide the detailed configurations of the experiments in the Appendix.Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [Yes]Justification: We report statistical significance for key findings in our study. Since we do not claim our method is superior to other models, our experiments are serving to justify our theoretical claims, therefore, we selectively provide statistical results where they are most relevant.Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.• The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. Justification: We report the compute resources for benchmarking our method.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines? 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [NA] Justification: The research is mostly theoretical. We do not see obvious social impact.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • 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). 
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: No obvious such risks.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: [NA] Justification: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. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 
new assets: 
 Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?Answer: [NA] Justification:Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. • At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file. 
crowdsourcing and research with human subjects: 
 Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?Answer: [NA] Justification:Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector. 
institutional review board (irb) approvals or equivalent for research with human subjects: 
 Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?Answer: [NA] Justification:Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.• We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: Jacobian Matrix The Jacobian matrix J of the system is crucial for applying Newton's method:J = J AA J AS J SA J SSwith components:J AA = -D A ∆ + diag(-2SA + µ + ρ), J AS = diag(-A 2 ), J SA = diag(2SA), J SS = -D S ∆ + diag(A 2 + ρ).The numerical simulation of the Gray-Scott model was configured with the following parameters:• Grid Size: The simulation grid is square with N = 63 points on each side, leading to a total of N 2 grid points. This resolution was chosen to balance computational efficiency with spatial resolution sufficient to capture detailed patterns. The spacing between each grid point, h, is computed as h = 1.0 N -1 . This ensures that the domain is normalized to a unit square, which simplifies the analysis and scaling of diffusion rates.• Diffusion Coefficients: The diffusion coefficients for species A and S are set to D A = 2.5 × 10 -4 and D S = 5.0 × 10 -4 , respectively. These values determine the rate at which each species diffuses through the spatial domain.• Reaction Rates: The reaction rate µ and feed rate ρ are crucial parameters that govern the dynamics of the system. For this simulation, µ is set to 0.065 and ρ to 0.04, influencing the production and removal rates of the chemical species.Simulations The simulation utilizes a finite difference method for spatial discretization and Newton's method to solve the steady-state given the initial state. The algorithm is detailed in A.6.Figure 6 : Three examples depicting the evolution from the initial state to the steady state via Newton's method.Here is how we generate the supervised data samples:1.Step 1: We use a classical numerical solver to obtain single (multiple) solutions of nonlinear PDEs. For example, in Case 2, there exist four solutions u 1 , u 2 , u 3 , u 4 . We have one solution for Case 1 and 10 solutions for the Case 3: Gray-Scott model.Step 2: The supervised dataset is generated by sampling a perturbation around the solution u i ∼ N (0, (-∆) -3 ) on the chosen true solution u i . We then set u i 0 = u i p + u i and calculate the convergent sequence u i 0 , u i 1 , . . . , u i n using Newton's method, which follows the formula:u i k+1 = u i k -J f (u i k ) -1 f (u i k ). Each convergent sequence u i 0 , u i 1 , .. . , u i n constitutes one supervised data entry in the dataset. In this case, we consider the initial conditions as perturbed, with all perturbations applied around the true solution. A comparison between the traditional method and our proposed method is summarized in Table 1 .For the unsupervised data samples, we only sample the perturbed initial states.Discrete Newton's Loss In solving partial differential equations (PDEs) numerically on a regular grid, the Laplace operator and other differential terms can be efficiently computed using convolution. Here, we detail the method for calculating J(u)δu -F (u) where J(u) is the Jacobian matrix, δu is the Newton step, and F (u) is the function defining the PDE.Discretization Consider a discrete representation of a function u on a N × N grid. The function u and its perturbation δu are represented as matrices:u, δu ∈ R N ×NThe function F (u), which involves both linear and nonlinear terms, is similarly represented asF (u) ∈ R N ×N .Laplace Operator Regarding the representing the J(u) with N ×N grid function u, the discretized Laplace operator using a finite difference method can be expressed as a convolution:-∆u = 0 -1 0 -1 4 -1 0 -1 0 * uThis convolution computes the result of the Laplace operator applied to the grid function u. The boundary conditions can be further incorporated into the convolution with different padding modes. Dirichlet boundary condition corresponds to zeros padding while Neumann boundary condition corresponds to replicate padding.A variant of DeepONet is used in our Newton-informed neural operator. In the DeepONet, we introduce a hybrid architecture that combines convolutional layers with a static trunk basis, optimized for grid-based data inputs common in computational applications like computational biology and materials science.The branch network is designed to effectively downsample and process the spatial features through a series of convolutional layers:• A Conv2D layer with 128 filters (7x7, stride 2) initiates the feature extraction, reducing the input dimensionality while capturing coarse spatial features.• This is followed by additional Conv2D layers (128 filters, 5x5 kernel, stride 2 and subsequently 3x3 with padding, 1x1) which further refine and compact the feature representation.• The convolutional output is flattened and processed through two fully connected layers (256 units then down to branch features), using GELU activation.Trunk Network The trunk utilizes a static basis represented by the tensor V, incorporated as a non-trainable component: The tensor V is precomputed, using Proper Orthogonal Decomposition (POD) as in [28] , and is dimensionally compatible with the output of the branch network.Forward Pass During the forward computation, the branch network outputs are projected onto the trunk's static basis via matrix multiplication, resulting in a feature matrix that is reshaped into the grid dimensionality for output.The following table 2 Below we summarize the key configurations and parameters employed in the training for three cases:• Case 1: For method 1, we use 500 supervised data samples (with ground truth) while for method 2, we use 5000 unsupervised data samples (only with the initial state) and 500 supervised data samples.• Case 2: For method 1, we use 5000 supervised data while for method 2, we use 5000 unsupervised data samples and 5000 supervised data samples.• Case 3 (Gray-Scott model): We only perform method 2, with 10000 supervised data samples and 50000 unsupervised data samples.• Optimizer: Adam, with a learning rate of 1 × 10 -4 and a weight decay of 1 × 10 -6 .• Training Epochs: The model was trained over 1000 epochs to ensure convergence and we use Batch Size: 50.These settings underscore our commitment to precision and detailed examination of neural operator efficiency in computational tasks. Our architecture and optimization choices are particularly tailored to explore and exploit the capabilities of neural networks in processing complex systems simulations.A.6 Benchmarking Newton's Method and neural operator based method Experimental Setup The benchmark study was conducted to evaluate the performance of a GPUaccelerated implementation of Newton's method, designed to solve systems of linear equations derived from discretizing partial differential equations. The implementation utilized CuPy with CUDA to leverage the parallel processing capabilities of the NVIDIA A100 GPU. The hardware comprises an Intel Cascade Lake 2.5 GHz CPU, and an NVIDIA A100 GPU.The performance was assessed in terms of total execution time, which includes the setup of matrices and vectors, computation on the GPU, and synchronization of CUDA streams. Both methods leverage the parallel processing capabilities of the GPU. Specifically, the Newton solver explicitly uses 10 streams and CuPy with CUDA to parallelize the computation and fully utilize the GPU parallel processing capabilities, aiming to optimize execution efficiency. In contrast, the neural operator method is inherently parallelized, taking full advantage of the GPU architecture without the explicit use of multiple streams as indicated in the table. The computational times of both methods were evaluated under a common hardware configuration.Software Environment: Ubuntu 20.04 LTS. Python Version: 3.8. CUDA Version: 11.4.Newton's method was implemented to solve the Laplacian equation over a discretized domain.Multiple system solutions were computed in parallel using CUDA streams. The key parameters of the experiment are as follows: Data Type (dtype): Single precision floating point (float32). Number of Streams: 10 CUDA streams to process data in parallel. Number of Repeated Calculations:The Newton method was executed multiple times for 500/5000 Newton linear systems, respectively, distributed evenly across the streams. Function to Solve Systems: The CuPy's spsolve was used for solving the sparse matrix systems. The following algorithm A.6 summarizes the procedure to benchmark the time used for solving Newton's system for 5000 different initial states.Algorithm 1 Solve Newton's Systems on GPU 1: procedure SOLVESYSTEMSGPU(A, u, rhs_f , N ) ▷ Precompute RHS and diagonal for all systems 2:rhs_list ← rhs_f + u 2 -A × u 3: diag_list ← -2 × u ▷ Initialize solution storage 4:delta_u ← initialize zero matrix with shape of u T ▷ Solve each system 5:for i = 0 to num_sys -1 do 6: rhs ← rhs_list[i] 7: diag ← diag_list[i] 8:J ← A + diagonal matrix with diag on the main diagonal ▷ Solve the linear system 9:delta_u[i] ← spsolve(J, rhs) W n,p (Ω) := f ∈ L p (Ω) : D α f ∈ L p (Ω) for all α ∈ N d with |α| ≤ n with a norm ∥f ∥ W n,p (Ω) :=   0≤|α|≤n ∥D α f ∥ p L p (Ω)   1/p if p < ∞, and ∥f ∥ W n,∞ (Ω) := max 0≤|α|≤n ∥D α f ∥ L ∞ (Ω) . Furthermore, for f = (f 1 , . . . , f d ), f ∈ W 1,∞ (Ω, R d ) if and only if f i ∈ W 1,∞ (Ω) for each i = 1, 2, . . . , d and ∥f ∥ W 1,∞ (Ω,R d ) := max i=1,...,d {∥f i ∥ W 1,∞ (Ω) }. When p = 2, denote W n,2 (Ω) as H n (Ω) for n ∈ N + . Proposition 1 ([30]). Suppose σ is a is a continuous non-polynomial function and K is a compact in R d , then there are positive integers p, constants w k , ζ k for k = 1, . . . , p and bounded linear functionalsEQUATION)Proposition 2 ( [32, 44] ). Suppose σ is a continuous non-polynomial function and Ω is a compact subset of R d . For any Lipschitz-continuous function f , there exists a shallow neural network such thatEQUATIONwhere C depends on the Lipschitz constant but is independent of m. Lemma 1 ( [22] ). The ϵ-covering number of [-B, B] d , K(ϵ), satisfiesK(ϵ) ⩽ CB ϵ d ,for some constant C > 0, independent of ϵ, B, and d.Step 5: Now we estimateIn this subsection, we present the proof of Theorem 1, which describes the approximation ability of DeepONet.Proof of Theorem 1.Step 1: Firstly, we need to verify that the target operator G(u) is well-defined.Due to Assumption 1 (i), we know that for u ∈ X ⊂ H 2 (Ω), Equation (2) will have unique solutions. This means that G(u) is a well-defined operator for the input space u ∈ X .Step 2: Secondly, we aim to verify that G(u) is a Lipschitz-continuous operator in H 2 (Ω) for u ∈ X .Consider the following:EQUATIONwhere δu(x) is the solution of Eq.( 2) for the input u, and δ v u(x) is the solution of Eq.( 2) for the input u + v. Denote δ v u(x) -δu(x) =: ϵ(x).Therefore, we have:EQUATIONSince u and v are in H2 and ∂Ω is in C 2 (Assumption 1 (iii)), according to [6, Theorem 4 in Section 6.3], there exist constants C 2 and C such that:EQUATIONThe last inequality is due to the boundedness of f ′ (ξ 2 ) + f ′′ (ξ 1 )δu (Assumption 1 (ii)).Step 3: In the approximation, we first reduce the operator learning to functional learning.When the input function u(x) belongs to X ⊂ H 2 , the output function δu also belongs to H 2 , provided that ∂Ω is of class C 2 . The function G(u) = δu can be approximated by a two-layer network architected by the activation function σ(x), which is not a polynomial, in the following form by Proposition 1 [30] (given in Subsection 16):EQUATIONwherew k ∈ R d , ζ k ∈ R for k = 1, . . . , p, c k is a continuous functional, and C 1 is a constant independent of the parameters. Denote ϕ k (u) = c k [G(u)],which is a Lipschitz-continuous functional from H 2 (Ω) to R, which is due to G is a Lipschitz-continuous operator and c k is a linear functional. The remaining task in approximation is to approximate these functionals by neural networks.Step 4: In this step, we reduce the functional learning to function learning by applying the operator P as in Assumption 1 (iv).Based on ϕ k (u) being a Lipschitz-continuous functional in H 2 (Ω), we have|ϕ k (u) -ϕ k (Pu)| ≤ L k ∥u -Pu∥ H 2 (Ω) ≤ L k ϵ,where L k is the Lipschitz constant of ϕ k (u) for u ∈ X .Furthermore, since Pu is an n-dimensional term, i.e., it can be denoted by the n-dimensional vector Pu ∈ R n , we can rewrite ϕ k (Pu) as ψ k ( Pu), where ψ k : R n → R for l = 1, . . . , p. Furthermore, ψ k is a Lipschitz-continuous function since ϕ k is Lipschitz-continuous and P is a continuous linear operator.Step 5: In this step, we will approximate ψ k using shallow neural networks.Due to Proposition 2, we have that there is a shallow neural network such thatEQUATIONwhere Above all, we have that there is a neural network in Ξ p such thata ⊺ k ∈ R m , M k ∈ R m×n ,EQUATIONwhere C 1 is a constant independent of m, ϵ, and p, C 2 is a constant depended on p.Remark 3. We want to emphasize the reasonableness of our assumptions. For condition (i), we are essentially restricting our approximation efforts to local regions. This limitation is necessary because attempting to approximate the neural operator across the entire domain could lead to issues, particularly in cases where multiple solutions exist. Consider a scenario where the input function u lies between two distinct solutions. Even a small perturbation of u could result in the system converging to a completely different solution. Condition (i) ensures that Equation (2) has a unique solution, allowing us to focus our approximation efforts within localized domains.Conditions (ii) and (iii) serve to regularize the problem and ensure its tractability. These conditions are indeed straightforward to fulfill, contributing to the feasibility of the overall approach.For the embedding operator P in (iv), there are a lot of choices in DeepONet, such as finite element methods like Argyris elements [2] or embedding methods in [23, 18] . We will discuss more in the appendix. We omit the detailed discussion in the paper. Furthermore, for the differential neural network, this embedding may be different; for example, we can use Fourier expansion [43] or multigrid methods [16] to achieve this task.Here, we discuss more about the embedding operator P. One approach is to use the Argyris element [2] . This method involves considering the 21 degrees of freedom shown in Fig. 7 . In this figure, each• denotes evaluation at a point, the inner circle represents an evaluation of the gradient at the center, and the outer circle denotes evaluation of the three-second derivatives at the center. The arrows indicate the evaluation of the normal derivatives at the three midpoints.Another alternative approach to discretizing the input space is to use the bi-cubic Hermite finite element method [23, 18] .The proof of Theorem 2 is inspired by that in [22] .Proof of Theorem 2.Step 1: To begin with, we introduce a new term called the middle term, denoted as E Sm (θ), defined as follows:E Sm (θ) := 1 M u Mu j=1 Ω |G(u j )(x) -O(u j ; θ)(x)| 2 dx,This term represents the limit case of E S (θ) as the number of samples in the domain of the output space tends to infinity (M x → ∞).Then the error can be divided into two parts:EQUATIONStep 2: For |E(E Sm (θ) -E Sc (θ)|, this is the classical generalization error analysis, and the result can be obtained from [34, 42, 41] . We omit the details of this part, which can be expressed asEQUATIONwhere C is independent of the number of parameters d θ and the sample size M x . In the following steps, we are going to estimate |E(E S (θ)-E Sm (θ))|, which is the error that comes from the sampling of the input space of the operator.Step 3: DenoteS M θ := 1 M M j=1 Ω |G(u j )(x) -O(u j ; θ)(x)| 2 dx.We first estimate the gap between S M θ and S M θ ′ for any bounded parameters θ, θ ′ . Due to Assumption 2 (i) and (ii), we have thatEQUATIONStep 4: Based on Step 3, we are going to estimateE sup θ∈[-B,B] d θ S M θ -ES M θ p 1 pby covering the number of spaces.Set {θ 1 , . . . , θ K } is a ε-covering of [-B, B] d θ i.e. for any θ in [-B, B] d θ , there exists j with ∥θ -θ j ∥ ℓ∞ ⩽ ϵ. Then we haveEQUATION)For 8ϵE [|ΨΦ| p ] 1/p , it can be approximate by8ϵE [|ΨΦ| p ] 1/p ⩽ 8ϵE |Ψ| 2p 1/2p E |Φ| 2p 1/2p = 8ϵ∥Ψ∥ L 2p ∥Φ∥ L 2p . For E max j=1,...,K S M θj -E S M θj p 1/p, by applied the result in [37, 22] , we knowE max j=1,...,K S M θj -E S M θj p 1/p ≤ 16K 1/p √ p∥Ψ∥ 2 L 2p √ M .Step 5: Now we estimate |E(E S (θ) -E Sm (θ)|.Due to Assumption 2 and directly calculation, we have that∥Ψ∥ L 2p , ∥Φ∥ L 2p ⩽ C(1 + γκp) κ ,for constants C, γ > 0, depending only the measure µ and the constant C appearing in the upper bound (7) . For example,EQUATION)Based on Lemma 1, we have thatE sup θ∈[-B,B] d θ S Mu θ -E S Mu θ p 1/p ⩽ 16C 2 (1 + γκp) 2κ ϵ + CB ϵ d θ /p √ p √ M u ,for some constants C, γ > 0, independent of κ, µ, B, d θ , N, ϵ > 0 and p ⩾ 2. We now chooseϵ = 1 √ Mu , so that ϵ + CB ϵ d θ /p √ p √ M u = 1 √ M u 1 + (CB M u ) d θ /p √ p .Next, let p = d θ log(CB √ M u ). Then,(CB M u ) d θ /p √ p = exp log(CB √ M u )d θ p √ p = e d θ log(CB M u ),and thus we conclude thatϵ + CB ϵ d θ /p √ p √ M u ⩽ 1 √ M u1 + e d θ log(CB M u ). .On the other hand, we have(1 + γκp) 2κ = 1 + γκd θ log(CB M u ) 2κ .Increasing the constant C > 0, if necessary, we can further estimate1 + γκd θ log(CB M u ) 2κ 1 + e d θ log(CB M u ). ⩽ C 1 + d θ log(CB M u ) 2κ+1/2, where C > 0 depends on κ, γ, µ and the constant appearing in (7) , but is independent of d θ , B and N . We can express this dependence in the form C = C(µ, Ψ, Φ) > 0, as the constants κ and γ depend on the Gaussian tail of µ and the upper bound on Ψ, Φ.Therefore,EQUATION)Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material?Answer: [Yes] Justification: We provide the link to the datasets.Guidelines:• The answer NA means that paper does not include experiments requiring code. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted.Justification: We read the NeurIPS Code of Ethics and confirm our research conform the code.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction).In this paper, we consistently employ the symbol C as a constant, which may vary from line to line. 


Decision: This paper:
* studies an interesting new setting (solving nonlinear PDE systems with multiple solutions)
* proposes an interesting new technique (Newton-informed neural operators)
* supports the technique with theoretical analysis as well as challenging experiments
* was reviewed positively by several experts in scientific machine learning. 

Given the above points, this should have been a fairly clear case! However, the initial scores were low, and the reason might have been the presentation. The authors provided elaborate responses to the reviews which helped clarify several points. Overall, I feel that the presentation is a fixable issue, and overall the paper is of high quality. 

For future revisions, I would encourage the authors to significantly revise Section 3.1 of the manuscript (it's far too terse), and explain very clearly (perhaps with a figure?) exactly what the Newton-informed neural operator is doing. I would also perhaps clarify the "multiple solutions" part up front a bit more clearly; really, the method speeds up the inner loop in Newton's method, enabling more rapid exploration of the solution space.

Reviews: review: paper_summary: paper_summary: The paper proposes the Newton Informed Neural Operator, a novel method that leverages classical Newton methods and neural network techniques, to efficiently learn multiple solutions to nonlinear PDEs.

Classical Newton's method iteratively linearizes the nonlinear equation and solves the resulting linear system to find one solution. The proposed Newton Informed Neural Operator incorporates DeepOnet to learn the operator solution to these linearized equations. After training, the nth iteration of the trained neural operator will approximate one solution.

The paper demonstrates the effectiveness of the proposed Newton Informed Neural Operator through experiments on various nonlinear PDEs. Notably, it accurately captures multiple solutions with a significant reduction in computational cost compared to traditional methods.
strength_weakness: strength_weakness: 1. The paper may lack a comprehensive comparison with existing methods. It is stated that applying the classical solver to solve Eq. (1) multiple times can be time-consuming, however,  the proposed method requires date generation, NN training, and several iterations of the trained neural operator. Without a thorough comparison, it may be challenging for readers to gauge the novelty and effectiveness of the proposed method.

2. The symbols used in the text are somewhat confusing and require careful revision. 
2.1. In Section 3.1, $\tilde{u}=u+\delta u$ is obtained by solving Equation 2, however, $u^*$ is the solution of Equation 2 in Section 3.2.
2.2. In the fourth term of Assumption 1, how do we define the H^2 norm of u-Pu? What is $\bar{P}u$? The reviewer guesses that the paper may be describing something similar to the Schauder basis and canonical projections. But the wording here is not rigorous.
2.3 What is the definition of $f(u)$ and $f'(u)$. Is $f$ just a function like $u^2$ rather than an operator?

3. The method  only needs to learn the operator solution to the linearized equations. It seems that this is a task that can be completed by a standard DeepOnet. Is there any difference between the neural network part and the approximation and generalization analysis used in this paper and those of the existing DeepOnet?
questions: questions: 1. The paper highlights that the proposed method requires fewer computational times when solving multiple linear Newton systems (e.g., 500 or 5,000 systems). This is a significant advantage. What if it's just a single linear Newton system? Do we really need to solve multiple linear Newton systems simultaneously here? In this paper, the goal is to solve linear Newton systems sequentially, correct? Could you also explain the aim of Table 1? 

2. Can the paper demonstrate the computational time required to solve the 4.4 Gray-Scott model using classical Newton iteration and the method proposed in this paper? as well as their accuracy?

3. Section B.2: Is there any specific difference between the approximation ability of the DeepOnet proved in the text and the approximation ability of the general DeepOnet. It seems that the key steps 3, 4 and 5 do not involve specific operators.

4. It is stated that equation 18 is derived via proposition 1 in subsection 16 in [26] (H. Mhaskar, Neural networks for optimal approximation of smooth and analytic functions. Neural computation, 8(1):164–177, 1996).  Is this a typo? It seems that the reference mentioned does not have such a proposition.

5. Does the considered embedding operator using Argyris element satisfy Assumption 1.4? Could the paper provide a rigorous mathematical description?

6. Theorem 2: Does the symbol E refer to expectation? Could you write it out in full in integral form? Does the theorem 2 hold for any parameter $\theta$?

7. How does the paper obtain equation 25? Is there a typo?
limitations: limitations: Yes
paper_summary: paper_summary: In this paper, the authors propose a newton informed neural operator for solving PDEs. Based on the description, it seems the neural operator is designed to predict the solution to the PDE given an initial solution (equation below eq. 2 on pg 3). The neural operator is trained with 2 types of losses, mse loss and newton loss. In the newton loss, the linearized PDE as shown in Eq. 2 is minimized using the neural operator predictions. The method is demonstrated on several use cases.
strength_weakness: strength_weakness: The paper is hard to follow. The training and inference workflows of the DeepONet are unclear. The generation of the data is not clearly described. There are no baseline comparisons or ablation studies to understand and analyze the method in a better way.
questions: questions: - The authors mention 500 supervised samples are generated but corresponding to what conditions? How is the testing data different from the training data? What happens in the scenario where the initial solution is slightly perturbed, how does the prediction and converged solution from the newtons solver compare?
- Error metrics between the solution generated by the newtons solver and the operator predictions are need to be provided for the solution generated. How do you validate the correctness of the solutions generated from the neural operator?
- In the results reported in Table 1, how does the neural operator take only 1e-4 seconds to perform 500 or 5000 iterations? Is it because it is converging faster? The information presented in the table is unclear.
- Based on the time reported in Table 1 indicate that the neural operator can directly predict the solution of the PDE given any initial condition whereas the newton solver has to iteratively solve it. Or, is the neural operator also recursively predicting the solutions? If so is this recursive operation performed while computing the Newtons loss as well? How is the memory managed in that case? Many of these details are unclear in the paper.
- The newton loss requires the computation of the adjoint and solution vector product which can be highly memory and compute intensive as the size of the grid increases? Can the authors provide details on how their method scales with increasing grid sizes?
- What is deepPOD network in Figure 2? Does it refer to the POD calculation in the trunk net? How will this POD calculation be carried out for a grid with million cells? Would it blow up the size of the trunk net? Would this make the neural operator evaluation slower as the grid size increases?  
- How do results from method 1 and method 2 compare for the different use cases?
limitations: limitations: NA
paper_summary: paper_summary: The authors presented a new machine learning based technique called the Newton-informed Neural Operator for solving PDEs that have multiple solutions due to nonlinearities. This is an important problem, and the Newton informed neural operator is capable of obtaining these multiple solutions using a single, unified training process. The authors present experimental evidence showing that their method indeed is capable of recovering these multiple solutions.
strength_weakness: strength_weakness: Unfortunately, this paper lacks clarity of exposition and the Newton-informed Neural operator is never clearly written down nor is its architecture illustrated anywhere. It appears to be a standard DeepONet augmented with a specific loss function, but the derivation of that loss function is unclear. I also don't see clearly how Equation 2 represents Newton's method, since a proper linearization would produce Jacobian terms involving both L and f. A derivation would have been helpful, but this is missing. The authors also do not compare against HomPINNs, which are designed for learning multiple solutions; if these can be used only in a supervised setting, why not compare in that setting? Finally, in addition to these serious weaknesses, I also could imagine other ways of solving this problem (such as an interval Newton iteration on a standard DeepONet).
questions: questions: 1. Introduction on page 1: are the multiple solutions coming from different possible choices of the coefficient functions, or just from the nonlinearity?  The introduction says "singular" solution when they mean unique solutions. This is somewhat problematic, since "singular" also has a mathematical meaning. I would reword.
2. Line 37 page 2, "Newton methods provide well-defined locally": is there a missing word?
3. Line 42 page 2, "Our approach combines the classical Newton method, which...": combines it with what?
4. Lines 92-95, page 3: I'm struggling to understand what the authors are saying by "remains uncontrollable". Reword for clarity?
5. Lines 99-103, page 3: Is this true? If a solution exists, there must be an operator, surely? Citations and/or proof needed.
6. Equation 2: How is this Newton's method? I see no linearization of the nonlinear operator, no Jacobian term corresponding to the nonlinearity, and only an f'(u). Perhaps I'm missing something. 
7. Page 4, Line 126: "In this paper, we will use DeepONet to approximate the Newton operator". What exactly is the Newton operator? None of the operators written out so far have been called that.
8. Line 130, page 4: "Furthermore, MgNO is replaced by...". This whole sentence is odd. Are the authors trying to say that different choices of W lead to different operator learning methods? If so, this makes sense but needs rewording.
9. Page 4, remark 1: DeepONet hasn't even been described yet. How are general readers supposed to know how this fits into this paper?
10. Page 5, Section 3.3: I still haven't seen an architecture or description of the Newton-informed Neural operator, but we're already talking about a loss function.
11. Page 6, Section 3.3.2: We finally get to the "Newton loss". It looks like a Newton-informed operator is a DeepONet trained with the Newton loss? But where did the loss in Equation 9 come from, how was it derived? I don't see a clear connection to the previous exposition.
limitations: limitations: The authors do not explicitly have a limitations section, but they briefly address their limitations in their conclusion.
paper_summary: paper_summary: The manuscript concerns an operator learning technique for elliptic partial differential equations with nonlinear, solution-dependent forcing terms. These PDE  can admit multiple solutions, which is a challenge for many existing PDE solution techniques, because they usually only result in a single solution. The authors propose to approximate the action of Newton's method, which is usually used to solve PDEs of this type, by a Deep-Operator-Network like architecture, so that the Newton steps can be performed more efficiently. They demonstrate the efficiency on multiple examples, and show results on convergence to the correct Newton iteration.
strength_weakness: strength_weakness: Major:
1) The paper does not address the issue that "multiple solutions as output of the operator" effectively means "the input is mapped to a distribution of solutions". This framework has not been addressed adequately in the PDE learning community, but of course a very large portion of machine learning concerns exactly this problem: generative modeling (i.e., diffusion models, score based learning, image generation in general, etc., all try to "learn a distribution"). Effectively, a diffusion model is a "nonlinear operator" that maps an initial distribution to a final distribution (usually, in image space). The same is true even for simple variational autoencoders.

2) An even bigger issue is that there is no explanation as to why approximating the result of Newton's method leads to a "better" (in whatever way) approximation of PDEs with *multiple* solutions. Newton's method is deterministic, and there is no guarantee that it will eventually (for some reason) "sample" all solutions. The introduced "Newton Operator Network" architecture seems to approximate the action of Newton's method on a given state, so that the operation can be performed more efficiently after the operator is trained. It does not seem to help to find multiple solutions.

3) The approximation of Newton's method has been addressed in earlier work already, which is not cited.

 [A] Doncevic, D.T., Mitsos, A., Guo, Y., Li, Q., Dietrich, F., Dahmen, M., Kevrekidis, I.G., 2024. A Recursively Recurrent Neural Network (R2N2) Architecture for Learning Iterative Algorithms. SIAM J. Sci. Comput. 46, A719–A743. https://doi.org/10.1137/22M1535310

 [B] Chevalier, S., Stiasny, J., Chatzivasileiadis, S., 2022. Accelerating Dynamical System Simulations with Contracting and Physics-Projected Neural-Newton Solvers, in: Proceedings of The 4th Annual Learning for Dynamics and Control Conference. Presented at the Learning for Dynamics and Control Conference, PMLR, pp. 803–816.



Minor:
1) l.30: it is also important to mention that operator learning techniques typically require a training data set, while methods that directly solve the PDE typically do not.
2) l.37: "Newton methods provide well-defined locally" is missing a word
3) l.50: "solve the inverse problem" is not clear at this point. Which inverse problem? 
4) The writing could be improved in general, with some expressions like "As mentioned earlier" (l.41) and "The following paper is organized as follows" (l51) being unnecessarily complicated. Section 2 "Backgrounds and Relative Works" should probably be "Background and related work". I do not list all language mistakes here.
5) Section 2.1 is not necessary, the information can be incorporated more concisely in the introduction.
6) l106,107: the index "i" is not used in the equation, same with "u_0" and "u dash".
7) l140 P_dash is not explained.
8) l172: this is not a complete sentence.
9) l244: L(u) cannot be $-\Delta u-u^2$, because it has to be a linear operator (as defined in line 20 of the manuscript). The example is still admissible, because the term $u^2$ can be absorbed into f(u).
10) l246-l248 contain multiple mistakes in the sentence structure.
questions: questions: 1) How does the approximation of the action of Newton iterations lead to a sampling of "all" (or at least more than one) solution of the PDE?

2) What are other methods that can be used to obtain multiple solutions to the given class of PDE, and what would it entail to learn them? Are they slower than the solution proposed here?

3) In case a specific PDE has to be solved (ideally, with several of its solutions as a result), does it make sense to train the Newton operator first, or would it be more reasonable to just use the classical Newton's method, because it can be started immediately without training the network?
limitations: limitations: The main paper does not contain any statements on specific limitations, or a broader discussion on potential negative impact of the work. This must be addressed.

Comments: title: Official comment
Official comment
Comments on rebuttal
comment: We appreciate your acknowledgment of the potential of our method. We will revise our paper based on our discussion. Thank you once again for your valuable insights and constructive feedback.
Thanks for the clarification. Please make the necessary revisions to the paper to improve the presentation of the ideas and results. I have increased my score to 5.
Thank you for your comments.

In Figure 3, the training loss and the absolute L2 test error (which is large because the magnitude of the true solution is large) are not directly comparable, as they represent different metrics. The training loss reflects the rescaled Newton's loss, as we have demonstrated. Therefore, their magnitudes should not be compared directly. The reason we plot them together is to illustrate that their trends during the training process are consistent. This consistency suggests that the unsupervised data, when trained with Newton's loss, contributes to the improvement of the model's performance. We will revise the figure to avoid any confusion.

Regarding convergence, it is important to note that we cannot guarantee it, as even the exact Newton's method does not ensure convergence for every unsupervised initial condition. For other unsupervised initial conditions, if the unsupervised data contains the solution and has enough data points, we can achieve convergence by effectively training with the Newton loss. However, if there is not enough unsupervised data for certain types of solutions, we cannot guarantee the same level of success. Nonetheless, our experiments demonstrate that adding unsupervised data expands the exploration range of initial conditions.

We will also enhance the clarity of Table 1 in the revised version.
Thanks for the clarification.

For example, if we look at Figure 3, the training loss is around 0.01 while the testing is of the order 1.0. I completely understand that the testing error improves with the addition of the Newton loss but my main question is whether that improvement is sufficient considering this discrepancy between the train and test loss? If the testing errors are large then how can this method be used to explore a range of initial conditions as the authors claim? I understand that the ring pattern shown in Fig. 3 is completely unsupervised. Are there other unsupervised initial conditions that have been analyzed? What guarantees that the method will always converge and the testing errors will remain bounded for any unsupervised initial condition?

I think the text in Table 1 needs to be modified because it is confusing. It gave me the impression that you are solving 1 initial condition but performing 500 and 5000 linear iterations in the Newton's solver. However, it seems like you are solving 500 or 5000 initial conditions until convergence.
We sincerely thank the reviewer for the insightful feedback.

1. The supervised training data includes perturbations on 10 solutions, representing only a small portion of the entire training dataset. The unsupervised data, on the other hand, is sampled from a much larger space with greater variance. This intentional setup is to demonstrate that only supervised data is not enough (also related with your second question "why test error is much larger than training error").  Our model is not limited to small perturbations around known solutions but is also effective in larger space. This effectiveness is due to the unsupervised data being sampled and trained using Newton's loss function economically. This is the primary motivation for using Newton's loss. In lines 287-290 and Figure 3(a), we show that our model, trained with unsupervised data, can successfully handle pattern formations that are not present in the supervised dataset. In that figure, the test case (a ring pattern) differs from the 10 patterns considered in the supervised dataset, but rather close to unsupervised training dataset.

2. Your questions align precisely with what we aim to demonstrate in Figure 1 and lines 237-241. In Figure 1(a), the significantly larger testing error compared to the training error suggests that using only MSE loss is insufficient (model is only valid at small pertubation around 10 solutions in this case). Figure 1(b) illustrates that increasing the number of unsupervised samples leads to a substantial improvement in the test L2 error, highlighting the necessity of introducing Newton's loss. The testing samples are drawn from the same distribution as the unsupervised training samples.

3. These operations are parallelized, meaning that we predict all 500 solutions simultaneously. The parallelization of solving 500/5000 Newton steps, whether using a classical Newton solver or a neural operator, is reasonable because we are addressing 500/5000 problems with 500/5000 different initial conditions. For a fair comparison, the classical Newton solver is also parallelized using CUDA implementation on a GPU. The speedup of the neural operator arises from its ability to benefit more from batch processing, as it naturally handles large batch sizes efficiently during the inference phase. We discuss this in more detail in Appendix A.5. This leads back to our primary motivation: enhancing the computational efficiency of finding multiple solutions by enabling extensive exploration of initial conditions. This significant speedup allows for practical exploration of a much larger space of initial conditions within a reasonable timeframe. We further discuss our motivation and methodology in the global response and response to Reviewer zTaB.

We thank the reviewer again for the valuable discussion and will revise the paper accordingly to improve the clarity.
I would like to thank the authors for their revision. I understand that training the neural operator with a Newton's loss is novel and the main goal of the paper is to demonstrate that neural operators can be trained for multiple solutions but I have a few followup questions.

If the training and testing data set is generated by considering small perturbations on 10 solutions then, is the validity of the trained model only within these bounds of the initial solutions? The authors claim that " This significant speed-up allows for the practical exploration of a much larger space of initial conditions in a reasonable timeframe." Where is this claim verified in the paper? Is the testing carried out on initial conditions + perturbations different from the 10 considered in the Gary Scott use case?

The errors reported in Figs. 1 and 3 are about 1-2 orders of magnitude larger than the training errors. Is this considered as a reasonable error? Why is that the case when the testing and training samples are from the same distribution? Do the testing results improve if you increase the number of training samples?

The neural operator evaluation is recursive meaning a solution at "n+1" is calculated based on the solution at "n". Table 3 reports that 500 evaluation steps require 1.1e-4 seconds and for 5000 evaluation steps the time increases to 1.4e-4. Why is that the case? Why doesn't the time increase 10x like the Newtons solver? The time taken by each step is constant and proportional to the number of parameters in the NN and by that logic it should be more for 5000 steps vs 500 steps. Or are these operations batched meaning you are predicting all the 500 solutions at once? If that is the case then it is not a correct representation because the solutions at future iterations are unknown during inference. Am I missing something?
Thanks for the reviewer's positive response. Here we want to mention that although we use DeepONet to learn the operator form of linear PDEs, there are still some differences between our work and the standard approach.

First of all, in our work, the loss function introduces the residual error from PDEs, whereas the standard DeepONet is purely a supervised learning method as described by Lu Lu, Pengzhan Jin, and George Em Karniadakis in DeepONet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators (arXiv preprint arXiv:1910.03193, 2019).

Secondly, in the proof of approximation, our input includes not only $ u $ itself but also the second derivative of $ u $. This requires that the encoding in our input part should be smoother. This is why we need an encoding like Argyris elements (Assumption 1 (iv)) to add the extra derivative information to achieve the approximation. The well-definition of the operator needs to be checked since this is not a smooth map from $ L^2 \to H^2 $ like in the standard DeepONet task, where regularity increases. In our case, we are dealing with $H^2 \to H^2$, which requires more assumptions and discussions to ensure the operator is well-defined. This is addressed in Steps 1 and 2 of the approximation part. Therefore, we need methods like Argyris elements for input and need to ensure the smoothness of $ f $ and the eigenvalue distribution of $ \mathcal{L} $.

Lastly, regarding generalization error, we consider not only the generalization error for $ L^2 $-loss (Theorem 1) but also for $H^2 $ (Corollary 1). The generalization error for $H^2$ requires a similar proof framework but with more stringent smoothness assumptions for the operator and input space compared to the generalization error for $ L^2 $-loss (Theorem 1). The standard DeepONet only considers $L^2 $-loss and does not include this part.

In summary, our paper uses DeepONet incorporating derivative information, unlike the standard DeepONet. We will discuss these differences and provide more details in the revised paper.
Dear Reviewer,

As the discussion period draws to a close, we kindly invite you to reassess your review in light of our responses. Please let us know if you have any remaining questions or concerns.

Best wishes,
Authors
This makes a lot of things much clearer. Assuming that the authors commit to revising the paper with the details above, I raise my score.
Thank you very much for the  response. However, the third point of weaknesses is not fully addressed. I understand the contribution of designing Newton-based loss function to address nonlinear pde. However,  the final task is to learn the operator solution to the linearized equations.   Therefore I would like to see the difference between the neural network part and the approximation and generalization analysis used in this paper and those of the existing DeepOnet. Having said that, I feel this paper is a valuable work in the field of operator learning and the broader scientific ML, and I maintain my original positive score as is.
We are pleased that the revisions have clarified the novelty of our approach. We appreciate your recognition of the potential of our method. Thank you once again for your valuable insights and constructive discussion.
We would like to further clarify the novelty of our method:

1. We introduced a general computational framework (as illustrated in the global response) for learning the Newton's nonlinear solver using neural operators integrated with the Newton-informed loss, with DeepONet/DeepPOD (architecture detailed in Appendix A.3) as a representative example. This computational framework can be used to identify multiple solutions of nonlinear PDEs.
    
2. After a thorough review of the HomPINNs paper, we identified a key distinction between our method and theirs. In HomPINNs, networks are primarily used to approximate/parametrize the solution. In contrast, our method leverages networks to approximate the Newton solver—specifically, the mapping from $u_k$ to $\delta u$,  $u_{k+1} = u_k+ \delta u$. While HomPINNs are designed to compute multiple solutions starting from simple initial functions, typically on coarse grids, by combining the homotopy approach with PINNs to learn multiple solutions on finer grids, our Newton-informed operator learning method directly trains the nonlinear solver on fine grids. This approach offers greater efficiency and reduces the risk of overlooking solutions that might be missed by HomPINNs.
I thank the authors for the careful revision. The explanation in the general comment was very helpful as well - it seems the major issue was the presentation of the results, not the idea or results themselves. I read the paper again and I am now more excited about the idea and results. Indeed, learning Newton's method in this context makes a lot of sense - if the classical approach requires millions of calls to a Newton solver (for the millions of initial conditions), then replacing this inner loop with a learned model is great. I updated my score to 6.