Newton Informed Neural Operator for Solving Nonlinear Partial Differential Equations
Authors: Wenrui Hao
Xinliang Liu
Yahong Yang
Venue: NeurIPS 2024 poster (2024)
Keywords: operator Learning
nonlinear partial differential equations
multiple solutions
Newton method
Primary area: machine_learning_for_physical_sciences
TL;DR: We develop a Newton-informed neural operator for solving nonlinear problems with multiple solutions, requiring less data and offering faster convergence than traditional Newton methods.

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.

Newton Informed Neural Operator for Solving Nonlinear Partial Differential Equations
Wenrui Hao"
Department of Mathematics
The Pennsylvania State University, University Park,
State College, PA, USA
wxh640psu.cc edu
Xinliang Liu
School of Mathematical Sciences, Shenzhen University, Shenzhen, China King Abdullah University of Science and Technology. Thuwal, Saudi Arabia xinliang. liu@kaust... edu. sa
Yahong Yang Department of Mathematics
The Pennsylvania State University, University Park, State College, PA, USA yxy54980psu.cm edu
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 meth- ods, 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 each iteration. This approach allows to compute multiple solutions single learning process while requiring fewer supervised data points than existing neural network methods.
Introduction
Neural networks have been extensively applied to solve partial differential equations (PDEs) in various fields, including biology, physics, and materials science (2. While much of the existing work focuses on PDEs with unique solution, nonlinear PDEs with multiple solutions pose a significant
challenge (39 10 but are widely encountered applications such 3T 36 14 [3] In this
paper, we aim to solve the following nonlinear PDEs with multiple solutions:
(L))(a)) = f(u(x)), x E Ω [u(x)) 0.
(1)
E UN,
"Corresponingg Authors
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

Here, Q is the domain of equation, f(u) is nonlinear function in R. u Rd + R and L is second-order elliptic operator given by Lu Eij= a's (z)zzz Σ (()) c(z)u, for given coefficient functions a²) "2,,b,,cii, c(i.j = d.....d) with Σij=1 a (z)(1E) 2 X|E?, for a constant À 2 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 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 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 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 [). 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 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 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.
2

---

2 Backgrounds and Relative Works
2.1 Nonlinear PDEs with multiple solutions
Significant mathematical models depicting natural phenomena in biology, physics, and materials science 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, fraction-diffusion systems in chemistry, and pattern formation in biology (20) [2 However, experimental techniques like synchrotronic and laser methods can only observe subset of these multiple solutions. Thus, there is 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 targe-scale systems encountered in real-world applications. response to these challenges [1], we propose an operator learning approach based on Newton's 's method to efficiently solve nonlinear PDEs.
2.2 Related works
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.
Function learning methods In function learning, a commonly employed method for addressing this problem involves the use of Infoiee-I 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 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 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, present operator learning approach combined with Newton's '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 [4] based on multilevel methods, and transformer-basee 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

---

3 Newton Informed Neural Operator
3.1 Review of Newton Methods Solve Nonlinear Partial Differential Equations
To tackle this problem Eq. ① we employ Newton's 's method by linearizing the equation.
For the Newton method applied operator, if aim to find the solution of F(u) 0, the iteration can be written as:
F (un) un+1 F (un) un Fuu)) = F (un) ou - (un),,
where du Un+11 Un.
In this context, F'(u)o is the (Fréchet) derivative of the operator, which is a linear operator to defined as follows: To find F") in X, for any X.
F(u + v) F(u) F')eel =0,
lim -||
|u|
where | denotes the norm in X.
given any initial guess uz) for = 1,2,.... M, in the i-th iteration of Newton's method, we have ū(x) u + 6u(x)) by solving
For solving Eq. I
[C - (l))) - Lu + f(u), du 0,
a E
(2)
E UΩ,
which is based on the fact that the (Fréchet) derivative of f(.) at unique solution, then by solving Eq. (2) and repeating the process M times, we will obtain one of the
f'(u). Eq. 2 has
solutions the nonlinear equation I Denoting the mapping for u and 8u, the solution of Eq.
with parameter u, Guu) ou, we know that
lim (9 (d)())(u)))))))))))))))))))))( = u", 12 22000
where u" is one of the solutions of Eq. T.. For different initial conditions, this process will converge to different solutions of Eq. (I), making this method suitable for finding multiple solutions. Furthermore, the Newton method is well well-posed, meaning that each initial condition un will converge to a single solution of Eq. ( under appropriate assumptions (see Assumption T). This approach helps to address the ill ill-posedness encountered when using PINNs directly to solve Eq. (T) However, repeatedly solving Eq. ①1) can be computationally expensive, especially in high-dimensional cases when large number of discrete points are involved. In this paper, we tackle these challenges by employing neural networks.
3.2 Neural Operator Structures
In this section, we introduce the structure of the neural operator to approximate the operator locally in the Newton methods from (2, i.e., du := ((u), where du is the solution of Eq. which depends on u. If we can learn the operator (uu) well using the neural operator (u; 0). then for an initial
function assume the n-th iteration will approximate one solution, i.e., (g + Id) (uo) z u". Thus,
(O + Id)( (m)))) (wo) (g + Id)(")(o)) 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 neurons for operators from X' to ) as
(aa; 0) = Σ Aso (Waa + B.) Ya E X
(3)
where W; Here, L(X,) denotes the set of all bounded (continuous) linear operators between A' and y and H+ defines the nonlinear wise activation.
CO,,), and 0 denote all the parameters in (W,, A;, B, --

---

In this paper, we will use shallow DeepONet 229 22 to approximate the Newton operator. To provide more precise description, the shallow neural network, W; 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 multigrid operator [3 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 C ( of the operator and fu) in Eq. (T). The definition of the Sobolev space can be found in Appendix B.1 Assumption 1. (i): For any u E X, we have that the linear equation
(c- f' -Lu + f(u)
is well-posed for solving 8u.
(ii): There exists a constant F such that [[(())ll F.
(ii:): All coefficients in L. are C and an C2.
(iv): has a Schauder basis {b we define the canonical projection operator Pr based on this basis. The operator Pn. projects any element E X onto the fiiteedimeesional subspace spanned by the first n. basis elements {bi, b2....,b,ll b.). Specifically, for u E X, u = Ek=o Ork by. let Pn(u)) = -
where og are the coefficients in the expansion of u with respect to the basis {bn ). The canonical projection operator Pn, is linear bounded operator on X. According to the properties of Schauder bases, these projections Pn, are uniformly bounded by some constant C. Furthermore, assume, for any X, 0, there exists such that
Iu P, P (2)))) S c,
for all X.
More discussion about the assumption shown the appendix. The sketch of the proof is illustrated in Fig.
g
X'
X'
Encoding Step 4
Reconstruction Step 3:
Functional Approximation
- Wa
(a,(0)) Σ (xu)o(wx + ()
12
IF
NN Approximations Step
(n) = A; a(M) u+b
Figure 1: The sketch of proof for Theorem T
Theorem 1. Suppose X = c (2) and Assumption 1 holds. Then, there exists neural network
Ou: 0) E E defined as
Aio (Wiu b.) a (wi IM L(X, R"),, E RM A₁ R Ix
(4)
such that
p Σ
A (W.u + bi) (w, x x (1) G(u) ((w)
S Cim + (5)
sup ue
(P(())
where a smooth non-pollnomial activation function, n. is shown in Assumption Aadd contained in W., C₁ is a constant independent of m, and C2 is a constant depended on p. n and F (see
in Assumption 7 is the scale of the P in Assumption 7 And depends on n. n and e are defined in Assumption 7
The approximation results of DeepONet in Sobolev training can be found in 44
5

---

3.3
Loss Functions of Newton Informed Neural Operator
3.3.1
Mean Square Loss
The Mean Square Error loss function is defined as:
M M2
(s( = M M.
99 (4,) (xx) O (uy; (xx))
(6)
where 11, .2..... UM. μ are independently and identically distributed (i.i.d) samples in X, and 3i, 22........ M,, are uniformly i.i.d samples in f2.
However, using only the Mean Squared Error loss function is not sufficient for training learn M. the Newton method, especially since in most cases, do not have enough data {uy: 9 (ug) (=)) j=1- Furthermore, there are situations where we do not know how many solutions exist for the nonlinear equation ( If the data is sparse around of the solutions, it becomes impossible effectively learn the Newton method around that solution.
Given that (O) can be viewed as the finite data formula of (O), where
lim ( M..MM.
(O(
The smallness of Ese can be inferred from Theorem I To understand the gap between (g) and (O)) we can rely on the following theorem. Before the proof, we need some assumptions about the data in (0):
Assumption 2. (i) Boundedness: For any neural network with bounded parameters, characterized
by a bound B and dimension de, there exists function V (22) [0. oo) such that
[2()) < (u),
sup
(x))) ((u), sup Ic0(u; ()())) < Vuu) 0E[-B.B]*e
BE(-B.B|d0
for all u E X Ω and there exist constants C, O, such that (u) C1 ||u||||llll m²)"
(ii) Lipschitz continuity:
There exists a function + (22 - [0, 20), such that (w: (x) Ou; (ax) < ((u) ||0 O'
(8)
for all E X, E St and (uu) < C (1 (1+ (a))) for the same constants C, 0 as in Eq. 77
(iii) Finite measure: There exists > O, such that
J é allul du((u) < 20.
Theorem 2. If Asssupption 2 holds, then the generalization error is bounded by
((Es(0) Es)( C
1 VM,
(1+Cdp L,,
do log Mg
sup 0E|-B,B|
VM,
where C is a constant independent of B, do. Mg. and M,. The parameter is is specified in 7 Here, B represents the bound of parameters, and dg is the number of parameters.
The proof of Theorem is presented Appendix B.3
Remark 1. Assumption 2 is easily satisfied if we consider X as the local function set around the solution, which typically the case in Newton's 's methods. This aligns with our approach and the working region in the approximation part (see Remark 3). The error supee]-B.Be
(0) Esc(O)) suggests that the network can perform well based on the loss function (())). The reasoning is as follows: let Os = arg mines -B.B (e( and Os. arg mineei-B,.BB (o) We aim for (Os) to be small, which can be written as:
(0s) < (0s.) + ((Es(0s) - Esc(Os)) I Eso(Os.) +
sup B.
IEEs (s)()))
where (0s.) is small, demonstrated by Theorem 1 when B is sufficiently large.

---

3.3.2 Newton Loss
As we have mentioned, relying solely on the MSE loss function can require 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 referred to here the Network loss function.
The Newton loss function defined as:
No N₂
1
(0) = Nu. j=1k=1
L (C- (( (ax) + f(u))
(9)
where u1, ....... UNN are independently and identically distributed (i.i.d) samples in N, and 21, -.,..... ,N, are uniformly i.i.d samples in Ω
Given that (0) can be viewed as the finite data formula of (0) where (0)
(0) Nu.Na-to
lim
To understand the gap between (0) and yn (0), we can rely on the following theorem: Corollary 1. If Assumption hoolds, then the generalization error is bounded by
1 LVN.
(1+ (oCCooogCBB//.)) 2xxt122) de davlog N₂
sup BE|-B,B|*
(0) Ex)(O)) C
VN²
where C' is a constant independent of B, de. Nz. and Nu. The parameter h. is specified in 7 Here, B represents the bound of parameters, and dg is the number of parameters.
The proof of Corollary is similar to that of Theorem 2 therefore, it will be omitted from the paper.
Remark 2. If we only utilize (0) as our loss function, as demonstrated Theorem 2 we require both My and M2 to be large, posing a significant challenge when dealing with complex nonlinear equations. Obtaining sufficient data becomes a criticul issue in such cases. In this paper, we integrate Newton's 's information into the loss function, defining it as follows:
{() == XEs(0) Ex
(10)
where (0) represents the cost associated with the unsupervised learning data. If we lack sufficient data for (0). can adjust the parameters by selecting small X and increasing No and Nu. This strategy enables effective learning even when data for (0) is limited. We refer to this neural operator, which incorporates Newton information, the Newton Informed Neural Operator.
In the following experiment, function in Eq. (10) to learn one step of the Newton method locally, i.c., the map between the input and the solution du in eq. (2). If we have a large dataset, we can choose large À in (()) (10: if we have a small dataset, we will use a small À to ensure the generalization of the operator minimized. After learning one step of the Newton method using the operator neural networks, quickly obtain the solution by the initial condition of the nonlinear PDEs T and find new solutions not present in the datasets.
will use the neural operator established in Eq. (3) and the loss
easily and
4 Experiments
4.1 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 '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 practical scenario with limited data availability, facilitating a comparison between these training strategies to evaluate their efficacy 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
7

---

4.2 Case 1: Convex problem
We consider 2D convex problem f()) = =0, where L(u) = f(u) -u² sin 5 ( and u = 0 on 0R2. 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 with À = 0.01. Please refer| A.1. for the dataset generation of convex case. For the detailed experimental settings, refer to Appendix
100
L2 test error deepONet-MSE
H1 test error of deepONet-MSS MSE
L2 train arror of deepONet- MSE
L2 test error deepONet-Neerinn
MSE logarithmic scale
L2 test despONet-MSE
H1 test error NetpONet-Newtor
100
10
0.1
200
400
600
800
200
400
600
800
Epoch
Epoch
(a) Training and testing errors using Method (b) Comparison of models trained using Method 1 (deepONet-MSEE) and Method 2 (deepONet-Nowton)
1 (MSE loss for 500 supervised data)
Figure Training and testing performance of DeepONet under different conditions.
MSE Loss Training (Fig. 2(a)) In method 1, Effective training is observed but exhibits poor generalization. The significantly larger testing error compared the training error suggests that using only MSE loss is insufficient. Performance Comparison (Fig, 2(b):: DeepONet-Newtor model (Method exhibits superior performance in both L2 and H1 error metrics, highlighting its enhanced generalization accuracy. This study shows the advantages of using Newton's 's loss for training DeepONet models, illustrating that increasing the number of unsupervised samples via introducing Newton's 's loss leads to substantial improvement in the test L2 error and H₃ error.
4.3
Case 2: Non-convex problem with multiple solutions
We consider 2D Non-convex problem,
- Au(x, y) u u²(x,, v) sin( in 2 u(( y) = 0, in an
(11)
where Q (0, (0, 1) B. In this case, L(w) == -(w), f(u) = n² sin(Tra) sin(Ty) 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 will saturate 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 Figure 3). On the other hand, the model trained via method is less affected by the limited supervised data since the utilization of Newton's '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 's method. Both methods parallelize operations to solve 500/5000 Newton lincar systems simultancously, each with distinct initial states. The key advantage of the
8

---

L2 teal desyPCD 12
10%
111 test ermor deepPOC-2
L2 test enor deepPDD Newton
test erron SexpPOC-Newiinn
100g
200
200
400
600
830
(a) Solutions of 2D Non-convex problem
(b) Method 1 VS Method 2
Figure 3: Solutions of 2D Non-convex problem
neural operator lies its ability to batch the computation these independent systems efficiently. By efficiently batching and sampling 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 rapid computation and extensive initial condition sampling enhances the chances of identifying multiple solutions.
For 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 A.6
Parameter
Newton's Method
NINO
Number of Streams
10
Data Type
float32
float32
Execution Time for 500 linear Newton systems (s)
31.52
1.1B-4
Execution Time for 5000 linear Newton systems (s)
321.15
1.4E-4
Table 1: BeBnhhmrrking the efficiency of Newton Informed Neural Operator. Computational Time Comparison for Solving 500 and 5000 Initial Conditions.
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 NIÑO, solving 5000 independent Newton linear systems scales up minimally compared to solving 500 systems, while the classical solver experiences 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.
4.4 Case The Gray-Scott model
The Gray-Scott model 31
describes the reaction and diffusion of two chemical species, A and S. governed by the following equations:
11
0A
Øt D SA² + (μ+ p)
0S
DsAS + SA² - p(1 - S),
Øt
where DA and D3 are the diffusion coefficients, and μ and p are rate constants.
Newton's Method for Steady-State Solutions Newton's 's method is employed to find steady-state
solutions AA = 0 and
dt = 0) by solving the nonlinear system:
dr
DAAA 5A²2 (1+ p)
(12)
DsAS + SA² p(1 S).

---

The Gray-Scon model is highly sensitive to initial conditions, where even minor perturbations can lead to vastly different emergent patterns. Please refer Figure some examples of the patters.. This sensitivity reflects the model's s complex, non-linear dynamics that can evolve into multitude of possible steady states based on the initial setup. Consequently, training neural operator to map initial conditions directly to their respective steady states presents significant challenges. Such model must learn from a vast functional space, capturing the underlying dynamics that dictate the transition from any given initial state its final pattern. This complexity and diversity of potential outcomes is the inherent difficulty training neural operators effectively for systems as complex the Gray-Scon model. One refer to A.1 for detailed discussion the Gray-Scott model. We employ Neural Operator as a substitute for the Newton solver in the Gray-Scott model, which recurrently maps the initial state to the steady state.
RZerrow
(a) An example demonstrating how the neural operator maps the initial state to the steady state in iterative manner
0.012
meen L2 errar
Training lows
L.2 error
0.01
crro
0.006
0.004
0.002
The number Neural Operator iteration
(b) Average convergence rate of NINO.
(c) Training via method
Figure 4: The convergence behavior of the Neural Operator based solver.
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 '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 's loss, can effectively approximate the mapping of the initial solution to its correct steady state. we further test our neural operator, utilizing as a surrogate for Newton's method address nonlinear problems with initial state drawn from the test dataset. The curve shows the average convergence rate of ||u - uxll across the test dataset, where и; represents the prediction at the i-th step by the neural operator. In subfigure (c), we compare the Training Loss (Rescaled Newton's 's Loss) and Absolute I.2 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.
5 Conclusion
In this paper, develop neural operators learn the Newton's solver related nonlinear PDEs (Eq. (T) 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 Newton- informed neural operator. We provide 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 ad advantages of our network mentioned earlier.
10

---

Reproducibility Statement Code Availability; The code used in our experiments can be ac- cessed via /ttttt://ggtthii com/x11ii2011/eaar aeeaarn_ewttoaaee and also the supplementary mate- rial. 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.
Acknowledgments and Disclosure of F Funding
Y.Y. and W.H. were supported by the National Institute of General Medical Sciences through grant IR35GM146894, The work of X.L. was partially supported by the KAUST Baseline Research Fund.
References
[1] H. Amann and P. Hess. A multiplicity result for class of elliptic boundary value problems. Proceedings of the Royal Society of Edinburgh Section A: Mathematics, 84(1-2):145-135 1979.
[2] S. Brenner. The mathematical theory of finite element methods. Springer, 2008.
[3] B. Breuer, McKenna, and M. Plum. Multiple solutions for semilinear boundary value problem: computational multiplicity proof. Journal of Differential Equations, 195(1):243 269, 2003.
[4] Shuhao Cao. Choose a transformer: Fourier or galerkin. Advances in Neural Information Processing Systems, 34, 2021.
[5] Mark C Cross and Pierre C Hohenberg. Pattern formation outside of equilibrium. Reviews of modern physics, 65(3):851, 1993.
[6] L. Evans. Partial differential equations, volume 19. American Mathematical Society, 2022.
[7] Somdatta Goswami, Aniruddha Bora, Yue Yu, and George Em Karniadakis. Physics-informed neural operators. 2022 arXiv preprint arXiv:2207.05744 2022.
[8] R. Guo, S. Cao, and L. Chen. Transformer meets boundary value inverse problems. In The Eleventh International Conference on Learning Representations, 2022.
[9] J. Han, A. Jentzen, and W. E. Solving high-dimensional partial differential equations using deep learning. Proceedings of the National Academy of Sciences, 115(34))8505-851) 8510, 2018.
[10] W. Hao, S. Lee, X. Xu, and Z. Xu. Stability and robustness of time-discretization schemes for the allen-cahn equation via bifurcation and perturbation analysis. arXiv preprint arXiv22406.18393, 2024.
[11] W. Hao, C. Liu, Y. Wang, and Y. Yang. On pattern formation in the thermodynamically- consistent variational gray-scott model. arXiv preprint arXiv:2409. 04666, 2024.
[12] W. Hao and C. Xue. Spatial pattern formation in reaction reaction-diffusion models: computational approach. Journal of mathematical biology, 80:521-543, 2020.
[13] Wenrui Hao, Jonathan D Hauenstein, Bei Hu, and Andrew J Sommese. A bootstrapping approach for computing multiple solutions of differential equations. Journal of Computational and Applied Mathematics. 258:181 1--10, 2014.
[14] Wenrui Hao, Jan Hesthaven, Guang Lin, and Bin Zheng. A homotopy method with adaptive basis selection for computing multiple solutions of differential equations. Journal of Scientific Computing, 82(1):19, 2020.
[15] Wenrui Hao, Sun Lee, and Young Ju Lee. Companion-based multi-level finite element method for computing multiple solutions of nonlinear differential equations. Computers & Mathematics with Applications, 168: 2024.

---

[16] J. He. X. Liu, and J. Xu. Mgno: Efficient parameterization of linear operators via multigrid. In The Twelfth International Conference on Learning Representations, 2023.
[17] R. Hoyle. Pattern formation: an introduction to methods. Cambridge University Press, 2006.
[18] J. Hu and S. Zhang. The minimal conforming H* finite element spaces R" rectangular grids. Mathematics of Computation, 84(292):563-5777 579, 2013.
[19] Y. Huang, W. Hao, and G. Lin. Hompinns: Homotopy physics-informed neural networks for learning multiple solutions of nonlinear elliptic differential equations. Computers & Mathemat- ics with Applications, 121:62-73, 2022.
[20] Kondo and T. Miura. Reaction-diffussin model as a framework for understanding biological pattern formation. science, 329(5999):1116-16)) 2010.
[21] I. Lagaris, A. Likas, and D. Fotiadis. Artificial neural networks for solving ordinary and partial differential equations. IEEE transactions on neural networks, 95))987-11000, 1998.
[22] Samuel Lanthaler, Siddhartha Mishra, and George E Karniadakis, Error estimates for deeponets: deep learning framework in infinite dimensions. Transactions of Mathematics and Its Applications, 6(1 :tna:00,, 2022.
[23] Z. Li and N. Yan. New error estimates of bi-cubic hermite finite element methods for biharmonic equations. Journal of computational and applied mathematics, 142(2))251-255 2002.
[24] Z. Li, Zheng, N. Kovachki, D. Jin, H. Chen, B. Liu, Arizzadenesheli, and Anand- kumar. Physics-informed neural operator for learning partial differential equations. CoRR, abs/2111.03799 2021.
[25] Zongyi Li, Nikola Borislavov Kovachki, Kamyar Azizzadenesheli, Kaushik Bhattacharya, An- drew Stuart, and Anima Anandkumar. Fourier NeuralOperator for Parametric Partial Differential Equations. In International Conference on Learning Representations, 2020.
[26] B. Lin, Z. Mao, Z. Wang, and G. Karniadakis. Operator learning enhanced physics-informed neural networks for solving partial differential equations characterized by sharp solutions. arXiv preprint arXiv:23 19590, 2023.
[27] X. Liu, B. Xu, S. Cao, and I.. Zhang. Mitigating spectral bias for the multiscale operator learning, Journal of Computational Physics, 506: 12944, 2024.
[28] L. Lu, X. Meng, S. Cai, Z. Mao, S. Goswami, Z. Zhang, and G. Karniadakis. A comprehensive and fair comparison of two neural operators (with practical extensions) based fair data. Computer Methods in Applied Mechanics and Engineering. 393: 2022.
[29] Lu Lu, Pengzhan Jin, and George Em Karniadakis. Deeponet: Learning nonlinear operators for identifying differential equations based on the universal approximation thcorem of operators. arXiv preprint arXiv: 1910.03190 2019.
[30] H. Mhaskar. Neural networks for optimal approximation of smooth and analytic functions. Neural computation, 8(1):164-17) 177, 1996.
[31] Pearson. Complex patterns in a simple system. Science, 261 1993.
[32] T. Poggio, H. Mhaskar, I.. Rosasco, B. Miranda, and Q. Liao. Why and when can deep-but not shallow-networks avoid the curse of dimensionality: review. International Journal of Automation and Computing, 14(5):503-519, 2017.
[33] Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics Physics-ifformed informed neural networks: deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686-707, 2019.
[34] A. Schmidt-Hieber. Nonparametric regression using deep neural networks with relu activation function. Annals of statistics, 48(4): 2020.
12

---

[35] M. Ulbrich. Semismooth newton methods for operator equations in function spaces. SIAM Journal on Optimization, 13(3):805-841, 2002.
[36] Yingwei Wang, Wenrui Hao, and Guang Lin. Two-level spectral methods for nonlinear elliptic equations with multiple solutions. SIAM Journal Scientific Computing, 40(4):B (144)BB180--D 1180-B12052 1205, 2018.
[37] T. Welti. High-dimensional stochastic approximation: algorithms and convergence rates. PhD thesis, ETH Zurich, 2020.
[38] Xu. Two-grid discretization techniques for linear and nonlinear pdes. SIAM journal on numerical analysis, 33(5):155--177,, 1996.
[39] Jinchao Xu and Xiaofeng Xu. Lack of robustness and accuracy of many numerical schemes for phase-field simulations. Mathematical Models and Methods in Applied Sciences, 33(08): 1721 1746, 2023.
[40] Y. Yang. Deeponet for solving PDEs: Generalization analysis in Sobolev training. arXiv preprint arXiv:2410.04344 2024.
[41] Y. Yang and He. Deeper or wider: A perspective from optimal generalization error with sobolev loss. International Conference on Machine Learning, 2024.
[42] Y. Yang, Y. Wu, H. Yang, and Y. Xiang. Nearly optimal approximation rates for deep super relu networks on sobolev spaces. arXiv preprint arXiv:2310. 10766, 2023.
[43] Y. Yang and Y. Xiang. Approximation of functionals by neural network without curse of dimensionality. arXiv preprint arXiv:2205. 14421, 2022.
[44] Y. Yang, Yang, and Y. Xiang. Nearly optimal VC-dimension and pseudo-dimension bounds for deep neural network derivatives. NuerIPS 2023, 2023.
[45] B. Yu and W. E. The deep ritz method: deep learning based numerical algorithm for solving variational problems. Communications in Mathematics and Statistics, 6(1):1-12, 2018.
13

---

Figure 5: Examples of steady states of the Gray Scott model
A Experimental settings
A.1 Background on the PDEs and generation of datasets
A.1.1 Case convex problem
Function and Jacobian
The function F(u) might typically be defined as:
Fu)) -u -u++u"
The Jacobian J(),, for the given function Fuu). involves the derivative of F with respect to 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 ug N A -) and then calculate the con- vergent sequence (uo, "i,..., by Newton's 's method. Each convergent sequence {uo, U1 ... Un one data entry the dataset.
The analysis of function and Jacobian for the non convec problem (case 2) is similar to the convex problem except that its Jacobian J(u) A 2u such that Newton's 's system not positive definite.
A.1.2 Gray Scott model
Jacobian Matrix
The Jacobian matrix J of the system is crucial for applying Newton's method:
JAA
JAS
=
JSA
Jss
with components:
JAA DAA A + diag(- -25 + # p),
JAS = diag(- -A²),
JSA diag(
Jss Ds diag(A² + p).
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 each side, leading to total of 2 grid points. This resolution was chosen to balance computational efficiency with spatial resolution sufficient capture detailed patterns. The spacing between each grid point, h, is computed as = 1.0
N 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 DA = 2.5 x 10 and Ds 5.0 x 10 respectively. These values determine the rate which each species diffuses through the spatial domain.
The reaction rate μ and feed rate are crucial parameters that govern the dynamics of the system. For this simulation, ja set to 0.065 and p to 0.04, influencing the production and removal rates of the chemical species.
Reaction Rates:
Simulations
The simulation utilizes finite difference method for spatial discretization and Newton's 's method to solve the steady-state given the initial state. The algorithm is detailed in A.6
14

---

Figure 6: Three examples depicting the evolution from the initial state to the steady state via Newton's 's method.
A.2
Data generation
Here is how generate the supervised data samples:
Step 1:
We use a classical numerical solver to obtain single (multiple) solutions nonlinear PDEs. For example, in Case 2, there exist four solutions u' u. u², wt. We have one solution for Case 1 and 10 solutions for the Case Gray-Scott model.
2. Step 2: The supervised dataset generated by sampling perturbation around the solution u' ~ (0, (-))-3) on the chosen true solution u'. We then set ub u, + u² and calculate the convergent sequence uny uj ................................ u using Newton's 's method, which follows the formula:
w = ut -)((t))1
Each convergent sequence u ui ... nt, constitutes supervised data entry in the dataset. 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 summarized in Table
For the unsupervised data samples, only sample the perturbed initial states.
A.3 Implementations of loss functions
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 (u) the Jacobian matrix, du is the Newton step, and F(u) is the function defining the PDE.
Discretization Consider discrete representation of function a N N grid. The function and its perturbation du are represented as matrices:
ut 8u E RNNN
The function F(u), which involves both linear and nonlinear terms, is similarly represented as F(u) e N
Laplace Operator Regarding the representing the [u) with N x N grid function u, the discretized Laplace operator using finite difference method can be expressed convolution:
0
Au
-1
4
-1
*
-1
0
This convolution computes the result of the Laplace operator applied to the grid function w. 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.
15

---

A.4 Architecture of DeepONet
variant DeepONet used in our Newton-informed neural operator. In DeepONet, we introduce a hybrid architecture that combines convolutional layers with static trunk basis, optimized for grid-based data inputs common in computational applications like computational biology and materials science.
Branch Network 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 followed by additional Cony2D layers (128 filters, 5x5 kernel, stride and subse- quently 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 28]. and is dimensionally compatible with the output of the branch network.
During the forward computation, the branch network outputs are projected onto the trunk's s static basis via matrix multiplication, resulting in feature matrix that is reshaped into the grid dimensionality for output.
Forward Pass
Hyperparameters
The following table summarizes the key hyperparameters used in the DeepONet architecture:
Parameter
Value
Number of Conv2D layers
Filters in Conv2D layers
128, 128, 128, 128, 256
Kernel sizes in Cony2D layers
5x5, 3x3, 1x1, 5x5
Strides in Conv2D layers
2, 1, 1,2 2
Fully Connected Layer Sizes
256, branch features
GELU Table Hyperparameters of the DeepONet architecture
Activation Function
A.5
Training settings
Below we summarize the key configurations and parameters employed in the training for three cases:
Dataset
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 1: For method
Case 2: For method we use 5000 supervised data while for method 2, 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.
Optimization Technique
Optimizer: Adam, with learning rate of x 10-1 and a weight decay of 6
16

---

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 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 GPU- accelerated implementation of Newton's method, designed 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 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 optimize execution efficiency. contrast, the neural operator method is inherently parallelized, taking full advantage of the GPU architecture without the explicit use of multiple streams indicated in the table. The computational times of both methods were evaluated under 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 used for solving the sparse matrix systems. The following algorithm summmrissssssss the procedure to benchmark the time used for solving Newton's 's system for 5000 different initial states.
Algorithm 1 Solve Newton's Systems on GPU 1: procedure SOLVESYSTTEMSGPUAA u, rhs f,N)
> Precompute RHS and diagonal for all systems
2: rhs rhs_f f + u² Axu x diag + -2 u
3:
> Initialize solution storage
4:
della_n +- initialize zero matrix with shape of T
Solve each system
5:
for := to numm_sys do
6:
rhs rhs _isis
7:
diag t diag_list list
8:
A diagonal matrix with diag on the main diagonal
>> Solve the linear system
delta_ui uld spsolve rhs) 10: end for
9:
return transpose (destose(deela__.) (delta_))
12:
end procedure
B Supplemental material for proof
B.1 Preliminaries
Definition (Sobolev Spaces 6). Let Q2 be [0, 1]d and let D be the operator of the weak derivative
of single variable function and D D D ... D be the partial derivative where a
17

---

(01, a....... adl and D; is the derivative in the i-th variable. Let n E N and $ S cc. Then we
define Sobolev spaces
w (S2) = (1 (P()) : e LP(9) for all N with |a| S 5n}
with norm
1/p
Σ
D | fll
OS|alSo
i < co, and [I]ll maxos|alsn | 1
Furthermore, for f (Ji,...,Sa), 11.00 (2, (4) and only (L. for each i = 1.2...... and
[f == max il..d
(
When p 2, denote ( as (Q2) for E
Proposition 1 ((30), Suppose is is continuous non non-polynomial function and K is compact in TR.,, then there are positive integers constants Wks Çk for 1,...,P and bounded linear functionals Cg H (K) IR such that for any H (K),,
(u(o) (wy, E + (k)
cp cp r/pp [1/l [|]|)
(13)
(22(K)
Proposition 2 (122 Suppose a is a continuous non-polynomial function and Q is compact subset of Rd. For any Lipschitz function there exists a shallow neural network such that
Σ ago (wy x + bb)
S Cm -1/d 1/d
(14)
where C' depends on the Lipschitz constant but independent of m. Lemma ((22). The e-covering number of B, B]d, K()), satisfies
CB
K(e) &
for some constant C > 0, independent of B, and d.
Step 5: Now we estimate
B.2 Proof of Theorem 1
In this subsection, we present the proof of Theorem [ which describes the approximation ability of DeepONet.
Proof of Theorem Step 1: Firstly, we need to verify that the target operator (uu) is well-defined.
Due to Assumption (i), know that for X H²2),, Equation 2 will have unique solutions. This means that (u) is well-defined operator for the input space u E
Step 2: Secondly, we aim to verify that (uu) is Lipschitz-comtinuuuu operator in H for E X.
Consider the following:
((u+o)=f(c f vf
[(++v)= = ((u) + +
8,u(2) = (u(x) <
(15)
where ōu(zz) is the solution of Eq. input u + v. Denote
for the input u, and ),,((z)) is the solution of Eq. for the
8,u(x) 6u(x) e(x))
18

---

Therefore, we have:
[CC fu (<x)) =0
u)){(z) ())+())
Ao + ( x
(16)
x: E 0Ω.
Since u and are in 2 and 05 is in C (Assumption (ii)), according Theorem 4 in Section 6.3], there exist constants C aadd C such that:
(((=))) S Cl (E2) f (6)
S <C(lo(0))l
(17)
The last inequality is due to the boundedness of f( (" (Eu))u (Assumption (ii)).
Step 3: In the approximation,
first reduce the operator learning functional learning.
When the input function u(xr) belongs to X C H2² the output function also belongs to H provided that an is of class C?2 The function ((u) du can be approximated by two-layer network architected by the activation function (2x), which not a polynomial, the following form by Proposition
[30] (given in Subsection I6]:
((w)(a) -
cx cxlG(u)lo (wy a + (t)
S C.p [|(()))
(18)
(P(0)
where wy E Rd. Ço E for I,....p Cj continuous functional, and C3 constant independent of the parameters.
Denote (u) c& [Q(u),, which is Lipschitz-contimuous functional from (S2) to R. which is operator and C is a linear functional. The remaining task in
due to G is a approximation is to approximate these functionals by neural networks.
Step this step, we reduce the functional learning to function learning by applying the operator P as in Assumption (iv).
Based on (u) being Lipschitz- functional in H2)),, we have
(Øx(u) dx (Pu)) < |x Pull S Lgc,
where La is the Lipschitz constant of (u) for E X.
Furthermore, since Pu an re-dimensional term, i.e., can be denoted by the n-dimensional vector rewrite @ (Pu) vix. (Pu), where the R" RR for 1,..., Furthermore,
Pu E RR"
4/ is Lipschitz-comtinuuus function since dę is Lipschitz-continuuus and P is a continuous linear operator.
Step this step, we will approximate 4 using shallow neural networks.
Due to Proposition 2 we have that there is shallow neural network such that
(Pu)) Ago (M₂ Pu + bb) S Cm-
(19)
where a. E R" M2 E Rmxn and by E R". For the simplicity notations, we can replace My by an operator Wk.
P
Above all, we have that there neural network in Ep such that
Σ Ago (Wau+ by) (wy, x + (x) ((u)
+ +Cge+ (20)
(2((0)
where Ci is constant independent of m, e, and p. C2 is constant depended p.
II this paper, consistently employ the symbol C constant, which may vary from line line.
19

---

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 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 lies between two distinct solutions. Even a small perturbation of 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 approximation efforts within localized domains.
Conditions (ii) and (iii) serve to regularize the problem and ensure tractability. These conditions are indeed strai trfightforwwad to fulfill, contributing to the feasibility of the overall approach.
For the embedding operator P in (iv), there are 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 []6] achieve this task.
Here, we discuss more about the embedding operator P. One approach is to use the Argyris element 12 This method involves considering the 21 degrees of freedom shown in Fig. 7 In this figure, each denotes evaluation at 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.
Z3
122
n1
m
m
1
L2
L
m
L3
21
n3
Z2
Figure 7: Argyris method
Another alternative approach discretizing the input space is the bi-cubic Hermite finite
element method [23 18].
B.3
Proof of Theorem
The proof of Theorem inspired by that in 22
Proof of Theorem 2 Step 1: To begin with, we introduce a new term called the middle term, denoted
Esm (0), defined follows:
M Mr
Σ / (()))(())
Esm (0)
Ouy; (x)) dæ,
j=1 a
This term represents the limit case of (e) as the number of samples in the domain of the output space tends to infinity (M, + oo).
Then the error can be divided into two parts:
(0)s(0) (0) < [B()((0) Esm (0) [B(Esm (s)
(21)
20

---

Step 2: For [B(Esm (0) E [ this is the classical generalization error analysis, and the result
can be obtained from 34 42 We omit the details of this part, which can be expressed as
[B{ (0) - {sc( Cde Vog M
(22)
VM2
where C' is independent of the number of parameters dg and the sample size M.. In the following steps, we are going to estimate [E(Es(0) - Esm (f)) which is the error that comes from the sampling of the input space of the operator.
Step Denote
M S
Σ
/
Oluy: ()x)) dæ. (()))())
M
52
We first estimate the gap between Sta and Sa for any bounded parameters 0, 0'. Due to Assumption and (ii), we have that
S
1
M
L ((())(()) ((;; (z) Ouu;; ((z;)')(z) dæ
M
==1
1 M
Σ
(2g(u)) +O(u; 0(u; -O(,; - Ou,;; da
M 1=1
4
Σ
(u,) ||0 0'll p0
(23)
M
Step 4: Based on Step 3, we are going to estimate
E
sup [B]
SMM
E.S P
by covering the number of spaces.
Set 10....., 0k} is a E-covering of |- B, B3de i.e. for any 0 [-B, B]de there exists with
||0 1, Then we have
1/p
E
sup
SM 50
E
S P
Dc B, B,B
<TR
sup
S8
SM
S8,
SM
+
-0.
SM
E
[0,
SM
+ B
s.M
EE [ |)"]
1/p
0El B,
M
m
<I
max
SM 50
E
50;
+ M
8e
| (u,))
+bb (uz)
1.....K K
1/p
< V [u;)) ((u))
max
SM S0,
S.M
(24)
For 8cl // 1/p it can be approximate by
ScE 4 | 1/r
<
ScE
[r]
2r] 1/2p
E
|
1//2
1/20
8e||
|| ||
(d 1/P
For E max mmxx==l.... Il....
SM
E
S 00
by applied the result in
37
22
we know
1/p
16K
12
SM
S
max
M
21

---

Step 5: Now we estimate [8)) Esm (
Due to Assumption 2 and directly calculation, have that
[| | & C(1 +ypp)"
for constants C, > 0, depending only the measure μ and the constant C appearing in the upper bound 7 For example,
1 S (/x CLL+ (1+1llull ( pp dux
(/s exp [2pre (1 [11)l - al|u|| : e
+
<0
KP (1+
2np
e
$ C(1 yxxp)"
(25)
x
Based Lemma
we have that
1/p
CB
dopr
sup B3
S
E
< 160 + 2x 71DD)))
6+
M
DE[
for some constants C,7 independent of K. p. B, do, N, 0 and p 2. We now choose
c MM
so that
CB
dopr
(1+
(CB do (CBBMM))oor vp)
E+
Ma
VM,
Next, let p dg (BB(CBBV /w)) Then, (CBVM,) do/p vp exp
Log(CBvMMddd vp = ev do log(CB Log(CBVMM Ma).
p
and thus we conclude that
CB
dalp
Vp MM
(+ do lo (B))
E +
VM
On the other hand, we have
(1 (1+p))) = (1 yndo ( Mw)) 2x
Increasing the constant C' > 0, if necesssary, we can further estimate
(1 (1+ (1+7 log(CBB M.)) 2x
(1
< < do log(CB Mw))
2x 1/2
where C > 0 depends K, 7, /1 and the constant appearing in (7) but is independent of de, B and
N. We can express this dependence in the form C' C(1, V. +) > 0, as the constants si and ? depend on the Gaussian tail of μ and the upper bound on V. Tb Therefore,
[E(Es(0) (0) Esm (0) E E
sup
IS
IE ScMM
S CIE
sup 0E(-B, B)d
SMM.
EE
5g
S.M
DE
B,B
C
<
[+ Cdo log My) 2 200112)
(26)
VMa
22

---

NeurIPS Paper Checklist
Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The main claims made in the abstract and introduction accurately reflect our contributions and scope.
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
be expected to generalize other settings.
It is fine to include aspirational goals as motivation long it is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: Our paper discusses the limitations of our work in the section where establish and prove our network's approximation ability.
Guidelines:
The answer NA means that the paper has 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 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 sell-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 few datasets or with 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, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-tes system might not be used reliably to provide closed captions for online lectures because it fails 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 faimess..
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 impor- tant role in developing norms that preserve the integrity of the community, Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof? Answer: [Yes]
23

---

For each theoretical result, we provide the full set of assumptions and complete (and correct) proof.
Justification:
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: provide the code and dataset used the experiments as described in reproducing statement,
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary either make it possible for others replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of large language model), releasing of a model checkpoint, or other means that appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome deseribe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
24

---

Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
Justification:
We provide the link to the datasets.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips... cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for a new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
/ ccippi.cc guwdes/ddeeeee CooeeSumeeeeeeeeeeeeee ssi onyylicy) for more details.
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data,
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental Setting/Details
Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Question:
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 the core the paper to a level 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, as supplemental material.
7. Experiment Statistical Significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
We report statistical significance for key findings in our study. Since we do not claim our method is superior other models, our experiments are serving to justify our theoretical claims, therefore, we selectively provide statistical results where they are most relevant.
Justification:
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by emor bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
25

---

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 the standard error of the mean.
It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not show in tables figures symmetric error bars that would yield results that are out range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
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 compute required for each of the individual experimental runs as well 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 into the paper).
9. Code Of Ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics https: http://urrrpp.cppppeeeeeeeeeeeee cttp:////cbccceiiiccceceeellllll
Answer: [Yes]
Justification:
We read the NeurIPS Code of Ethics and confirm our research conform the
code.
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (c.g.. if there special consid- eration due to laws or regulations in their jurisdiction..
10. 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 obvious social impact.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
26

---

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 direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that 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 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 attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
No obvious such risks.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have 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 challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g.. code, data, models), used in the paper, properly credited and the license and terms of use explicitly mentioned and properly respected?
Question:
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 dataset.
The authors should state which version of the asset used and, possible, include 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.
27

---

If assets are released, the license, copyright information, and terms of use in the thcode. has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
package should be provided. For popular datasets, paperswi
latasets.paperssi
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
If this information not available online, the authors are encouraged to reach out to the asset's s creators.
13. 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/moode as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent 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.
14. Crowdsoureing 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 details about compensation (if any)?
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the paper doess not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion 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.
15. 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 equivalent approval/review based the requirements of your country institution) were obtained?
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Depending 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.
28

---

We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors adhere 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.
29


Reviews:
[{"review_id": "ggvokSOOCV", "reviewer": "Reviewer_6H5E", "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.\n\nClassical 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.\n\nThe 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.", "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? \n\n2. 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?\n\n3. 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.\n\n4. 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.\n\n5. Does the considered embedding operator using Argyris element satisfy Assumption 1.4? Could the paper provide a rigorous mathematical description?\n\n6. 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$?\n\n7. How does the paper obtain equation 25? Is there a typo?", "limitations": "Yes", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper introduces a novel methodology that combines classical Newton methods with neural network techniques to address the challenge of computing multiple solutions of nonlinear PDEs.\n\nThe paper provides complete theoretic results, including approximation and generalization error analysis.\n\nThe methodology presented in the paper has wide applicability in fields such as physics, biology, and engineering, where nonlinear PDEs with multiple solutions are common. By efficiently handling such problems, this approach can contribute to advancements in various scientific disciplines.", "weaknesses": "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.\n\n2. The symbols used in the text are somewhat confusing and require careful revision. \n2.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.\n2.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.\n2.3 What is the definition of $f(u)$ and $f'(u)$. Is $f$ just a function like $u^2$ rather than an operator?\n\n3. 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?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "## Weaknesses\n\n> Comparison:\n\n- The idea behind our proposed Newton Informed Operator Learning is rooted in the computational data generated by classical Newton's methods. When computing multiple solutions of nonlinear PDEs in pattern formation, the initial guesses can be in the millions, especially with fine grid solutions. In such cases, we utilize a subset of the data generated by the classical Newton's method to train the Newton Informed Operator. This training process is performed offline. Once the operator is well-trained, it can be applied to replace the classical solver, significantly speeding up the computation.\n\n- We agree that a more comprehensive comparison with existing methods could further strengthen our paper. However, our primary goal was to highlight the novel integration of classical Newton methods with neural network techniques for solving nonlinear PDEs. We have demonstrated the speedup by comparing it with traditional Newton methods. To the best of our knowledge, no existing methods specifically address the acceleration of nonlinear PDE computation in the way we propose.\n\n> Clarity of Notation:\n\n- 2.1 The use of symbols will be clarified in the revised manuscript to ensure consistency and clarity. Specifically, $u_*$ should be $\\delta u$ as in equation 2, i.e., $\\delta u := \\mathcal{G}(u)$, where $\\delta u$ is the solution of Eq. 2. The notation $u^*$ is the solution of nonlinear PDEs, namely, for an initial function $u_0$, assume the $n$-th iteration will approximate one solution, i.e., $(\\mathcal{G} + \\text{Id})^{(n)}(u_0) \\approx u^*$.\n\n- 2.2 The definition of the Sobolev space can be found in Appendix B.1.  We rephrase the assumption here. Assuming that $X$ has a Schauder basis $\\{b_n\\}$, we define the canonical projection operator $P_n$ based on this basis. The operator $P_n$ projects any element $u \\in X$ onto the finite-dimensional subspace spanned by the first $n$ basis elements $\\{b_1, b_2, \\ldots, b_n\\}$. Specifically, for $u \\in X$, $u = \\sum_{k=0}^\\infty \\alpha_k b_k$, let\n$$\nP_n(u) = \\sum_{k=0}^n \\alpha_k b_k,\n$$\nwhere $\\alpha_k$ are the coefficients in the expansion of $u$ with respect to the basis $\\{b_n\\}$.\nThe 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 \\in X$, $\\epsilon > 0$, there exists an $n$ such that\n$$\n\\|u - P_n u\\|_{H^2(\\Omega)} \\le \\epsilon, \\quad \\text{for all } u \\in X.\n$$\nThis ensures that $P_n$ effectively approximates elements of $X$ within the desired accuracy.\n\n- 2.3. $f(u)$ is a nonlinear function in $\\mathbb{R}$ and $f'(u)$ is the derivative with respect to $u$, $u:\\mathbb{R}^d \\to \\mathbb{R}$. For example, $f(u) := u^2$ and $f'(u)=2u$. \n\n > DeepONet:\n\n- We understand the concern regarding the distinction between our method and standard DeepONet. The key difference is our integration of a Newton-based loss function, which is specifically designed to address the unique challenges posed by nonlinear PDEs. While existing convergence rate analyses typically focus on elliptic equations and linear advection–diffusion–reaction equations, our analysis in this paper is centered on the newly derived loss function tailored for nonlinear PDEs.\n\n## Questions\n\n> Computational Efficiency:\n\n- For single linear Newton systems, our method is also significantly faster than traditional linear solvers. Furthermore, its primary advantage becomes more pronounced when solving multiple systems simultaneously. In this case, the learned operator greatly accelerates the solution process. This is because solving multiple nonlinear PDEs often involves computing many solutions. By generating solution data from various initial states, we can train the Newton operator on this diverse dataset. Efficiently sampling a wide range of initial states enhances the likelihood of converging to different solutions, which is crucial in complex solution landscapes with multiple basins of attraction.\n\n> Computational Time and Accuracy:\n\n- For the Gray-Scott model, we test both methods with the same set of initial conditions. Newton's method requires approximately 0.26 seconds to achieve an $L^2$ norm of the residual below $1 \\times 10^{-4}$, while the neural operator only requires 0.00028 seconds to achieve an $L^2$ norm of the residual around $1 \\times 10^{-1}$.\n\n> Approximation:\n\n- The approximation abilities discussed are based on specific assumptions and setups unique to our approach, particularly in handling multiple solutions. For Steps 3, 4, and 5, there is a proof for the generalization of DeepONet, which can be directly applied to other approximation proofs if the task is to use DeepONet to solve the problem. The proofs in Steps 1 and 2 aim to clarify the approximation in our setup and ensure that the operator we tend to approximate is well-defined.\n\n>Referencing:\n\n- Thanks for the careful reading. We cite Theorem 2.1 in the reference paper. In that paper, they consider function approximation, whereas our paper focuses on operator learning. If you consider the output of the operator learning, it should be a function, allowing us to use that theorem directly.\n\n> Clarifications:\n- The Argyris element, which satisfies $H^2$ approximation, will be utilized. The proof of the error associated with this element is shown in _S. Brenner, The Mathematical Theory of Finite Element Methods, Springer, 2008_. This choice is because this approximation includes second-order derivatives information in the approximation.\n\n> Theorem 2\n- $E$ is the expectation of sample points since we randomly choose sample points. The theorem holds for all $\\theta$ as we consider the uniform error. Non-uniform generalization error remains an open question in learning theory, even for function learning. \n\n> Eq.~(25)\n- Eq.~(25) is not a typo; it is based on the norm embedding, which we will discuss further in the paper.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "JfHjoa7Z6S", "reviewer": "Reviewer_tGuq", "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.", "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?\n- 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?\n- 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.\n- 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.\n- 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?\n- 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?  \n- How do results from method 1 and method 2 compare for the different use cases?", "limitations": "NA", "rating": 5, "confidence": 2, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "The method claims to learn multiple PDE solutions given any initial solution and uses a newton based loss for improved learning.", "weaknesses": "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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "## Weaknesses\n\n> Data Generation:\n- 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.\n\nHere is how we generate the supervised data samples:\n\n1. **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 **Gray-Scott model**.\n\n2. **Step 2**: The supervised dataset is generated by sampling a perturbation around the solution $u^i_{p} \\sim \\mathcal{N}(0, (-\\Delta)^{-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_0^i, u_1^i, \\ldots, u_n^i\\}$ using Newton's method, which follows the formula:\n   $$\n   u_{k+1}^i = u_k^i - J_f(u_k^i)^{-1} f(u_k^i).\n   $$\n   Each convergent sequence $\\{u_0^i, u_1^i, \\ldots, u_n^i\\}$ constitutes one supervised data entry in the dataset.\nIn 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.\n\n> Baseline Comparisons and Ablation Studies:\n- We acknowledge the reviewer's concern regarding the lack of baseline comparisons and ablation studies. However, the primary goal of our work is to establish that neural operators can effectively solve partial differential equations (PDEs) with multiple solutions. Our aim is to demonstrate the feasibility and potential of this approach rather than to provide an exhaustive comparison with existing neural operators.\n\n- To this end, we provided a general framework for finding multiple solutions using neural operators and used the DeepPOD as an example to showcase the applicability of our method to several PDE problems with multiple solutions. Our focus was on illustrating that various neural operators, such as the Fourier Neural Operator, can be seamlessly integrated into our framework, replacing the DeepOnet/DeepPOD. This demonstrates the versatility and robustness of our approach.\n\n\n\n## Questions\n\n>Data Generation and Testing Conditions:\n- We have addressed the procedures in the weakness regarding how to generate 500 supervised samples. As for the test data and training data, their initial states are both sampled from the same distribution $u^i_{p} \\sim \\mathcal{N}(0, (-\\Delta)^{-3})$.\n\n- If the initial state is only slightly perturbed, the Newton informed neural operator will converge quickly as the classical Newton's method, since the Newton informed neural operator approximates each Newton's step and recursively predicts the true solution. We compare the error for each step with Newton's method in Figure 1 and the error for the final prediction in Figure 3.\n\n> Error Metrics and Validation:\n- The error is measured by $L^2$ and $H^1$ norm, for example in Figure 1 and Figure 3. The definitions are provided in Appendix B.1.\n\n>The Information Presented in the Table is Unclear:\n- The time reported in Table 1 reflects the efficiency of the neural operator compared with classical direct solvers for Newton's linear systems. It is not because the neural operator converges faster in the sense of fewer iterations. Rather, the time is measured for solving 500 or 5000 independent Newton's linear systems, each with different initial states.\n- This efficiency gain is due to the neural operator's ability to quickly approximate the solution to each linear system, leveraging the pre-trained network to perform these tasks much faster than traditional solvers. The details of how we benchmarked the efficiency are provided in the paragraph **Efficiency** and in Appendix A.5 **Benchmarking Newton’s Method and Neural Operator-Based Method**. These sections explain the experimental setup and the metrics used to evaluate the performance.\n\n> Does the Neural Operator Recursively Predict the Solution? How does the Method Scale with Grid Size.\n- Neural operator recursively predict the solution, since the Newton informed neural operator approximates the Newton's iteration step at given state. Since the architecture is simple, the memory is only of  $\\mathcal{O}(N)$  in terms of the degrees of freedom $N$.\n- At the evaluation phase of the Newton informed neural operator-based solver, the computational cost is $\\mathcal{O}(N)$ while the classical direct solver for solving the Newton's step is $\\mathcal{O}(N^3)$ for dense systems and less for sparse systems, but still significant. \n- The POD calculation does indeed scale with the grid size, and for grids with millions of cells, this could lead to increased computational costs. In such cases, alternative approaches, such as using a pre-defined mesh-less basis like kernel basis functions, can be utilized to manage computational complexity and maintain efficiency.\n\n>  DeepPOD Network Architecture?\n- Yes, the deepPOD network in Figure 2 refers to the neural operator where the trunk net is replaced by pre-calculated basis functions obtained via the Proper Orthogonal Decomposition (POD) method. The discussion regarding the architecture is discussed in Appendix A.3.\n\nWe want to emphasize that our framework is not tied to a specific neural operator architecture. While we demonstrated our method using the deepPOD as an example, our primary claim is that neural operators, in general, can effectively find multiple solutions for partial differential equations. For different PDEs or larger grid sizes, more advanced neural operator architectures could be employed to achieve optimal performance.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "iS63W4VHI4", "reviewer": "Reviewer_3WGd", "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.", "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.\n2. Line 37 page 2, \"Newton methods provide well-defined locally\": is there a missing word?\n3. Line 42 page 2, \"Our approach combines the classical Newton method, which...\": combines it with what?\n4. Lines 92-95, page 3: I'm struggling to understand what the authors are saying by \"remains uncontrollable\". Reword for clarity?\n5. Lines 99-103, page 3: Is this true? If a solution exists, there must be an operator, surely? Citations and/or proof needed.\n6. 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. \n7. 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.\n8. 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.\n9. Page 4, remark 1: DeepONet hasn't even been described yet. How are general readers supposed to know how this fits into this paper?\n10. 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.\n11. 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": "The authors do not explicitly have a limitations section, but they briefly address their limitations in their conclusion.", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 4, "strengths": "The fundamental idea behind this work is timely, given the rise of ML-based techniques to solve PDEs, including but not limited to physics-informed neural nets (PINNs) and operator learning. ML is a promising way of obtaining multiple solutions to PDEs, and the Newton informed operator appears to do this quite well.", "weaknesses": "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).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "## Weaknesses\n\n> Clarity of Exposition:\n- We will include a clear architectural diagram of the Newton-informed Neural Operator and provide a step-by-step derivation of the Newton loss function as follows:\n\nFor the nonlinear PDEs:\n$$\n\\begin{cases}\n\\mathcal{L} u(\\mathbf{x}) = f(u), & \\mathbf{x} \\in \\Omega \\\\\nu(\\mathbf{x}) = 0, & \\mathbf{x} \\in \\partial \\Omega,\n\\end{cases}\n$$\nwe can represent the PDEs as $ F(u) = 0 $, where $ F $ is an operator defined as $ F(u) = \\mathcal{L}(u) - f(u) $. In the context of the Newton method for operators, we start with an initial condition $ u_0 $ and iteratively obtain a solution. For each iteration, $ u_{n+1} = u_n + \\delta u_n $, where $ u_n $ is the solution from the previous iteration, and $ \\delta u_n $ is obtained from the iteration equation $ F'(u) \\delta u = F(u) $. Here, $ F'(u) $ is the (Fréchet) derivative of the operator, defined as follows:\n\nTo find $ F'(u) $ such that for any $ v \\in H^2(\\Omega) $,\n$$\n\\lim_{\\|v\\| \\to 0} \\frac{\\|F(u + v) - F(u) - F'(u)v\\|}{\\|v\\|} = 0,\n$$\nwhere $ \\|\\cdot\\| $ is the norm in the Sobolev space $ H^2(\\Omega) $. In our case,\n$$\nF'(u)v + \\mathcal{O}(v) := \\mathcal{L}(u + v) - \\mathcal{L}(u) - (f(u + v) - f(u)) = \\mathcal{L}(v) - f'(u)v + \\mathcal{O}(v).\n$$\n\nTherefore, the Newton linear PDE we use in this paper is given by:\n$$\n\\begin{cases}\n(\\mathcal{L} - f'(u)) \\delta u(\\mathbf{x}) = \\mathcal{L} u - f(u), & \\mathbf{x} \\in \\Omega \\\\\n\\delta u(\\mathbf{x}) = 0, & \\mathbf{x} \\in \\partial \\Omega.\n\\end{cases}\n$$\nThis represents the linear PDE for $ \\delta u $. Based on the assumptions in our paper, we can prove that this has a unique solution.\n\nIn this paper, we aim to solve nonlinear PDEs with any given initial guess. Therefore, we need to learn the Newton solver, which is an operator between $ u $ and $ \\delta u $. The way to learn this is by employing DeepONet with a Newton information loss. The structure of DeepONet is provided in the appendix, and the Newton information loss is derived from the above equation. The Newton loss function is defined as:\n$\nE_{N}(\\theta) := \\frac{1}{N_u \\cdot N_x} \\sum_{j=1}^{N_u} \\sum_{k=1}^{N_x} \\left|(\\mathcal{L} - f'(u_j(x_k))) \\mathcal{O}(u_j; \\theta)( x_k) - \\mathcal{L} u_j(x_k) - f(u_j(x_k))\\right|^2\n$\nwhere $ u_1, u_2, \\ldots, u_{N_u} \\sim \\nu $ are independently and identically distributed (i.i.d.) samples in $ \\mathcal{X} $, and $ x_1, x_2, \\ldots, x_{N_x} $ are uniformly i.i.d. samples in $ \\Omega $. As you can see, this loss function directly follows from the Newton iteration.\n\n> Comparison with HomPINNs:\nHomPINNs are developed to compute multiple solutions of nonlinear PDEs. In this paper, we present an operator learning approach combined with Newton's method to learn the nonlinear solver. While this approach is not specifically designed for computing multiple solutions, it can accelerate the nonlinear solver process if initial guesses are provided.\n\n## Questions\n\n> Multiple Solutions and Nonlinearity:\n- In this paper, we focus on the multiple solutions arising from nonlinear terms in nonlinear PDEs. Analyzing the energy formula reveals that such terms may contain multiple local minima, leading to multiple solutions, commonly referred to as pattern formation. In contrast, for linear PDEs with multiple solutions, the differences are typically in the form of a constant $C$ or other simple polynomial functions. These simpler cases are not the focus of our study.\n\n> Methodological Clarifications:\n- Thank you for the careful reading. We will correct the typos in the paper and the details of our method as shown above. Additionally, we will move the structure of DeepONet from the appendix into the introduction.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "aziOpwoF0V", "reviewer": "Reviewer_zTaB", "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.", "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?\n\n2) 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?\n\n3) 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": "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.", "rating": 6, "confidence": 4, "soundness": 2, "presentation": 1, "contribution": 3, "strengths": "In general, the topic addressed in the manuscript is very interesting. Many machine learning techniques rely on the assumption that (at least...) there exists a unique solution to a problem (like a PDE here), and then approximate it as accurately as possible. The authors consider a much more challenging case: what if the problem has multiple solutions? They then proceed to use a classical numerical scheme (Newton's method) that is typically used to obtain one of the solutions of the given PDE, and approximate its action through a neural operator. This is interesting even away from the particular problem addressed here, because Newton's method is being used in many other applications as well.", "weaknesses": "Major:\n1) 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.\n\n2) 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.\n\n3) The approximation of Newton's method has been addressed in earlier work already, which is not cited.\n\n [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\n\n [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.\n\n\n\nMinor:\n1) 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.\n2) l.37: \"Newton methods provide well-defined locally\" is missing a word\n3) l.50: \"solve the inverse problem\" is not clear at this point. Which inverse problem? \n4) 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.\n5) Section 2.1 is not necessary, the information can be incorporated more concisely in the introduction.\n6) l106,107: the index \"i\" is not used in the equation, same with \"u_0\" and \"u dash\".\n7) l140 P_dash is not explained.\n8) l172: this is not a complete sentence.\n9) 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).\n10) l246-l248 contain multiple mistakes in the sentence structure.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "## Weaknesses\n\n> Approximation of Newton's Method and Finding Multiple Solutions\n\n- Thank you for your insightful comments.\nWe acknowledge that a large portion of machine learning research indeed focuses on learning distributions, as seen in generative models like diffusion models, score-based learning, and variational autoencoders. These models typically aim to map an initial distribution to a final distribution, often in a high-dimensional space like images.\n\n- However, we would like to clarify that our method, the Newton Informed Neural Operator, differs from these generative modeling approaches. Our approach is not designed to learn or output a distribution of solutions per se. Instead, it focuses on efficiently approximating the deterministic action of Newton's method on a given initial state, specifically in the context of solving nonlinear PDEs.\n\n- Our primary goal is to enhance the computational efficiency of finding multiple solutions by allowing for extensive exploration of initial conditions. The primary advantage of our approach lies in its efficiency; This efficiency gain is due to the neural operator's ability to quickly approximate the solution to each linear system, leveraging the pre-trained network to perform these tasks much faster than traditional solvers. The details of how we benchmarked the efficiency are provided in the paragraph **Efficiency** and in Appendix A.5 **Benchmarking Newton’s Method and Neural Operator-Based Method**. This significant speed-up allows for the practical exploration of a much larger space of initial conditions in a reasonable timeframe. By efficiently sampling a wide variety of initial states, our method increases the likelihood of finding multiple solutions of a PDE, especially in cases where the solution landscape is complex. \n\n>  Prior Work on Newton's Method Approximation\n\nWe appreciate the reviewer's insight into prior work that also addresses the approximation of Newton's method using neural networks. We will ensure to include the following references in the revised manuscript:\n\n- While these cited works focus on ODEs and propose neural network architectures like R2N2 to emulate iterative methods such as Runge-Kutta and Krylov, our paper extends these concepts to handle partial differential equations (PDEs). This extension involves significant additional considerations including spatial discretization, boundary conditions, and more complex solution structures.\n\n- Our contribution specifically introduces the **Newton Informed Neural Operator**, which incorporates a **Newton informed loss function**. This function leverages classical Newton methods to refine the approximation capabilities of neural networks for PDEs. This approach is distinct in that it is tailored to capture multiple solutions and managing spatial complexity.\n\n- We believe that while the foundational ideas may overlap, our work provides a novel application and significant advancements in the context of PDEs. We will clearly delineate these differences in our revision and appropriately contextualize our contributions in relation to the existing literature.\n\n> Typos and Expressions\n\n Thank you for your careful reading. We will correct the typos and inappropriate expressions in the manuscript. Specifically:\n  - In line 172, the term $P_{dash}$ will be replaced with a rigorous expression. Please also refer to our response to Reviewer 6H5E for additional details.\n  - In line 244, $L(u)$ should be corrected to $-\\Delta u$.\n\n## Questions\n\nWe believe all three questions are discussing the rationale for training a neural operator to approximate Newton's method. We list it in the following:\n\n- There exists no general method with a guarantee that all the solutions for complex nonlinear problems can be found, for example, the Gray Scott model. Newton's method is a very efficient method but suffers from the computational cost, especially when extensively sampling from various initial conditions as we discussed in the weaknesses section.\n\n- Once the neural operator is trained, this efficiency gain is due to the neural operator's ability to quickly approximate the solution to each linear system, leveraging the pre-trained network to perform these tasks much faster than traditional solvers. The details of how we benchmarked the efficiency are provided in the paragraph **Efficiency** and in Appendix A.5 **Benchmarking Newton’s Method and Neural Operator-Based Method**.\n\n- For a simple problem such as the convex problem in our **Case 1**, the extra time for data generation and training makes the neural operator-based method less useful, but for a complex case like the Gray Scott model, the efficiency gain of a pre-trained neural operator compensates for the cost of training.\n\n## Limitations and Broader Impact:\nWe admit that limitations are not discussed enough. The limitations include two main parts:\n  - We did not conduct baseline comparisons or ablation studies, as our primary intent was to prove the concept that neural operators can be effectively used to solve PDEs with multiple solutions. Future work could explore these aspects in greater detail, comparing different neural operator architectures and their performance relative to traditional numerical methods.\n  - Although the Newton informed loss function alleviates the cost for data generation, the data efficiency of the method still needs to be improved.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
We thank all the reviewers for their careful reading and insightful comments. We would like to clarify our motivation and summarize the pipeline of our method as follows:

**Motivation** 

The Newton Informed Neural Operator, focuses on efficiently approximating the iteration of Newton's method for a given state, specifically in the context of solving nonlinear PDEs with **Pattern Formation** (isolated multiple solutions). The primary approach of our method is to enhance the computational efficiency of finding multiple solutions by enabling extensive exploration of initial conditions. Let us briefly explain the motivation and methodology of our approach.

Newton's method for finding the roots of a function $ f(x) = 0 $ iteratively updates an initial guess $ x_0 $ using the formula:
$$
x_{k+1} = x_k - J_f(x_k)^{-1} f(x_k),
$$
where $ J_f(x_k) $ is the Jacobian matrix of $ f $ at $ x_k $. This method converges to a solution $ x^* $ if the initial guess $ x_0 $ is sufficiently close to $ x^* $.

In the context of nonlinear partial differential equations (PDEs), solving the nonlinear system often requires multiple Newton iterations, which can be computationally expensive, particularly when computing multiple solutions. Traditional numerical solvers for Newton's method can become a bottleneck due to their high computational cost. Our Newton Informed Neural Operator addresses this challenge by learning the operator defined within Newton's iterations, thereby speeding up the process and improving efficiency.

Our method approximates one iteration of Newton's method on a given state using a neural operator. The neural operator can perform these Newton's iterations orders of magnitude faster than using traditional solvers for Newton's iteration, as demonstrated in Table 1 of the manuscript. This significant speed-up allows for the practical exploration of a much larger space of initial conditions in a reasonable timeframe.

By efficiently sampling a wide variety of initial states, our method increases the likelihood of converging to different solutions of nonlinear PDEs, especially in cases where the solution landscape is complex and includes multiple basins of attraction. 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.

**Pipeline of the method**

To illustrate the logic of our algorithm, consider the following steps:

1. **Initial Sampling**: Generate a diverse set of initial conditions $ \{x_0^i\} $ to explore multiple solutions.
2. **Neural Network Approximation**: Use the trained neural operator $ G_f $ to approximate the Newton step for each initial condition:
   $$
   x_{k+1}^i = x_k^i - G_f(f(x_k^i)),
   $$
   where the trained neural operator $ G_f $ approximates the action of $ J_f(x_k)^{-1} $.
3. **Iteration**: Repeat the neural operator-based Newton steps until convergence for each initial condition.
4. **Solution Identification**: Collect and analyze the converged solutions to identify distinct solutions of the PDE.

We will expand on this explanation in the revised manuscript to make the benefits and rationale of our approach more explicit. Thank you for highlighting this important aspect, which allows us to refine the presentation of our work.

Meta-review: 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.
Decision: Accept (poster)