Newton Informed Neural Operator for Solving Nonlinear Partial Differential Equations

Weaknesses

Clarity of Exposition:

• 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:

For the nonlinear PDEs:

$$\begin{cases} \mathcal{L} u(\mathbf{x}) = f(u), & \mathbf{x} \in \Omega \\ u(\mathbf{x}) = 0, & \mathbf{x} \in \partial \Omega, \end{cases}$$

we 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:

To find $F'(u)$ such that for any $v \in H^2(\Omega)$,

$$\lim_{\|v\| \to 0} \frac{\|F(u + v) - F(u) - F'(u)v\|}{\|v\|} = 0,$$

where $\|\cdot\|$ is the norm in the Sobolev space $H^2(\Omega)$. In our case,

$$F'(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).$$

Therefore, the Newton linear PDE we use in this paper is given by:

$$\begin{cases} (\mathcal{L} - f'(u)) \delta u(\mathbf{x}) = \mathcal{L} u - f(u), & \mathbf{x} \in \Omega \\ \delta u(\mathbf{x}) = 0, & \mathbf{x} \in \partial \Omega. \end{cases}$$

This represents the linear PDE for $\delta u$. Based on the assumptions in our paper, we can prove that this has a unique solution.

In 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: $E_{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$ where $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.

Comparison with HomPINNs: HomPINNs 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.

Questions

Multiple Solutions and Nonlinearity:

• 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.

Methodological Clarifications:

• 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.

We would like to further clarify the novelty of our method:

1. We introduced a general computational framework (as illustrated in the global response) for learning the Newton's nonlinear solver using neural operators integrated with the Newton-informed loss, with DeepONet/DeepPOD (architecture detailed in Appendix A.3) as a representative example. This computational framework can be used to identify multiple solutions of nonlinear PDEs.

2. After a thorough review of the HomPINNs paper, we identified a key distinction between our method and theirs. In HomPINNs, networks are primarily used to approximate/parametrize the solution. In contrast, our method leverages networks to approximate the Newton solver—specifically, the mapping from $u_k$ to $\delta u$, $u_{k+1} = u_k+ \delta u$. While HomPINNs are designed to compute multiple solutions starting from simple initial functions, typically on coarse grids, by combining the homotopy approach with PINNs to learn multiple solutions on finer grids, our Newton-informed operator learning method directly trains the nonlinear solver on fine grids. This approach offers greater efficiency and reduces the risk of overlooking solutions that might be missed by HomPINNs.

Weaknesses

Approximation of Newton's Method and Finding Multiple Solutions

• Thank you for your insightful comments. We 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.

• 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.

• 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.

Prior Work on Newton's Method Approximation

We 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:

• 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.

• 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.

• 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.

Typos and Expressions

Thank you for your careful reading. We will correct the typos and inappropriate expressions in the manuscript. Specifically:

• 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.
• In line 244, $L(u)$ should be corrected to $-\Delta u$.

Questions

We believe all three questions are discussing the rationale for training a neural operator to approximate Newton's method. We list it in the following:

• 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.

• 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.

• 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.

Limitations and Broader Impact:

We admit that limitations are not discussed enough. The limitations include two main parts:

• 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.
• Although the Newton informed loss function alleviates the cost for data generation, the data efficiency of the method still needs to be improved.

Weaknesses

Data Generation:

• 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.

Here is how we generate the supervised data samples:

1. Step 1: We use a classical numerical solver to obtain single (multiple) solutions of nonlinear PDEs. For example, in Case 2, there exist four solutions $u^1, u^2, u^3, u^4$. We have one solution for Case 1 and 10 solutions for the Gray-Scott model.

2. 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:

   $$u_{k+1}^i = u_k^i - J_f(u_k^i)^{-1} f(u_k^i).$$

   Each convergent sequence $\{u_0^i, u_1^i, \ldots, u_n^i\}$ constitutes one supervised data entry in the dataset. In this case, we consider the initial conditions as perturbed, with all perturbations applied around the true solution. A comparison between the traditional method and our proposed method is summarized in Table 1.

Baseline Comparisons and Ablation Studies:

• 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.

• 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.

Questions

Data Generation and Testing Conditions:

• 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})$.

• 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.

Error Metrics and Validation:

• 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.

The Information Presented in the Table is Unclear:

• 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.
• 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.

Does the Neural Operator Recursively Predict the Solution? How does the Method Scale with Grid Size.

• 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$.
• 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.
• 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.

DeepPOD Network Architecture?

• 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.

We 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.

Dear Reviewer,

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

Best wishes, Authors

Weaknesses

Comparison:

• 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.

• 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.

Clarity of Notation:

• 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^*$.

• 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

$$P_n(u) = \sum_{k=0}^n \alpha_k b_k,$$

where $\alpha_k$ are the coefficients in the expansion of $u$ with respect to the basis $\{b_n\}$. The canonical projection operator $P_n$ is a linear bounded operator on $X$. According to the properties of Schauder bases, these projections $P_n$ are uniformly bounded by some constant $C$. Furthermore, we assume, for any $u \in X$, $\epsilon > 0$, there exists an $n$ such that

$$\|u - P_n u\|_{H^2(\Omega)} \le \epsilon, \quad \text{for all } u \in X.$$

This ensures that $P_n$ effectively approximates elements of $X$ within the desired accuracy.

• 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$.

DeepONet:

• 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.

Questions

Computational Efficiency:

• 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.

Computational Time and Accuracy:

• 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}$.

Approximation:

• 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.

Referencing:

• 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.

Clarifications:

• 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.

Theorem 2

• $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.

Eq.~(25)

• Eq.~(25) is not a typo; it is based on the norm embedding, which we will discuss further in the paper.