Derivative-enhanced Deep Operator Network

Authors: yuan  qiu
nolan  bridges
peng  chen
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 06:54:16

Abstract: The deep operator networks (DeepONet), a class of neural operators that learn mappings between function spaces, have recently been developed as surrogate models for parametric partial differential equations (PDEs). In this work we propose a derivative-enhanced deep operator network (DE-DeepONet), which leverages derivative information to enhance the solution prediction accuracy and provides a more accurate approximation of solution-to-parameter derivatives, especially when training data are limited. DE-DeepONet explicitly incorporates linear dimension reduction of high dimensional parameter input into DeepONet to reduce training cost and adds derivative loss in the loss function to reduce the number of required parameter-solution pairs. We further demonstrate that the use of derivative loss can be extended to enhance other neural operators, such as the Fourier neural operator (FNO). Numerical experiments validate the effectiveness of our approach.

Main paper: 
 introduction: 
 Using neural networks to approximate the maps between functions spaces governed by parametric PDEs can be very beneficial in solving many-query problems, typically arising from Bayesian inference, optimization under uncertainty, and Bayesian optimal experimental design. Indeed, once pre-trained on a dataset, neural networks are extremely fast to evaluate given unseen inputs, compared to traditional numerical methods like the finite element method. Recently various neural operators are proposed to enhance the learning capacity, with two prominent examples deep operator network (DeepONet) [1] and Fourier neural operator (FNO) [2] , which are shown to be inclusive of each other in their more general settings [3, 4] , see also their variants and other related neural operator architectures in [5, 6, 7, 8, 9, 10] . Though these work demonstrate to be successful in approximating the output function, they do not necessarily provide accurate approximation of the derivative of the output with respect to the input, which are often needed for many downstream tasks such as PDEconstrained optimization problems for control, inference, and experimental design [11, 12, 13, 14] .In this paper, we propose to enhance the performance of DeepONet through derivative-based dimension reduction for the function input inspired by [15, 16, 17, 18] and the incorporation of derivative information in the training to learn both the output and its derivative with respect to the input inspired by [19, 20] . These two derivative-enhanced approaches can significantly improve DeepONet's approximation accuracy for the output function and its directional derivative with respect to the input function, especially when the training samples are limited. We provide details on the computation of derivative labels of the solution of PDEs in a general form as well as the derivative-based dimension reduction to largely reduce the computational cost. We demonstrate the effectiveness of our proposed method (DE-DeepONet) compared to three other neural operators, including DeepONet, FNO, and derivative-informed neural operator (DINO) [20] , in terms of both test errors and computational cost. In addition, we apply derivative learning to train the FNO and also compare its performance with other methods. The code for data generation, model training and inference, as well as configurations to reproduce the results in this paper can be found at https://github.com/qy849/DE-DeepONet. This section presents the problem setup, high-fidelity approximation using finite element for finite dimensional discretization, and the DeepONet architecture in learning the solution operator. 
problem setup: 
 Let Ω ⊂ R d denote an open and bounded domain with boundary ∂Ω ⊂ R d-1 , where the dimension d = 1, 2, 3. We consider a PDE of the general form defined in Ω as R(m, u) = 0, (1) prescribed with proper boundary conditions. Here m ∈ V in is an input parameter function defined in a separable Banach space V in with probability measure ν and u ∈ V out is the output as the solution of the PDE defined in a separable Banach space V out . Our goal is to construct a parametric model û(m; θ) to approximate the solution operator that maps the parameter m to the solution u.Once constructed, the parametric model û(m; θ) should be much more computationally efficient to evaluate compared to solving the PDE with high fidelity approximation.For the high fidelity approximation of the solution, we consider using a finite element method [21] in this work. We partition the domain Ω into a finite set of subregions, called cells or elements. Collectively, these cells form a mesh of the domain Ω. Let h represent the diameter of the largest cell. We denote V in h indexed by h as the finite element space for the approximation of the input space V in with Lagrange basis {ϕ in 1 ,• • • , ϕ in N in h} of dimension N in h such that ϕ in i (x (j) ) = δ ij at the finite element node x (j) , with δ ij being the Kronecker delta function. Similarly, we denote V out h as the finite element space for the approximation of the solution space V out with basis{ϕ out 1 , • • • , ϕ out N out h } of dimension N out h .Note that for the approximation to be high fidelity, N in h and N out h are often very large. To this end, we can write the high fidelity approximation of the input and output functions asm h (x) = N in h i=1 m i ϕ in i (x) and u h (x) = N out h i=1 u i ϕ out i (x), with coefficient vectors m = (m 1 , • • • , m N in h ) T ∈ R N in h and u = (u 1 , • • • , u N out h ) T ∈ R N out h, whose entries are the nodal values of m and u at the corresponding nodes.We briefly review the DeepONet architecture [1] with a focus on learning the solution operator of the PDE in Equation (1) . To predict the evaluation of solution function u at any point x ∈ Ω ∪ ∂Ω for any given input function m, [1] design a network architecture that comprises two separate neural networks: a trunk net t(• ; θ t ), which takes the coordinate values of the point x at which we want to evaluate the solution function as inputs, and a branch net b(• ; θ b ), which takes the vector m encoding the parameter function m as inputs. In [1] , the vector m is the function evaluations at a finite set of fixed points {x (j) } N in h j=1 ⊆ Ω ∪ ∂Ω, that is, m = (m(x (1) ), • • • , m(x (N in h ) )) T , which corresponds to coefficient vector in the finite element approximation with Lagrange basis at the same nodes. If the solution function is scalar-valued, then both neural networks output vectors of the same dimensions. The prediction is obtained by taking the standard inner product between these two vectors and (optionally) adding a real-valued bias parameter, i.e, û(m; θ)(x) = ⟨b(m; θ b ), t(x; θ t )⟩ + θ bias , with θ = (θ b , θ t , θ bias ). If the solution u is vector-valued of N u components, i.e., u = (u (1) , . . . , u (Nu) ), as in our experiments, we can use one of the four approaches in [3] to construct the DeepONet. Specifically, for each solution component, we use the same outputs of branch net with dimension N b and different corresponding groups of outputs of trunk net to compute the inner product. More precisely, the solution u (i) of component i = 1, . . . , N u , is approximated byEQUATIONwheret (i) (x; θ t ) = t(x; θ t )[(i -1)N b + 1 : iN b ],the vector slice of t(x; θ t ) with indices ranging from (i -1)N b + 1 to iN b , i = 1, . . . , N u .For this construction, the outputs of the branch net can be interpreted as the coefficients of the basis learned through the trunk net. By partitioning the outputs of the trunk net into different groups, we essentially partition the basis functions used for predicting different components of the solution.The DeepONet is trained using dataset D = {(m (i) , u (i) )} N D i=1 with N D samples, where m (i) are random samples independently drawn from ν and u (i) are the solution of the PDE (with slight abuse of notation from the vector-valued solution) at m (i) , i = 1, . . . , N D . 
de-deeponet: 
 The DeepONet uses the input parameter and output solution pairs as the labels for the model training. The approximation for the derivative of the solution with respect to the parameter (and the coordinate) are not necessarily accurate. However, in many downstream tasks such as Bayesian inference and experimental design, the derivatives are required. We consider incorporating the the Fréchet derivative du(m; •) for the supervised training of the DeepONet, which we call derivative-enhanced DeepONet (DE-DeepONet). By doing this we hope the optimization process can improve the neural network's ability to predict the derivative of the output function with respect to the input function. Let θ denote the trainable parameters in the DE-DeepONet. We propose to minimize the loss functionEQUATIONwhere the L 2 (Ω) norm of a square integrable function f is defined as ∥f ∥ 2L 2 (Ω) = ( Ω ∥f (x)∥ 2 2 dx) 1/2, the Hilbert-Schmidt norm of an operator T : H → H that acts on a Hilbert space H is defined as ∥T ∥ 2 HS = i∈I ∥T e i ∥ 2 H , where {e i } i∈I is an orthonormal basis of H. Here, λ 1 , λ 2 > 0 are hyperparameters that balance different loss terms.The main challenge of minimizing the loss function (3) in practice is that with high fidelity approximation of the functions m and u using high dimensional vectors m and u, the term∥du(m; •) -dû(m; •)∥ 2 HS (approximately) becomes ∥∇ m u -∇ m û∥ 2 F, where the Frobenius norm of a matrix M ∈ R m×n is defined as∥M ∥ 2 F = m i=1 n j=1 M 2 ij .It is a critical challenge to both compute and store the Jacobian matrix at each sample as well as to load and use it for training since it is often very large with sizeN out h × N in h , with N in h , N out h ≫ 1.To tackle this challenge, we employ dimension reduction for the high dimensional input vector m. The reduced representation of m is given by projecting m into a low dimensional linear space spanned by a basis ψ in 1 , . . . ,ψ in r ∈ R N in h , with r ≪ N in h , that is, m = (⟨ψ in 1 , m⟩, . . . , ⟨ψ in r , m⟩) T ∈ R r .To better leverage the information of both the input probability measure ν and the map u : m → u(m), we consider the basis generated by active subspace method (ASM) using derivatives [22] . ASM identifies directions in the input space that significantly affect the output variance, or in which the output is most sensitive. See Section 3.3 for the detailed construction. In this case, the termE m∼ν ∥du(m; •) - dû(m; •)∥ 2HS can be approximated by1 N1N2 N1 i=1 N2 j=1 ∥∇ m u(m (i) )(x (j) ) -∇ m û(m (i) )(x (j) )∥ 2 2with a small amount of functions m (i) sampled from input probability measure ν and points x (j) in the domain Ω. Note that ∇ m û(m (i) )(x (j) ) is vector of size r, which is computationally feasible. For comparison, we also conduct experiments if the basis is the most commonly used Karhunen-Loève Expansion (KLE) basis.We next formally introduces our DE-DeepONet, which uses dimension reduction for the input of the branch net and incorporates its output-input directional derivative labels as additional soft constraints into the loss function.We incorporate dimension reduction into DeepONet to construct a parametric model for approximating the solution operator. Specifically, if the solution is scalar-valued (the vector-valued case can be constructed similar to (2)), the prediction is given byEQUATIONwhere θ = (θ b , θ t , θ bias ) are the parameters to be learned. The branch net b(•; θ b ) and trunk net t(•; θ t ) can be chosen as an MLP, ResNet, etc. Note that the branch net takes a small vector of the projected parameter as input. We also apply the Fourier feature embeddings [5] , defined as γ(x) = [cos(Bx), sin(Bx)], to the trunk net, where each entry in B ∈ R m×d is sampled from a Gaussian distribution N (0, σ 2 ) and m ∈ N + , σ ∈ R + are hyper-parameters.In practical training of the DE-DeepONet, we formulate the loss function as followsEQUATIONwhere {x (j) } Nx j=1 are the nodes of the mesh,Ψ in = (ψ in 1 | • • • |ψ in r )is the matrix collecting the nodal values of the reduced basis of the input function space, and Φ out (x) = (ϕ out 1 (x), . . . , ϕ outN out h (x))is the vector-valued function consisting of the finite element basis functions of the output function space. The err({(a (i) , b (i) )} n i=1 ) denotes the relative error between any two groups of vectors a (i) and b (i) , computed aserr({(a (i) , b (i) )} n i=1 ) = ( n i=1 ∥a (i) -b (i) ∥ 2 2 ) 1/2 ε + ( n i=1 ∥b (i) ∥ 2 2 ) 1/2, where ε > 0 is some small positive number to prevent the fraction dividing by zero.In the following, we explain how to compute different loss terms in Equation (5) • The first term is for matching the prediction of the parametric model û(m (i) ; θ) evaluated at any set of points {x (j) } Nx j=1 ⊆ Ω ∪ ∂Ω with the high fidelity solution u(m (i) ) evaluated at the same points. The prediction û(m (i) ; θ)(x (j) ) is straightforward to compute using Equation (4). This involves passing the reduced branch inputs m (i) and the coordinates of point x (j) into the branch net and trunk net, respectively. The label u(m (i) )(x (j) ) is obtained using finite element method solvers.• The second term is for learning the directional derivative of the evaluation u(x (j) ) with respect to the input function m, in the direction of the reduced basis ψ in 1 , . . . , ψ in r . It can be shown that∇ m û(m (i) ; θ)(x (j) )Ψ in = ∇ m û(m (i) ; θ)(x (j) ).Thus, the derivative of the outputs of the parametric model can be computed as the partial derivatives of the output with respect to the input of the branch net via automatic differentiation. On the other hand, the derivative labelsΦ out (x (j) )(∇ m u(m (i) )Ψ in ) = (du(m (i) ; ψ in 1 )(x (j) ), . . . , du(m (i) ; ψ in r )(x (j))) are obtained by first computing the Gateaux derivatives du(m (i) ; ψ in 1 ), . . . , du(m (i) ; ψ in r ) and then evaluating them at x (j) . See Appendix B.3 for details about the computation of the Gâteaux derivative du(m; ψ).• We initialize the loss weights λ 1 = λ 2 = 1 and choose a loss balancing algorithm called the self-adaptive learning rate annealing algorithm [23] to update them at a certain frequency. This ensures that the gradient of each loss term have similar magnitudes, thereby enabling the neural network to learn all these labels simultaneously. Remark. For the training, the computational cost of the second term of the loss function (and its gradients) largely depends on the number of points used in each iteration. To reduce the computational and especially memory cost, we can use a subset of points, N batch x = αN x , where α is small number between 0 and 1 (e.g., α = 0.1 in our experiments), in a batch of functions, though their locations could vary among different batches in one epoch. We find that this approach has little impact on the prediction accuracy of the model when N x is large enough.Throughout the paper, we assume that the parameter functions m are independently drawn from some Gaussian random field [24] . In particular, we consider the case where the covariance function is the Whittle-Matérn covariance function, that is, m i.i.d.∼ N ( m, C), where m is the (deterministic) mean function and C = (δI -γ∆) -2 (I is the identity and ∆ is the Laplacian) is an operator such that the square root of its inverse, C -1 2 , maps random function (m -m) to Gaussian white noise with unit variance [25] . The parameters δ, γ ∈ R + jointly control the marginal variance and correlation length.We consider two linear projection bases for dimension reduction of the parameter function.The KLE basis is optimal in the sense that the mean-square error resulting from a finite representation of the random field m is minimized [26] . It consists of eigenfunctions determined by the covariance function of the random field. Specifically, an eigenfunction ψ of the covariance operator C = (δI -γ∆) -2 satisfies the differential equationEQUATIONWhen solved using the finite element method, Equation ( 6) is equivalent to the following linear system (See Appendix A.2 for the derivation)EQUATIONwhere the (i, j)-entries of matices A in and M in are given byA in ij = δ⟨ϕ in j , ϕ in i ⟩ + γ⟨∇ϕ in j , ∇ϕ in i ⟩, M in ij = ⟨ϕ in j , ϕ in i ⟩.Here, we recall that V in h is the finite element function space of input function, ϕ in i , i = 1, . . . , N in h for the finite element basis of V in h , and ⟨•, •⟩ for the L 2 (Ω)-inner product. We select the first r (typically r ≪ N in h ) eigenfunctions ψ in 1 , . . . , ψ in r corresponding to the r largest eigenvalues for the dimension reduction. LetΨ in = (ψ in 1 | • • • |ψ in r )denote the corresponding nodal values of these eigenfunctions. Since the eigenvectors ψ in i are M in -orthogonal (or equivalently, eigenfunctions ψ in i are L 2 (Ω)-orthogonal), the reduced representation of m can be computed as m = (Ψ in ) T M in m, that is, the coefficients of the low rank approximation of m in the linear subspace spanned by the eigenfunctions ψ in 1 , . . . , ψ in r .Active Subspace Method (ASM) basis. The active subspace method is a gradient-based dimension reduction method that looks for directions in the input space contributing most significantly to the output variability [27] . In contrast to the KLE basis, the ASM basis is more computationally expensive. However, since the ASM basis captures sensitivity information in the input-output map rather than solely the variability of the input space, it typically achieves higher accuracy in predicting the output than KLE basis. We consider the case where the output is a multidimensional vector [22] , representing the nodal values of the output function u. The ASM basis ψ i , i = 1, . . . , r, are the eigenfunctions corresponding to the r largest eigenvalues of the generalized eigenvalue problemEQUATIONwhere the action of operator H on function ψ is given byEQUATIONHere, du(m; ψ) is the Gâteaux derivative of u at m ∈ V in h in the direction of ψ ∈ V in h , defined as lim ε→0 (u(m + εψ) -u(m))/ε, and d * u(m; •) is the adjoint of the operator du(m; •). When solved using finite element method, Equation ( 8) is equivalent to the following linear system (See Appendix A.3 for the derivation)EQUATIONwhere the action of matrix H on vector ψ is given byEQUATIONand the action of matrix C -1 on vector ψ is given byEQUATIONHere, M out denotes the mass matrix of the output function space, i.e., M out ij = ⟨ϕ out j , ϕ out i ⟩. In practice, when solving Equation (10) , we obtain its left hand side through computing(⟨Hψ, ϕ in 1 ⟩, . . . , ⟨Hψ, ϕ in N in h ⟩) Tand its right hand side through the matrix-vector multiplication in Equation ( 12) (See Appendix A.3 for details). Similar to the KLE case, let Ψ in denote the nodal values of r dominant eigenfunctions. Since the eigenvectors ψ in i are C -1 -orthogonal, the reduced representation of m can be computed as m = (Ψ in ) T C -1 m. We use a scalable double pass randomized algorithm [28] implemented in hIPPYlib to solve the generalized eigenproblems Equation (7) and Equation (10) .To this end, we present the computation of the derivative label and the action of H as follows.Theorem 1. Suppose the PDE in the general form of (1) is well-posed with a unique solution map from the input function m ∈ V in to the output function u ∈ V out with dual (V out ) ′ . Suppose the PDE operator R :V in × V out → (V out ) ′ is differentiable with derivatives ∂ m R : V in → (V out ) ′ and ∂ u R : V out → (V out ) ′ , and in addition ∂ u R is invertible with invertible adjoint (∂ u R) * .Then the directional derivative p = du(m; ψ) for any function ψ ∈ V in , and an auxillary variable q such that d * u(m; p) = -(∂ m R) * q can be obtained as the solution of the linearized PDEsEQUATIONEQUATION)Proof sketch. We perturb m with εψ for any small ε > 0 and obtain R(m + εψ, u(m + εψ)) = 0.Using Taylor expansion to expand it to the first order and letting ε approach 0, we obtain Equation ( 13), where p = du(m; ψ). Next we compute d * u(m; p). By Equation ( 13), we havedu(m; ψ) = -(∂ u R) -1 (∂ m R)ψ. Thus, d * u(m; ψ) = -(∂ m R) * (∂ u R) - * p.Then we can first solve for q = (∂ u R) - * p and then compute -(∂ m R) * q. See Appendix A.1 for a full proof. 
experiments: 
 In this section, we present experiment results to demonstrate the derivative-enhanced accuracy and cost of our method on two test problems of nonlinear vector-valued PDEs in comparison with DeepONet, FNO, and DINO. Details about the data generation and training can be found in Appendix B.In all test cases, we assume that the input functions m (i) are i.i.d. samples drawn from a Gaussian random field with mean function m and covariance operator (δI -γ∆) -α . We take α = 2 in two space dimensions so the covariance operator is of trace class [29] . It is worth noting that the parameters γ and δ jointly control the marginal variance and correlation length, for which we take values to have large variation of the input samples that lead to large variation of the output PDE solutions. In such cases, especially when the mapping from the input to output is highly nonlinear as in our test examples, a vanilla neural network approximations tend to result in relatively large errors, particularly with a limited number of training data. To generate samples from this Gaussian random field, we use a scalable (in mesh resolution) sampling algorithm implemented in hIPPYlib [28] .We consider two nonlinear PDE examples, including a nonlinear vector-valued hyperelasticity equation with one state (displacement), and a nonlinear vector-valued Navier-Stokes equations with multiple states (velocity and pressure). For the hyperelasticity equation, we consider an experimental scenario where a square of hyperelasticity material is secured along its left edge while a fixed upwardright force is applied to its right edge [30] . Our goal is to learn the map between (the logarithm of) the material's Young's modulus and its displacement. We also consider the Navier-Stokes equations that describes viscous, incompressible creeping flow. In particular, we consider the lid-driven cavity case where a square cavity consisting of three rigid walls with no-slip conditions and a lid moving with a tangential unit velocity. We consider that the uncertainty comes from the viscosity term and aim to predict the velocity field. See Appendix B.1 for details.We use the following three metrics to evaluate the performance of different methods. For the approximations of the solution u(m), we compute the relative error in the L 2 (Ω) norm and H 1 (Ω) norm on the test datasetD test = {(m (i) , u (i) )} Ntest i=1 , that is, 1 N test m (i) ∈Dtest ∥û(m (i) ; θ) -u(m (i) )∥ L 2 (Ω) ∥u(m (i) )∥ L 2 (Ω) , 1 N test m (i) ∈Dtest ∥û(m (i) ; θ) -u(m (i) )∥ H 1 (Ω) ∥u(m (i) )∥ H 1 (Ω) ,where∥u∥ L 2 (Ω) = Ω ∥u(x)∥ 2 2 dx = u T M out u, ∥u∥ H 1 (Ω) = (∥u∥ 2 L 2 (Ω) + ∥∇u∥ 2 L 2 (Ω) ) 1/2 .For the approximation of the Jacobian du(m; •), we compute the relative error (for the discrete Jacobian) in the Frobenius norm onD test along random directions ω = {ω i } Ndir i=1 , that is, 1 |D test | m (i) ∈Dtest ∥dû(m (i) ; ω) -du(m (i) ; ω)∥ F ∥du(m (i) ; ω)∥ F ,where ω i are samples drawn from the same distribution as m.We compare the prediction errors measured in the above three evaluation metrics and computational cost in data generation and neural network training for four neural operator architectures, including DeepONet [1] , FNO [2] , DINO [20] , and our DE-DeepONet. We also add experiments to demonstrate the performance of the FNO trained with the derivative loss (DE-FNO) and the DeepONet trained with input dimension reduction but without the derivative loss (DeepONet (ASM) or DeepONet (KLE)). For FNO, we use additional position embedding [4] that improves its approximation accuracy in our test cases. For DINO, we use ASM basis (v.s. KLE basis) as the input reduced basis and POD basis (v.s. output ASM basis [20] ) as the output reduced basis, which gives the best approximation accuracy. For the input reduced basis of DE-DeepONet, we also test and present the results for both KLE basis and ASM basis. Test errors. In Figure 1 , we show a comparison of different methods for the hyperelasticity and Navier-Stokes problem in predicting the solution (in L 2 norm, left, and H 1 norm, middle) and its derivative (in Frobenius norm, right) with respect to parameter in the direction of N dir = 128 test samples {ω i } Ndir i=1 . First, we can observe significant improvement of the approximation accuracy of DE-DeepONet compared to the vanilla DeepONet for all three metrics in all cases of training samples. Second, we can see that FNO leads to larger approximation errors than DINO and DE-DeepONet for all metrics except when number of training samples is 1024. We believe that the reason why FNO performs better than DE-DeepONet and DINO when training samples are large enough is mainly due to the use of input dimensionality reduction in DE-DeepONet and DINO (where the linear reduction error cannot be eliminated by increasing training samples) whereas in FNO we use full inputs. We also see that DE-FNO performs the best among all models when the training samples are sufficient (256 or 1024), although in the compensation of much longer training time shown in Figure 2 . Third, we can see that the approximation accuracy of DINO is similar to DE-DeepONet (ASM) but requires much longer inference time as shown in Table 5 (see Section 5 and Appendix B.3.2 for reasons). In DINO, the output reduced basis dimension is set to be smaller than or equal to the number of training samples as the output POD basis are computed from these samples, i.e., 16 for 16 samples and 64 for ≥ 64 samples. Increasing the output dimension beyond 64 does not lead to smaller errors in our test. Finally, we can observe that DE-DeepONet using ASM basis leads to smaller errors than using KLE basis, especially for the Navier-Stokes problem. Moreover, we present the output reconstruction error due to the input dimension reduction using KLE or ASM basis in Table 1 .The errors provide the lower bound of the relative errors in L 2 (Ω) norm of DINO and DE-DeepONet. We can see that using the ASM basis results in a lower output reconstruction error than the KLE basis (more significant difference observed in the more nonlinear Navier-Stokes equations). See Appendix B.6 for the decay of the reconstruction error with increasing number of reduced basis.In addition, we provide full visualizations of the ground truth and prediction of both solution and derivatives in Appendix B.7. Visually, DE-DeepONet (ASM) consistently provides the best estimate (in terms of pattern and smoothness similarity with the ground truth) for both the solution and its derivative when the training data is limited (16 or 64 samples).Data generation computational cost. We use MPI and the finite element library FEniCS [31] to distribute the computational load of offline data generation to 64 processes for the PDE models considered in this work. See Table 2 for the wall clock time in generating the samples of Gaussian random fields (GRFs), solving the PDEs, computing the r = 16 reduced basis functions (KLE or ASM) corresponding to the 16 largest eigenvalues, generating the derivative labels, respectively. In Table 3 , We also provide the total wall clock time of data generation of DE-DeepONet (ASM) (we only includes the major parts -computing high fidelity solution, ASM basis and dm labels [16 directions]) when N train = 16, 64, 256, 1024 using 16 CPU processors. 1 The inference time of DINO is dominated by the time required to compute evaluations of all finite element basis functions at the grid points using FEniCS (which may not be the most efficient, see Appendix B.3.2).Even though these grid points overlap with parts of the finite element nodes-allowing us to skip evaluations by extracting the relevant nodes-for a fairer comparison with DE-DeepONet (in terms of its ability to evaluate at any arbitrary point), we assume they are arbitrary points requiring explicit evaluation.Model training computational cost. We present comparisons of the wall clock time of each optimization iteration (with batch size 8) of different methods in Table 4 and convergence plot (error versus training time) in Figure 2 and the figures in Appendix B.5. We find that incorporating derivative loss leads to longer training time as expected. However, when the training data are limited, the increased computation cost is compensated for a significant reduction of errors. We note that there are potential ways to further reduce the training cost, e.g., by training the model with additional derivative loss only during the later stage of training, or by using fewer points for computing the derivative losses in each iteration. Additionally, thanks to the dimension reduction of the input, we can define a relatively small neural network and thus are able to efficiently compute the derivatives using automatic differentiation. 
related work: 
 Our work is related to Sobolev training for neural networks [32] , which was found to be effective in their application to model distillation/compression and meta-optimization. In the domain of surrogate models for parametric partial differential equations, our work is more closely related to derivativeinformed neural operator (DINO) [20] which is based on a derivative-informed projected neural network (DIPNet) [17] , and presents an extension to enhance the performance of the DeepONet. Compared to DINO, although the DeepONet architecture (and its formulation of dm loss) requires longer training time, it offers the following advantages: (1) Potentially shorter inference time. The additional trunk net (which receives spatial coordinates) allows us to quickly query the sensitivity of output function at any point when input function is perturbed in any direction. While DINO can only provide the derivative of the output coefficients respect to the input coefficients (we call reduced dm), in order to compute the sensitivity at a batch of points, we need to post process the reduced dm by querying the finite element basis on these points and computing large matrix multiplications; (2) Greater flexibility and potential for improvements. Although both DeepONet and DINO approximate solution by a linear combination of a small set of functions, these functions together in DeepONet is essentially the trunk net, which is "optimized" via model training, whereas in DINO, they are POD or derivative-informed basis precomputed on training samples. When using DINO, if we encounter a case where the training samples not enough to accurately compute the output basis, the large approximation error between the linear subspace and solution manifold will greatly restrict the model prediction accuracy (see Figure 2 when N train = 16). And the reduced dm labels only supports linear reduction of output. However, it is possible that we can further improve DeepONet by, e.g., adding physical losses (to enhance generalization performance) and Fourier feature embeddings (to learn high-frequency components more effectively) on the trunk net [5] and replacing the inner product of the outputs of two networks by more flexible operations [33, 9] (to enhance expressive power).The dm loss formulation of our work is broadly suitable any network architecture that has multiple subnetworks, where at least one of them receives high-dimensional inputs. 
discussion: 
 In this work, we proposed a new neural operator-Derivative-enhanced Deep Operator Network (DE-DeepONet) to address the limited accuracy of DeepONet in both function and derivative approximations. Specifically, DE-DeepONet employs a derivative-informed reduced representation of input function and incorporates additional loss into the loss function for the supervised learning of the derivative of the output with respect to the inputs of the branch net. Our experiments for nonlinear PDE problems with high variations in both input and output functions demonstrate that adding this loss term to the loss function greatly enhances the accuracy of both function and derivative approximations, especially when the training data are limited. We also demonstrate that the use of derivative loss can be extended to enhance other neural operators, such as the Fourier neural operator.We presented matrix-free computation of the derivative label and the derivative-informed dimension reduction for a general form of PDE problems by using randomized algorithms and linearized PDE solves. Thanks to this scalable approach, the computational cost in generating the derivative label data is shown to be only marginally higher than generating the input-output function pairs for the test problems, especially for the more complex Navier-Stokes equations which require more iterations in the nonlinear solves than the hyperelasticity equation.We require the derivative information in the training and dimension reduction using ASM, which may not be available if the explicit form of the PDE is unknown or if the simulation only provides input-output pairs from some legacy code. Another limitation is that dimension reduction of the input function plays a key role in scalable data generation and training, which may not be feasible or accurate for intrinsically very high-dimensional problems such as high frequency wave equations. Such problems are also very challenging and remain unsolved by other methods to our knowledge.A.1 Proof of Theorem 1We first show how to compute p given the PDE residual R(m, u) = 0. Since u is uniquely determined by m, we can write u = u(m) so that R(m, u(m)) ≡ 0 holds for any m ∈ V in . Thus, for any ε > 0 and ψ ∈ V in , we have R(m + εψ, u(m + εψ)) = 0. Using the Taylor expansion we obtainEQUATIONwhere δu = u(m + εψ) -u(m). Dividing both sides of Equation ( 15) by ε and letting ε approach 0 yields(∂ m R(m, u(m)))ψ + (∂ u R(m, u(m)))du(m; ψ) = 0. For ease of notation, we writeEQUATIONWe solve Equation ( 16) via its weak formEQUATION) where v is a test function in V out .Next we show how to compute d * u(m; p). By Equation ( 16), we have du(mEQUATIONLet q := (∂ u R) - * p. We can solve for q via the weak form ⟨(∂ u R) * q, v⟩ = ⟨p, v⟩, or equivalently, ⟨q,EQUATIONwhere v is a test function in V out .EQUATIONNote that we do not need to solve for w explicitly; we only compute ⟨w, v⟩ with v as the finite element basis functions ϕ in 1 , . . . , ϕ inN in h. The cost of computing the right hand side of Equation ( 20) arises from evaluating the directional derivative and its inner product with the finite element basis functions. Remark. By Equation (9), we use N grad samples {(m (i) , u (i) )} Ngrad i=1 to compute the Monte Carlo estimate of the action of operator H on any function ψ ∈ V in , that is,EQUATIONRemark. When using the double pass randomized algorithm to obtain the first r eigenpairs in Equation (10) , we need to compute the action of H on 2(r + s) random functions in V in (e.g., their nodal values are sampled from the standard Gaussian distribution), where s ∈ N + is typically a small oversampling parameter, often chosen between 5 and 20. To speed up the computation, we first compute the LU factorization of the matrices resulting from the discretization of the linear PDEs in Equation (17) and Equation (19) . Then the action of d * u(m (i) ; du(m (i) ; •)) on these random functions can be efficiently computed via the forward and backward substitution. Furthermore, the computational time can be significantly reduced by parallelizing the computation of the average value in Equation (21) across multiple processors.A.2 Proof of the equivalence of Equation ( 6) and Equation ( 7)By the definition of C = (δI -γ∆) -2 , Equation ( 6) is equivalent toEQUATIONWe first compute (δI -γ∆)ψ. To do this, let p = (δI -γ∆)ψ and multiply both sides of this equation by a test function v and integrateEQUATIONwhere, in the second equality, we use integration by parts and the assumption that the test function vanishes on the boundary. Then we substitute v with all of the finite element basis functions ϕ in 1 , . . . , ϕ in N in h and collect the corresponding linear equations for the nodal values of pEQUATIONwhere the (i, j)-entries of A in and M in are given by 22), function p satisfies ψ = λ(δI -γ∆)p. In the same manner we can see thatA in ij = δ⟨ϕ in j , ϕ in i ⟩ + γ⟨∇ϕ in j , ∇ϕ in i ⟩, M in ij = ⟨ϕ in j , ϕ in i ⟩. By Equation (EQUATION)Note that both matrices M in and A in are symmetric positive definite and thus nonsingular. Combining Equation (24) with Equation ( 25) yieldsEQUATIONor equivalently,EQUATIONA.3 Proof of the equivalence of Equation ( 8) and Equation (10) By Equation ( 8), for any test function v ∈ V in h , it holds that ⟨Hψ, v⟩ = ⟨λC -1 ψ, v⟩.In particular, for any m ∼ ν(m), as we let v go through all of the finite element basis functionsEQUATION)Indeed, by the definition of Gateaux derivative, we haveEQUATIONwhere ϕ out (x) = (ϕ out 1 (x), . . . , ϕ N out h (x)) are the finite element basis functions of output function space V out h . Then for any test function v ∈ V in h , it holds that ⟨d * u(m; p), v⟩ = ⟨p, du(m; v)⟩ (by the definition of adjoint operator)= ⟨ϕ out (x)(∇ m u)ψ, ϕ out (x)(∇ m u)v⟩ (by Equation (30)) = v T (∇ m u) T M out (∇ m u)ψ,where M out is the mass matrix of output function space V out h , i.e., M out ij = ⟨ϕ out j , ϕ out i ⟩ for 1 ≤ i, j ≤ N out h . Note that if we replace v by the i-th finite element basis functions ϕ in i , then v becomes the standard unit vector e i ∈ R N in h (with the k-th entry one and all others zero). Thus,⟨d * u(m; p), ϕ in i ⟩ = e T i (∇ m u) T M out (∇ m u)ψ, 1 ≤ i ≤ N in h .   ⟨d * u(m; p), ϕ in 1 ⟩ . . . ⟨d * u(m; p), ϕ in N in h ⟩    = (∇ m u) T M out (∇ m u)ψ.Next, we prove thatEQUATIONIndeed, if we let w = λC -1 ψ, then similar to the argument in Appendix A.2, we haveλA in (M in ) -1 A in ψ = M in w, Note that ⟨w, ϕ in i ⟩ = e T i M in w, 1 ≤ i ≤ N in h . Thus,    ⟨w, ϕ in 1 ⟩ . . . ⟨w, ϕ in N in h ⟩    = IM in w = λA in (M in ) -1 A in ψ.Combining ( 28), ( 29) and ( 31) yieldsE m∼ν(m) [(∇ m u) T M out (∇ m u)ψ] = λA in (M in ) -1 A in ψ.Hyperelasticity equation. We follow the problem setup in [30] . Write X (instead of x) for a material point in the domain Ω and u = u(X) : R 2 → R 2 for the displacement of the material point. Under the influence of internal and/or external forces, the material point is mapped to a spatial point x = x(X) = X + u(X) : R 2 → R 2 . Let F = ∇ X x = I + ∇ X u : R 2 → R 2×2 denote the deformation gradient. For a hyperelasticity material, the internal forces can be derived from a strain energy densityW (X, C) = µ(X) 2 (tr(C) -3) + λ(X) 2 (ln(J)) 2 -µ(X) ln(J).Here, C = F T F is the right Cauchy-Green strain tensor, tr(C) is the trace of matrix C, J is the determinant of matrix F , and µ(X), λ(X) : R 2 → R are the Lamé parameters which we assume to be related to Young's modulus of elasticity E(X) : R 2 → R and Poisson ratio ν ∈ Rµ(X) = E(X) 2(1 + ν) , λ(X) = νE(X) (1 + ν)(1 -2ν).We assume the randomness comes from the Young's modulus E(X) = e m(X) + 1. Let S = 2 ∂W ∂C denote the second Piola-Kirchhoff stress tensor. We consider the case where the left boundary of the material is fixed, and the right boundary is subjected to stretching by an external force t = t(X) : R 2 → R 2 . The strong form of the steady state PDE can be written as       ∇ X • (F S) = 0, X ∈ Ω, u = 0, X ∈ Γ left , F S • n = 0, X ∈ Γ top ∪ Γ bottom , F S • n = t, X ∈ Γ right ,where Γ left , Γ right , Γ top and Γ bottom denote the left, right, top, and bottom boundary of the material domain Ω, respectively, and n is the unit outward normal vector on the boundary. Our goal is to learn the operator that maps the parameter m to the displacement u. For demonstration, we choose m = 0, δ = 0.4, γ = 0.04, Poisson ratio ν = 0.4, and the external forcet(X) = 0.06 exp(-0.25|X 2 -0.5| 2 ), 0.03(1 + 0.1X 2 ) T .In practice, we solve the PDE by first formulating the energy W in the weak form Navier-Stokes equations. Let u = u(x) ∈ R 2 and p = p(x) ∈ R denote the velocity and pressure at point x ∈ Ω = (0, 1) 2 . The strong form of the PDE can be written asW = Ω W dX -       -∇ • e m ∇u + (u • ∇)u + ∇p = 0, x ∈ Ω, ∇ • u = 0, x ∈ Ω, u = (1, 0) T , x ∈ Γ top , u = (0, 0) T , x ∈ Γ \ Γ top ,where Γ top and Γ denote the left and whole boundary of the cavity domain Ω, respectively. Here, we assume that the randomness arises from the viscosity term e m . Our goal is learn the operator that maps the parameter m to the velocity u. For demonstration, we choose m = 6.908 (e m ≈ 10 3 , thus the viscosity term dominates), δ = 0.6, and γ = 0.06. See Figure 4 for the visualization for one parameter-solution pair of the Navier-Stokes equations. For all PDEs in this work, we use the class dolfin.UnitSquareMesh to create a triangular mesh of the 2D unit square with 64 cells in horizontal direction and 64 cells in vertical direction. For the Darcy flow equation and hyperelasticity equation, we set the direction of the diagonals as 'right', while for the Navier-Stokes equation, we set the direction of the diagonals as 'crossed'. See Figure 5 for a visualization of the unit square mesh with a 10 by 10 resolution. We use the class dolfin.FunctionSpace or dolfin.VectorFunctionSpace to create the finite element function space of input function and output function. For the finite element basis functions, we consider the Continuous Galerkin (CG) family (or the standard Lagrange family) with degree 1 or 2. See Table 6 for details. We generate N train = 1500 and N test = 500 input-output pairs (m (i) , u (i) ) for training and testing, respectively. We compute first r = 16 KLE basis and ASM basis using double pass randomized algorithm with an oversampling parameter s of 10. In the computation of ASM basis, we use N grad = 16 samples for the Monte Carlo estimate of the action of operator H in Equation ( 21). In our method, we formulate the different labels into arrays with the shape as follows• Evaluation labels: (N train , N x , N u ) • Derivative m labels: (N train , N x , r, N u )Recall that N train is the number of functions used for training, N x is the number of nodes of the mesh, N u is the dimension of output function, r is the number of reduced basis, and d is the dimension of the domain. Since the PDE residual R(m, u(m)) ≡ 0 holds for any m ∈ V in h , we have R(m + εψ, u(m + εψ)) = 0, ∀ε > 0, ψ ∈ V in h . By the Taylor expansion and R(m, u(m)) = 0, we obtainEQUATIONwhere δu = u(m + εψ) -u(m). Dividing both sides of Eq. ( 32) by ε and letting ε approach 0 yieldsEQUATIONwhere p := du(m; ψ). We solve Eq. ( 33) for p via its weak form⟨(∂ m R)ψ, v⟩ + ⟨(∂ u R)p, v⟩ = 0,where v is a test function in V out h .Example. Consider the (nonlinear) diffusion-reaction equation-∇ • (e m ∇u) + u 3 = 1, x ∈ Ω, u(x) = 0, x ∈ ∂Ω. Then • R = -∇ • (e m ∇u) + u 3 -1 • (∂ m R)ψ = -∇ • (e m ψ∇u) • (∂ u R)p = -∇ • (e m ∇p) + 3u 2 pThus, p satisfies the linear PDE ⟨e m ψ∇u, ∇v⟩ + ⟨e m ∇p, ∇v⟩ + ⟨3u 2 p, v⟩ = 0.Using FEniCS, we can easily derive Gâteaux derivative via automatic symbolic differentiation instead of by hand. In this case, the Python code for representing the weak form of the residual ⟨R, v⟩ and Gâteaux derivatives ⟨(∂ m R)ψ, v⟩ and ⟨(∂ u R)p, v⟩ can be written as import dolfin as dl û(m; θ)(x) = ⟨b((Ψ in ) T C -1 m; θ b ), t(x; θ t )⟩ + θ bias ,whereΨ in = (ψ in 1 | • • • |ψ in rin ) ∈ R N in h ×rinare the nodal values of input (ASM) reduced basis functions, C -1 is the inverse of the covariance matrix of the Gaussian random field ν where parameter m is sampled from (Recall that the ASM basis are orthonormal in the inner product with weight matrix C -1 .) Thus, by the chain rule of derivative, for any test direction ψ test , one has∇ m û(m; θ)(x)ψ test = ∇ m ⟨b( m), t(x; θ t )⟩(Ψ in ) T C -1 ψ test ,where m = (Ψ in ) T C -1 m. The Jacobian matrix ∇ m ⟨b( m), t(x; θ t )⟩ can be efficiently computed using, e.g., torch.func.jacrev, and further parallelized with torch.vmap. If ψ test is in Ψ in during model training, we can see that (Ψ in ) T C -1 ψ test becomes a unit vector which frees us the need to compute it, otherwise in the model inference stage when ψ test is (the nodal values of) a random function sampled from ν, we compute T = (Ψ in ) T C -1 and then T ψ test .û(m; θ)(x) = Φ out (x)Ψ out f θ ((Ψ in ) T C -1 m),where Φ out (x) = (ϕ out 1 (x), . . . , ϕ outN out h (x)) ∈ R 1×N out h denotes the output finite element basis functions evaluated on point x ∈ Ω, Ψ out = (ψ out 1 | • • • |ψ out rout ) ∈ R N out h ×routare the nodal values of output (POD) reduced basis functions, similarly for Ψ in denoting the input (ASM) reduced basis, and C -1 is the inverse of the covariance matrix of the Gaussian random field where parameter m is sampled from. Thus, by the chain rule of derivative, for any test direction ψ test , one has∇ m û(m; θ)(x)ψ test = Φ out (x)Ψ out ∇ m f θ ( m)(Ψ in ) T C -1 ψ test ,where m = (Ψ in ) T C -1 m. For fast evaluations given different x, m, and ψ test , we first compute T left = Φ out (x)Ψ out for all points x that need to be evaluated and T right = (Ψ in ) T C -1 . Next, we compute J = ∇ m f θ ( m) using, e.g., torch.func.jacrev, and finally compute T left JT right ψ test .DeepONet & FNO. Both models receive the full high dimensional parameter m, so we compute the directional derivative ∇ m û(m; θ)ψ test using the Jacobian-vector product torch.jvp instead of the full Jacobian in DE-DeepONet and DINO. For the coordinates x as additional inputs, we pad zeros to the tensor ψ test to match the dimension of the input tensor (m, {x (j) } Nx j=1 ).For the DeepONet, we parameterize the branch net using a CNN and the trunk net using a ResNet.For the FNO, we use the package neuraloperator. For the DINO, we use 16 ASM basis functions for the input dimension reduction and 64 POD basis functions for the output dimension reduction, and parameterize the neural network using a ResNet. For the DE-DeepONet, both the branch net and trunk net are parameterized using ResNets. We train each model for 32768 iterations (with the same batch size 8) using an AdamW optimizer [34] and a StepLR learning rate scheduler (We disable learning rate scheduler for DE-DeepONet). See Tables 7 to 10 for details. When the loss function comprises two terms, we apply the self-adaptive learning rate annealing algorithm from [23] , with an update frequency of 100 and a moving average parameter of 0.9, to automatically adjust the loss weights {λ i } 2 i=1 in Equation (5) . Additionally, we standardize the inputs and labels before training. To measure the error induced by the projection, we define the output reconstruction error as follows1 N N i=1 ∥u(P r m (i) ) -u(m (i) )∥ L 2 ∥u(m (i) )∥ L 2 ,where P r is the rank r linear projector. We provide the plots of the output reconstruction error vs number of reduced basis r using KLE basis and ASM basis in Figure 12 . We can see that using the ASM basis results in a lower output reconstruction error than the KLE basis. • It should be clear whether the error bar is the standard deviation or 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 a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments? Answer: [Yes] Justification: We provide the machine information and time of execution in the Tables 4 and 6 .Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines? Answer: [Yes] Justification: We have reviewed the code of ethics and checked that the research follows them.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 10. Broader Impacts Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed? Answer: [NA] Justification: Our work lies in the field of scientific machine learning. We do not notice any societal impact for now. Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.• The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: We think our research poses no such risk.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes] Justification: We cite the corresponding papers and respect the licenses when using the packages FEniCS, hIPPYlib, and neuraloperator to conduct experiments.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.38th Conference on Neural Information Processing Systems (NeurIPS 2024). 
 
 Appendix:  


Decision: The paper introduces a modified version of DeepONet, termed Derivative-Enhanced DeepONet, by incorporating derivative terms relative to the spatial domain into the loss function. The paper use KL expnasion/active subspace methods to approximate the derivative so that all the computation can be implemented in the function space. This paper implements solid experiments and showed the perfomance boost of the proposed loss. The newly introduced trunknet alows faster inference net. This is definitely a solid paper. Thus I vote for acceptness.

Reviews: review: paper_summary: paper_summary: The paper introduces a modified version of DeepONet, termed Derivative-Enhanced DeepONet, by incorporating derivative terms relative to the input function and spatial domain into the loss function. It outlines a practical method for calculating these derivatives, which serve as supplementary supervision terms in the training process. The authors demonstrate enhanced performance over other Neural Operator baselines with this approach, particularly in low-data scenarios, using two datasets centered on hyperelasticity and Navier-Stokes equations.
strength_weakness: strength_weakness: - As a reviewer with limited expertise in functional analysis, I found the sections detailing the supplementary supervision terms challenging to comprehend. A more introductory overview explaining how the derivative ground truths are derived would be beneficial for clarity and accessibility.
questions: questions: - Could you provide a higher-level explanation of how the derivative labels are calculated? Are they derived from the underlying equation?
- On line 325, you mention reducing training costs by introducing additional derivative losses in the later stages of training. Do you have any data or figures to support this claim?
- Why does supervision of the derivatives become less significant as more data becomes available?
limitations: limitations: Limiatations have been addressed
paper_summary: paper_summary: The paper proposes an extension of deep operator networks (DeepONets) enhanced by matching the derivatives of the output function with respect to the input function, for example, the derivatives of the PDE solution w.r.t. the input coefficient. To make the computation tractable, the dimensionality of the input function is reduced via a dimension reduction technique, active subspace method, avoiding expensive Jacobian evaluations.
strength_weakness: strength_weakness: - The results of the  empirical evaluation do not seem to be strong enough to make the proposed method look like an effective alternative to the existing method (FNO). The gain obtained at the cost of increased computation seems to be marginal. 

- Although KLE or ASM is a standard method for reducing dimensions of a field, the number of required bases depends on smoothness or regularity of the field and, thus, the proposed method could be benefitted in some specific scenarios. Some discussions on where this method could be benefitted and where this method would struggle are required. 

- Presenting wall time per epoch is informative, but it would provide a more complete picture if the authors could provide the wall time for training to achieve a certain level of accuracy. 

- Regarding the statement “when the training data are limited, the increased computation cost is compensated for a significant reduction of errors”: It would be informative if the authors could provide a summarizing figure presenting a result for varying number of training instances and show that the benefits of the proposed method is more pronounced in the data-scarce regime. 

- Although it is not in the main body of the paper, Figure 2 seems to miss the entire information. This omission undermines the completeness of the paper and gives the impression that the paper was prepared hastily.
questions: questions: - Eq (4). Does the branch net take the reduced input vector as an input? Could the branch net just take the original input and approximate the Frechet derivative using the reduced representation? 

- Can the authors describe how the Frechet derivative is computed (the derivative with respect to the reduced input) to collect the ground truth data? 

- Is the reason for using the relative error for measuring the loss for preventing numerical issues caused by different numerical scales?
limitations: limitations: - A thorough empirical investigation on in which scenarios this proposed method outperforms baseline (and also the opposite scenarios) would be needed.
paper_summary: paper_summary: This manuscript introduces a derivative-enhanced deep operator network that utilizes derivative information to improve prediction accuracy.
strength_weakness: strength_weakness: W1) The baselines appear limited as recent benchmarks like GNOT, LSM, ONO, Transsolver, etc., have ye
t to be compared. Numerous variations of FNO have been introduced, such as FFNO, UNO, CoNO, GFNO, UFNO, etc.

W2) The benchmarks dataset appears to be restricted to datasets for Navier-Stokes equations and hyperelasticity, which are not open-source.

W3) ASM is costly in computation, and derivative information is necessary when the PDE form is unknown. Additionally, it needs to scale better with dimensions.

W4) The proposed method builds upon the DINO paper by introducing additional informed losses.

W5) The new losses were found effective only for DeepOnet and lacked experimental evidence compared to FNO.
questions: questions: Q1) The primary distinction between DINO and the proposed approach lies in utilizing DeepOnet and spatial derivatives. However, spatial derivatives were also introduced in [1]. Also, have you used the approximation for spatial derivative as used in [1]?

[1] DEEP MULTI-SCALE VIDEO PREDICTION BEYOND MEAN SQUARE ERROR. Michael Mathieu, Camille Couprie & Yann LeCun.

Q2) What is the impact of proposed losses when combined with other operators like FNO? Do they enhance the performance of different operators?

Q3) Could approximating the Frobenius norm along a random direction be analogous to score-based diffusion models?

Q4) Is there a specific reason why adding new loss terms relative to L2 still results in poorer performance than FNO when dealing with hyperelasticity as training samples increase?

Q5) How does incorporating this noisy gradient version aid in training?
limitations: limitations: Yes, the authors adequately addressed the limitations.
paper_summary: paper_summary: The authors propose incorporating derivative information into the training of DeepONets to improve predictive accuracy on PDE problems. Instead using an encoder-decoder neural operator architecture as in prior work, the authors use the DeepONet architecture as well as incorporate spatial derivative information on top of functional derivatives previously used.  The authors provide numerical comparisons with other models using both types of derivative training as well as results for two types of input space dimensionality reduction methods.
strength_weakness: strength_weakness: -	The cost of generating the dm and dx labels with which to incorporate into training was 2-3 times greater than the cost of solving the PDE itself given Table 1. Of course, as with all neural operators the question is the tradeoff between offline data generation and training versus the benefit of quick online predictions so this cost could be mitigated given the accuracy benefit it provides. Nevertheless, its quite the increase and one may wonder if its worth doing at that point compared to just using the high fidelity FEM solver for N number of problems. Unless N is high enough, the tradeoff looks less appealing here. 

-	I would like to see the authors rephrase their findings in the context of the results in B.4 which show that the dx information provides minimal to no benefit compared to the dm training, and it some cases it even makes the accuracy worse than dm alone. This comparison of dx and dm alone should be stated in the main text as it appears to be critical information. Given that beside the solution accuracy, the dm accuracy is vastly more important that the dx accuracy to perform PDE-constrained optimization problems, it begs the question why dx regularization is done at all here. This then begs the question how novel this is compared to DINO which already trained on dm information. The authors claim Sobolev learning is novel on top of the DINO results but it does not appear impactful. The authors need to address this. 

-	Following onto the prior point, I am on the edge regarding impact and novelty here. That being said I still think the work is informative and of high enough quality to publish, albeit without as strong of an impact due to the main benefits shown to be in line with the existing method DINO. Looking at DE-DeepONet vs. DINO in Figure 1, there is not much difference, and in Table 2 DINO is shown to be much faster.  I would consider changing my rating if this was adequately addressed.
questions: questions: -	Could the authors please include total wall-clock time to Figure 2 in addition to per epoch. This would be helpful in comparing to the data cost in Figure 1 which is a total. Additionally, the authors may want to consider including a row for end-to-end data and training time such that the DeepONet and FNO include the PDE solution generation time and DINO and DE-DeepONet include the data for all three loss terms. This way the reader can very easily compare the total cost associated with making predictions with a DeepONet versus DE-DeepONet. I think this would benefit the manuscript since it mitigates the cost of the dm dx label generation. 

-	How can the epochs only be 1,000 in B.3? In the original DeepONet paper 50,000 – 500,000 iterations were used. Are all models converged to fairly compare them? What do the convergence plots look like? It would be helpful to see them. 

-	Is the vanilla DeepONet and FNO also trained with the dimensionality reduction technique on the input, and if so, which one? Figure 1 shows KLE and ASM for DE-DeepONet, and DINO is stated to use ASM, but what about the baseline models. How is it a fair comparison if the inputs are not the same?

-	The usage of CNN for the trunk and ResNet for the branch should not be hidden in the appendix. The vanilla DeepONet uses a fully-connected NN for both the trunk and branch and that would be assumed here, please mention it in 3.1. Additionally, what justification is there for those choices? The CNN for the trunk makes sense to construct the basis but I’m curious about the ResNet for m?

-	The authors should consider the following papers [2,3]. In [3] the authors train DeepONets the physics-informed (dataless) way using gradient information to obey the governing equation. In [2], the authors train a PINN with gradient-enhanced information which could also be done with the previously mentioned PI-DeepONet in [3]. How might DE-DeepONet methodology be incorporated into these models and in what ways is it distinct from them? It would be nice to see this discussion in the main or appendix text. 

[2] Yu, Jeremy, et al. "Gradient-enhanced physics-informed neural networks for forward and inverse PDE problems." Computer Methods in Applied Mechanics and Engineering 393 (2022): 114823.

[3] Wang, Sifan, et al. "Learning the solution operator of parametric partial differential equations with physics-informed DeepONets." Science advances 7.40 (2021): eabi8605.

-	How were the gradients computed using automatic-differentiation (AD)? I do not see the ML package mentioned, a package like Jax has more accurate and substantially faster AD than one like PyTorch, see [4] Table 6 and Figure 15. The authors may want to consider this to improve performance and reduce the overhead cost for dm and dx label generation. 

[4] Jagtap, Ameya D., et al. "How important are activation functions in regression and classification? A survey, performance comparison, and future directions." Journal of Machine Learning for Modeling and Computing 4.1 (2023).
limitations: limitations: -	Adequately described in manuscript.

Comments: comment: We thank the reviewer for the updated response!
I thamk the authors for their clarifications and additional experiments, which help me understand the paper better and answer my questions. I have raised my score from 5 to 6.
Thank you for your updated response. We're glad that the new results help address your concern to some extent. When revising the manuscript to include the new results, we will also keep your (and other reviewers') feedback in mind.
Thank you for providing the clarification and the experimental results. The new information on training/inference time is informative and clears my concern to some extent. I adjusted the score accordingly.
We greatly appreciate your positive feedback and support. Considering how dm loss performs with noisy datasets is indeed an interesting and meaningful direction to explore, but we believe it would be more natural to study this within the context of, say, inverse problems, which is beyond the scope of our manuscript (our work only considers enhancing surrogate models for forward operators governed by known PDEs). 

In terms of readlability, we notice some parts of the manuscript might be too dense or abstract for readers that are not from heavy math backgrounds. When revising the manuscript to includer the new results, we will also pay attention to adding more intuitive explanations (although sometimes at the expense of accuracy), concrete examples, or computational details to help a broader audience understand, particularly in how to derive both DE-DeepONet's (same for FNO+dm) and DINO's dm labels in practice.
Thank you for the authors' addressing of my comments, weaknesses, and questions. New experimental results were conducted to help the readers understand the proposed losses more appropriately. Concerning the last question, I wanted to see how proposed losses perform when we have noisy datasets. Does it improve the performance?  

Although the paper should be revised with the above information before publication, the author continues to polish the manuscript. Lastly, in the revised version, try to make it a little bit more readable for someone not from a purely mathematical background. So, I am willing to raise the score. Thank you for the responses. I have raised my score and hope to see the new changes incorporated into the revised version.
We sincerely appreciate your feedback and support. We will definitely incorporate the new results into the revised manuscript.
I thank the authors for their response to my review and particularly the new results provided in the overarching rebuttal. In the context of these items, a number of my points have been addressed (removing dx from training, gaining improvement over DINO with the new FNO + dm results, including convergence plots, etc.) and therefore I am raising my score from 6 $\rightarrow$ 7 under the condition these new (necessary) results make it into the revised manuscript before publication and don't die here.
Dear reviewers,

could you have a look at the authors response and comment on them if you have done so, yet.

thanks in advance

your area chair