Derivative-enhanced Deep Operator Network
Authors: Yuan Qiu
Nolan Bridges
Peng Chen
Venue: NeurIPS 2024 poster (2024)
Keywords: Neural operators
Operator learning
Derivative learning
DeepONet
Dimensionality reduction
Adjoint method
Primary area: deep_learning_architectures
TL;DR: In this work we propose a derivative-enhanced deep operator network (DE-DeepONet), which leverages the derivative information to enhance the prediction accuracy, and provide a more accurate approximation of the derivatives.

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.

Derivative-enhanced Deep Operator Network
Yuan Qiu, Nolan Bridges, Peng Chen
Georgia Institute of Technology, Atlanta, GA 30332 {yuan. qiu, bridges, pchen402/@gatech..
Abstract
The deep operator networks (DeepONet), a class neural operators that learn mappings between function spaces, have recently been developed as surrogate mod- els for parametric partial differential equations (PDEs). In this work we propose derivative-enhane deep operator network (DF-DeepONet), which leverages derivative information to enhance the solution prediction accuracy and provides a more accurate approximation of solution-to toraammeteeerr 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-solutee 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 our approach.
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 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, 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 the input, which are often needed for many downstream tasks such PDE- constrained optimization problems for control, inference, and experimental design [11, 12, 13, 14).
In this paper, propose enhance the performance of DeepONet through derivative- based dimen- sion 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-enhance proximation 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 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 reproduce the results this paper can be found https://github.. http://hhhhh.cc//8800001
approaches can significantly improve DeepONet's ss ap-
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

2 Preliminaries
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.
2.1
Problem setup
Let Ω C Rd denote an open and bounded domain with boundary an C R where the dimension 3. We consider PDE of the general form defined in f2 as
(m, u) 0,
(I)
prescribed with proper boundary conditions. Here m E Vin is an input parameter function defined in a separable Banach space V with probability measure and u / is the output as the solution of the PDE defined in separable Banach space Vout, Our goal is to construct parametric model i(;; 0) to approximate the solution operator that maps the parameter m to the solution u.
Once constructed, the parametric model ü m: 9) should be much more computationally efficient to evaluate compared solving the PDE with high fidelity approximation.
2.2 High fidelity approximation
For the high fidelity approximation the solution, we consider using a finite element method [21] in this work. We partition the domain Ω into finite set of subregions, called cells or elements. Collectively, these cells form mesh of the domain Ω2 Let represent the diameter of the largest cell. We denote Vi indexed by h as the finite element space for the approximation of the input space Vin with Lagrange basis {dif,.. Ne of dimension Ni. such that (f(()) = dig the finite element node 20)) with d2y being the Kronecker delta function. Similarly, we denote Vi as the finite element space for the approximation of the solution space yourr with basis { G
of dimension Na. Note that for the approximation to be high fidelity, Nj, and No are often very large. To this end, can write the high fidelity approximation of the input and output functions as Nu
N
my (x)
Σ
mi (x) and un (x)
Σ
u (
(x),
with coefficient vectors (m,," m Ny )T N and u (w1, (N ) R.V.". whose entries are the nodal values of and u at the corresponding nodes.
2.3 DeepONet
We briefly review the DeepONet architecture [1] with focus on learning the solution operator of the PDE in Equation (1). To predict the evaluation of solution function at any point E Ω U an for any given input function m, [1] design a network architecture that comprises two separate neural networks: trunk net E(- 0.). which takes the coordinate values of the point at which we want to evaluate the solution function as inputs, and a branch net b(- ; ;d)), which takes the vector m encoding the parameter function m inputs. In [1], the vector is the function evaluations at
finite set of fixed points {2(0) [== C Ω U On, that is, m (m) mmx((()))) ())))) 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 obtained by taking the standard inner product between these two vectors and (optionally) adding real-valued bias parameter, i.e.
a(m; (x) = (b(m;; 0p), (x;; 0z)) + Obias;
with (0b, Or, Opins) the solution u is vector-valued (u(1)... ...... (N.)), as in our experiments, we can use one of the four approaches in [3] con- struct the DeepONet. Specifically, for each solution component, we use the same outputs of branch net with dimension N₂ and different corresponding groups of outputs of trunk net to compute the inner product. More precisely, the solution ( of component 1,..., N.. approximated by
N₂, components, i.e.,
(m; (z) (b(m; 0s), t (z; 0,)) bias
(2)

---

where (1; 01) = t(x; 8,) 1).6 + 1: :N:, the vector slice of t(2; 01) with indices ranging from (1-1))N to iN, ...,,,,.,,, Nu..
For this construction, the outputs of the branch net can he interpreted as the coefficients of the basis learned through the trunk net. By partitioning the outputs of the trunk net into different groups, essentially partition the basis functions used for predicting different components of the solution. The DeepONet is trained using dataset D = {[m!), u J
No with ND samples, where m are random samples independently drawn from and (()) are the solution of the PDE (with slight abuse of notation from the vector-valued solution) at m(*), 1,...,ND
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 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-enhancee DeE-eepONee), 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 0 denote the trainable parameters the DE-DeepONet. We propose to minimize the loss function
DeepONet
L(0)- A1 Bmm-w//luee (l((m)) -(m)) + (|du(m;))) ddūmm;
(3)
where the (22) norm of square integrable function I is defined as II/lI(=
(a (x=)l dz)1/², the Hilbert-Schmid norm of an operator T H + H that acts on Hilbert space H is defined as || = Eier ||Te,|ly, where {eiher is an orthonormal basis of H. Here, \1, > 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 [[d((m;)) dm(; -) liås (approximately) becomes V Vmue V V lI. where the Frobenius norm of a matrix M E Rm is defined as M = > 2
m M2, It is a critical challenge to both
-i=1
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 size N x Njn, with Ni. Nou >> 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 low dimensional linear space spanned by basis w/in ,/.,.... 4/-i TR. Ni, with < Nia that is,
((15"" m....,, (1b"" m))" ER".
To better leverage the information of both the input probability measure v and the map u: m - u(m). consider the basis generated by active subspace method (ASM) using derivatives [22]. ASM iden- tifies 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 term Em [du[m,;))
dū(m;;)) Il can be approximated by N,No Σ Σjn || V (x)()) - Vmū(m (8) (x()) (
(0)
with small amount of functions ( sampledd from input probability measure v and points 201) in the domain Ω. Note that V, (2)(()))) is vector of size which is computationally feasible. For comparison, we also conduct experiments if the basis is the most commonly used Lorvenen-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.
3.1 Model architecture
We incorporate dimension reduction into DeepONet to construct a parametric model for approximat- ing the solution operator. Specifically, if the solution is scalar-valued (the vectur- case can be constructed similar to (2)), the prediction is given by
ä(m; (z) (b(m; 0,), f(x;; 0,)) + O
(4)
3

---

where 0 = (fb. B2, (uas) are the parameters to be learned. The branch net b(:;6) and trunk net 0,) can be chosen as an MLP, ResNet, etc. Note that the branch net takes small vector of the projected parameter as input. We also apply the Fourier feature embeddings [5], defined as ((x) = [cos{ Bæ), sin Rx)], to the trunk net, where each entry in B E Rmxdd sampled from Gaussian distribution N00, 2 and m E N+, IR+ are hyper-parameters.
3.2 Loss function
In practical training of the DE-DeepONer, formulate the loss function as follows
No
((0) =
enr cm(((á(m ())))) ((()))))))))mm)))))()))))))
ND
(5)
A2
e ((((((((m))) ()) (z())))) yin a (x (V um (1)))))))) )
N.
where (-(6) N are the nodes of the mesh - (wi") | the matrix collecting the nodal values of the reduced basis of the input function space, and क.nh (a:) = (x)) dout Npmi (z)) the
vector-valued function consisting of the finite element basis functions of the output function space.
The err( (((()) b())) denotes the relative error between any two groups of vectors (8) and ((),
computed
9 13)1/
err( (a (() 1i=1
2i=1 B 1/2
where 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 i(m (8); () evaluated at any set of points { (x())) S Z aQ with the high fidelity solution (f) evaluated at the same points. The prediction ú (0) (z(()) is straightforward to compute using Equa- tion (4). This involves passing the reduced branch inputs m and the coordinates of point z05) into the branch net and trunk net, respectively. The label u(m()) (2)(0) is obtained using finite element method solvers.
The second term is for learning the directional derivative of the evaluation u(x())) with respect to the input function m, in the direction of the reduced basis wh ...., with It can be shown that
Umu(m (1) (z) yp = Vn (m (1) (z)
Thus, the derivative of the outputs of the parametric model can be computed as the par- tial derivatives of the output with respect to the input of the branch net via automatic differentiation. On the other hand, the derivative labels
gou (z) (V m u(m (+) "
(du(m
v0th (2(0)))
are obtained by first computing the Gateaux derivatives du(m w/"" ...., dul (mm())))) with) and then evaluating them at (J), See Appendix B.3 for details about the computation of the Gâteaux derivative du(m; v).
We initialize the loss weights A₁- A2 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) largety depends on the number of points used in each iteration. To reduce the computational and especially memory cost, we can use a subsel of points, Nbatch aN,, where a is small number beiween O and 1 (e.g., = 0.1 in our experiments), in batch of functions, though their locations could vary among different batches one epoch. We find that this approach has little impact on the prediction accuracy of the model when N₂ is large enough.

---

3.3 Dimension reduction
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 t.ddd (Mm, C),, where m. is the (deterministic) mean function and C (8I ) (I is the identity and A is the Laplacian) is operator such that the square root of its inverse, C- à maps random function (m m) to Gaussian white noise with unit variance [25]. The parameters 6, y E R+ jointly control the marginal variance and correlation length.
We consider two linear projection bases for dimension reduction of the parameter function.
Karhumen-LLovv Expansion (KLE) basis. The KLE basis is optimal the sense that the mean-square error resulting from finite representation of the random field is minimized [26]. It consists of cigenfunctions determined by the covariance function of the random field. Specifically, an cigenfunction d of the covariance operator C = (6l-664)) -> satisfies the differential equation
Cy Xxp.
When solved using the finite element method, Equation (6) is equivalent to the following linear system (See Appendix A.2 for the derivation)
M (A - = AM
where the (enties of matices A and M are given by
1 (0)) in ( + VV in Vol" MI (1)), (0
Here, we recall that Vith is the finite element function space of input function, dft 1,..., Ni, for the finite element basis of Vi and ( for the -in product. We select the first r (typically < _iit)) eigenfunctions qian ...... dimension reduction. Let Trin (vii- ... (bb")) denote the corresponding nodal values of these eigenfunctions. Since the eigenvectors v are - "M"-orthogonal (or equivalently, eigenfunctions with are 1.2 (2)) that is, the coefficients of the low rank approximation of m in the linear subspace spanned by the eigenfunctions wi " viin
....,." 4/th corresponding the largest eigenvalues for the
the reduced representation of m can be computed as m = (T Mimm,
Active Subspace Method (ASM) basis. The active subspace method is 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, typically achieves higher accuracy in predicting the output than KLE basis. We consider the case where the output a multidimensional vector [22]. representing the nodal values of the output function The ASM basis vh, 1,..., r, are the eigenfunctions corresponding to the r largest eigenvalues of the generalized eigenvalue problem
Hb XC- v..
(8)
where the action of operator H. on function vb is given by
Hor (m)e(m)) (d u(m; du(m; x))]
Here, du(m; 1) is the Gâteaux derivative of u at m E V in the direction of a E Vin defined as lime-o(u(m + =œ) u/m))/E, and d' (((m; -) is the adjoint of the operator du(m;:). When solved using finite element method, Equation (8) equivalent to the following linear system
(See Appendix A.3 for the derivation)
Hap -C
(10)
where the action of matrix H on vector 2b is given by
v/ E Em-u(m) um) M (V u)
(11)
and the action of matrix C-1 on vector 4b is given by
~ (M
(12)
5

---

Here, M denotes the mass matrix of the output function space, i.e., M = (d)" , ) In practice,
when solving Equation (10), we obtain its left hand side through computing
(HHb, of....., (Ha,, o Ny
and its right hand side through the matrix- multiplication in Equation (12) (See Appendix A.3 for details). Similar to the KLE case, let in denote the nodal values of dominant eigenfunctions. Since the eigenvectors 2h!h are C- -othhooggonall, the reduced representation of m can be computed as m (1'")"'C We use a scalable double pass randomized algorithm [28] implemented in hIPPYlib 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 E V to the output function
Vou with dual (Vmily. Suppose the PDE operator R. : V x Vone - (Vruay' is differentiable with derivatives 0n, R. v - -Vorry and O.R: : Vou + (Vony, and in addition O. is invertible with invertible adjoint (0.R)). Then the directional derivative p = du(m; 1b) for any function ub E V and an auxillary variable such that d' u((;;)) = -0 (0m R)² can be obtained the solution of the linearized PDES
(0. R)p R)mp = 0,
(13)
(0. q = p.
(14)
Proof sketch. We perturb m with eú for any small e > 0 and obtain R(m + stb u(m + Et)) = 0.
Using Taylor expansion to expand it to the first order and letting e approach 0 we obtain Equation (13),
where = du(m; v5). Next we compute d' ummmppp) p). By Equation (13), we have du(m (m; 1/) =
R)) 11 (3m R) Thus, d' umm;; +) R)) ((,,R)
Then we can first solve for q
(0,R) p and then compute- - R)
R) q. See Appendix A.1 for full proof.
4 Experiments
In this section, we present experiment results to demonstrate the derivative-enhancee 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.
4.1 Input probability measure
In all test cases, we assume that the input functions (+) are i.i.d. samples drawn from a Gaussian random field with mean function vi and covariance operator (87 7 We take a = 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, vanilla neural network approximations tend 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].
4.2 Governing equations
We consider two nonlinear PDE examples, including nonlinear vector-valued hyperelasticity equation with one state (displacement), and 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 upward- right force is applied its right edge [30]. Our goal is to learn the map between (the logarithm of) the material's 's Young's 's modulus and its displacement. We also consider the Navier-Stokes equations that describes viscous, incompressible creeping flow. particular, we consider the lid-driven cavity case where square cavity consisting of three rigid walls with no-slip conditions and lid moving with tangential unit velocity. We consider that the uncertainty comes from the viscosity term and aim to predict the velocity field. See Appendix B.I for details.

---

4.3 Evaluation metrics
We use the following three metries evaluate the performance different methods. For the approximations of the solution u(m), we compute the relative error in the (2) norm and H' (S2) norm on the test dataset Diest = ((m(+) u () that is,
[ 0) - [lu(m(())
[ ;0) ((
Σ
N ED
New mm
[|u(m(()) (1
where
1 [m(x))] dx = u² Mo [ (llull220)) + [V
||
For the approximation of the Jacobian -), we compute the relative error (for the discrete Jacobian) in the Frobenius norm Diest along random directions = (w that is,
Idū(((()::) du(ml(;w)lle [[du(m())w)]]l
E
Desa
where w; are samples drawn from the same distribution as m.
4.4 Main results
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 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- DE-DeepONet, we also test and present the results for both KLEE basis and ASM basis.
100%
Extuciv avori Cito san bertiir 227681
200%
heldiim ena
32766)
200%
Ta ha kmviN bexiirn 22390
28%
120%
25%
309
30%
15%
20%
205
1044
Felttrr everts fenterra
Eate ever
279
37%
2000
39%
2004
27h-
33%
100%
220
5
305
75-
75c
58
Figure 1: Mean relative emror (# standard deviation) 5 random seeds of neural network training for a varying number of training samples for the [top: hyperelasticity: bottom: Navier-Stokes) equation using different methods. Relative errors in the L²( norm (left) and (I2 norm (middle) for the prediction of u (Mi, U2). Right: Relative error in the Frobenius (Fro) norm for the prediction of du(m;; w) (du, (du;(m:)) duz (u))m;)))
7

---

cmyy nor Pet
500%
200
100%
55
10%
30%
Flaiing Sin tarords
Flarin tine Jecores
tritie (H0E
Aahalee ever leet
50%
57h
270
15%
15
Figure 2:
Mean relative errors (#: standard deviation) over 5 random seeds versus model training time for the Navier-Stokes equations when the number of training samples is [top: 16; bottom: 256].
Test errors. In Figure 1, show comparison of different methods for the hyperelasticity and -tokersStokes problem in predicting the solution (in L norm, left, and H norm, middle) and its derivative (in Frobenius norm, right) with respect to parameter in the direction of Nur 128 test samples )2 Na, First, we 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 DB-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 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 Third, we can see that the approximation accuracy of DINO is similar to DE-DeepONet (ASM) but requires much longer inference time as shown Table (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.
Table 1: Output reconstruction error with 16 input reduced bases
Moreover, we present the output reconstruction error due to the imput dimension reduction using KLE or ASM basis in Table 1. The eтгоιs provide the lower bound of the relative errors in ( norm of DINO and DE-DeepONet. We can see that using the ASM basis results in lower output reconstruction error than the KLE basis (more significant difference observed in the more nonlinear equations). See Appendix B.6 for the decay of the reconstruction error with increasing number of reduced basis.
Relative L erroe
Dataset
KLE
ASM
Hypenelasticity
3.8%
2.7%
Navicr-Stokos
17.4%
5,8%
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).
8

---

Data generation computational cost. We use MPI and the finite element library FEniCS [31] to distribute the computational load of offlime 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 Nirin = 16, 64, 256, 1024 using 16 CPU processors.
Table 2: Wall clock time (in seconds) for data generation on 2 x AMD EPYC 7543 32-Core Processors
Dataset
GRFs (Na 2000)
PDEs (Na 2000)
KLE (r 16)
ASM (r
dm labels 2000, (64 procs)
(Na
(641 proes)
(64 proes)
(1 proes)
(16 proes)
Hyperelasticity
1.1
9.7
0.4
1.4
19.5
Navicr-Stokes
99.1
9.7
125.5
Table 3: Wall clock time (in seconds) for Table Wall clock time (seconds/iteration with batch size 8) for training single NVIDIA RTX A6000 GPU
data generation with different number of training samples using 16 CPU processors
PDEs + ASM husis dm lahels
Dataset
Dataset', Nean
16
64
256
1024
DeepONet
PNO
DINO
DE-INO
DE DeepONet
Hyperelasticity
2
5
16
61
Hyperelasticity Navier-Stokes
0.007
0.015
0.007
0.215
0.022
Navier-Stokes
17
38
124
470
0.007
0.015
0.007
0.216
Table 5: Total wall clock time (in seconds) for each model inferring on 500 test samples of both the solution and dm in 128 random directions, using single GPU and single CPU (except where specified)
Model
Dataset
DeepONet
FNO
DINO GPU CPU/16 CPUs
DE DE-DeepONet
Numerical solver 0 GPU + 16 CPUs
Hyperelasticity
166
Navier-Stokes
33
2152/151
18
1103
The inference time of DINO is dominated by the time required compute evaluations of all finite element basis functions at the grid points using FEniCS (which may not be the most efficient, Appendix B.3.2). Even though these grid points overlap with parts of the finite element nodes aaloowal extracting the relevant nodes-fit for a fairer comparison with DE-Deep DeepOOeee (in terms its ability to evaluate at any arbitrary point), we assume they are arbitrary points requiring explicit evaluation.
skip evaluations by
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 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 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 relatively small neural network and thus are able to efficiently compute the derivatives using automatic differentiation.
5
Related work
Our work related to Sobolev training for neural networks [32], which was found to be effective in their application to model distillation/comprrssion and meta-optimization. In the domain of surrogate

---

models for parametric partial differential equations, our work is more closely related to derivative- informed neural operator (DINO) [20] which based on a irtvatiie-informedd 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 the output coefficients respect to the input coefficients (we call reduced dm), in order to compute the sensitivity at batch of points, 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-informedd basis precomputed on training samples. When using DINO, if we encounter 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 Niriin 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-free freh-frey 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 new neural seratto--Derivatvveenaaaaeee Deep Operator Network (DE-DeepONet) to address the limited accuracy of DeepONet in both function and derivative approximations. Specifically, DE-DeepONet employs 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 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 be extended to enhance other neural operators, such as the Fourier neural operator.
We presented matrix-free free computation of the derivative label and the derivative-informee dimension reduction for general form of PDE problems by using randomized algorithms and lincarized PDE solves. Thanks to this scalable approach, the computational cost in generating the derivative label data is shown be only marginally higher than generating the input-output function pairs for the test problems, especially for the more complex Savier-Stoker equations which require more iterations in the nonlinear solves than the hyperelasticity equation.
Limitations: 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 if the simulation only provides input-output pairs from some legacy code. Another limitation is that dimension reduction of the input function plays 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 also very challenging and remain unsolved by other methods our knowledge.
10

---

Acknowledgements
We would like to thank Jinwoo Go, Dingcheng Luo and Lianghao Cao for insightful discussions and helpful feedback.
References
[1] L. Lu, P. Jin, Pang, Z. Zhang, and E. Karniadakis, "Learning nonlinear operators via Deep- onet based on the universal approximation theorem of operators," Nature machine intelligence, vol. 3, no. 3, pp. 218-229, 2021.
[2] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and Anand- kumar, "Fourier neural operator for parametric partial differential equations." arXiv preprint arXiv:2010.0889. 2020.
[3] I.. Lu, X. Meng. S. Cai, Z. Mao, S. Goswami, Z. Zhang, and G. E. Karniadakis, "A comprehen- sive and fair comparison of two neural operators (with practical extensions) based on fair data," Computer Methods in Applied Mechanics and Engineering. vol. 393, p. 114778, 2022.
[4] N. Kovachki, Z. Li, B. Liu, K. Azizzadenesheli, K. Bhattacharya, A. Stuart, and A. Anandkumar, "Neural operator: Learning maps between function spaces with applications to PDEs," Journal of Machine Learning Research, vol. 24, no. 89, pp. 1-97, 2023.
[5] S. Wang, H. Wang, and P. Perdikaris, "I.earning the solution operator of parametric partial differential equations with physics-informed deeponets," " Science advances, vol. 7, no. 40, eabi8605, 2021.
[6] S. Goswami, M. Yin, Y. Yu, and G. E. Karniadakis, "A physics-informed variational deeponet for predicting crack path in quasi-brittle materials," Computer Methods Applied Mechanics and Engineering, vol. 391, p. 1 115577, 2022.
[7] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, A. Stuart, K. Bhattacharya, and A. Anandkumar, "Multipole graph neural operator for parametric partial differential equations," Advances in Neural Information Processing Systems, vol. 33, pp. 6755-6766, 2020.
[8] B. Raonić, R. Molinaro, T. Rohner, S. Mishra, and E. de Bezenac, "Convolutional Neural 2023.
Operators," arXiv preprint arXiv:2302.01
[9] Z. Hao, Z. Wang, H. Su, C. Ying, Y. Dong, S. Liu, Z. Cheng, J. Song, and J. Zhu, "Gnot: general neural operator transformer for operator learning," in International Conference on Machine Learning, pp. 12556-12569, PMLR, 2023.
[10] A. Tran, A. Mathews, L. Xie, and C. S. Ong, "Factorized fourier neural operators," arXiv preprint arXiv:2111.13800 2021.
[11] D. Luo. T. O'Leary-Roseberry, Chen, and O. Ghattas, "Efficient PDE-constrained optimization under high-dimensional uncertainty using derivative-informed neural operators," arXiv preprint 2XXv:23355.2053. 2023.
[12] L. Cao, T. "Leary-Roseber and O. Ghattas, "Efficient geometric Markov chain Monte Carlo for nonlinear Bayesian inversion enabled by derivative-infin -intoivve---omeedd neural operators," arXiv preprint arXiv:2403.088 082200 2024.
[13] J. Go and P. Chen, "Accelerating Bayesian optimal experimental design with derivative-infommed neural operators," urXiv preprint arXiv:2312. 14810, 2023.
[14] Go and Chen, "Sequential derivative-informee latent attention neural operator," arXiv preprint arXiv:2409.091 2024.
Bayesian optimal experimental design with
[15] P. G. Constantine, E. Dow, and Q. Wang, "Active subspace methods in theory and prac- tice: applications to kriging surfaces," SIAM Journal on Scientific Computing, vol. 36, pp. A1500-A1524, 2014.
[16] Chen and O. Ghattas, "Projected Stein variational gradient descent," Advances in Neural Information Processing Systems, vol. 33, pp. 1947-1958, 2020.
[17] T. "'Leary-Rosebers U. Villa, P. Chen, and O. Ghattas, "Derivative-inormeed projected neural networks for dige-dimensional parametric maps governed by PDEs," " Computer Methods in Applied Mechanics and Engineering, vol. 388, p. 1 114199, 2022.

---

[18] O. Zahm, T. Cui, K. Law, A. Spantini, and Y. Marzouk, "Certified dimension reduction in nonlinear Bayesian inverse problems," Mathematics of Computation, vol. 91, no. 336, pp. 1789-
18 1835, 2022.
[19] H. Son, J. W. Jang, W. Han, and H. J. Hwang, "Sobolev training for the neural network solutions of PDES," 2020.
[20] T. "'Leary-Rsebenn" P. Chen, U. Villa, and O. Ghattas, "Derivative-informe neural opera- tor: an efficient framework for nigh-dimensional parametric derivative learning," " Journal of Computational Physics, vol. 496, p. 112555, 2024.
[21] J. T. Oden and J. N. Reddy, An introduction to the mathematical theory of finite elements. Courier Corporation, 2012.
[22] O. Zahm, P. G. Constantine, C. Prieur, and Y. M. Marzouk, "Gradient-based dimension reduction of multivariate vector-valued functions," SIAM Journal on Scientific Computing, vol. 42, no. 1, pp. A534-A558, 2020.
[23] S. Wang. S. Sankaran, H. Wang. and P. Perdikaris, "An expert's guide to training physics- informed neural networks," arXiv preprint 2rXiv22308.08468, 2023.
[24] R. J. Adler and E. Taylor, Random fields and geometry. Springer Science Business Media, 2009.
[25] F. Lindgren, H. Rue, and J. Lindström, "An explicit link between gaussian fields and Gaussian markov random fields: the stochastic partial differential equation approach," " Journal of the Royal Statistical Society Series B: Statistical Methodology, vol. 73, no. 4, pp. 423-498, 2011.
[26] R. G. Ghanem and P. D. Spanos, Stochastic finite elements: a spectral approach. Courier Corporation, 2003.
[27] P. G. Constantine, Active subspaces: Emerging ideas for dimension reduction in parameter studies. SIAM, 2015.
[28] U. Villa, N. Petra, and O. Ghattas, "HIPPYlib: an extensible software framework for large-scale inverse problems governed by PDEs: part I: deterministic inversion and linearized Bayesian inference," ACM Transactions on Mathematical Software (TOMS), vol. 47, no. 2, pp. 1-34, 2021.
[29] A. M. Stuart, "Inverse problems: a BByesiia perspective," " Acta numerica, vol. 19, pp. 451-559, 2010.
[30] L. Cao. T. O'Leary-Rosebens P. K. Jha, J. T. Oden, and O. Ghattas, "Residual-bare error correction for neural operator accelerated infinie-diimensinall bayesian inverse problems," Journal of Computational Physics, vol. 486, p. 1 112104, 2023.
[31] A. Logg, K.-. Mardal, and G. Wells, Automated solution of differential equations by the finite element method: The FEniCS book, vol. 84. Springer Science & Business Media, 2012.
[32] W. M. Czarnecki. S. Osindero, M. Jaderberg, G. Swirszez, and R. Pascanu, "Sobolev training for neural networks," Advances in neural information processing systems, vol. 30, 2017.
[33] S. Pan, L. Brunton, and J. N. Kutz, "Neural implicit flow: a mesh-agnostic dimensionality reduction paradigm of spatio-temporal data," Journal of Machine Learning Research, vol. 24, no. 41, pp. 1-60, 2023.
[34] Loshchilov and F. Hutter, "Decoupled weight decay regularization," arXiv preprint axXiv.171.005101, 2017.
12

---

A Proofs
A.1 Proof of Theorem
We first show how to compute given the PDE residual R(m, u) = 0. Since u is uniquely determined by m, write u(m) that R(m, (m)) = 0 holds for any E V Thus, for any 0 and to E Vin, we have
R(m + Ev, a(m + Ev))
Using the Taylor expansion obtain
(Om R(m, u(m))))= + (0.R(m, u(m)))ōu = 0,
(15)
where du u(m Ev) u(m). Dividing both sides of Equation (15) by and letting approach 0 yields
(0. R(m, (( (8,, (mm, u(m))))du(m; up)
For case of notation, write
8m = 6m R(m, u(m)),
O,R = 3 R,RRmm u(m)).
Then p is the solution the linear PDE
RM + (0,R)p 0.
(16)
We solve Equation (16) via its weak form
(O R)bb ((O,R)p, v) =0
(17)
where u is a test function in you
Next we show how compute d' d'u(m;p). By Equation (16), we have du(m; 4)) ((0.R) Thus,
w d'ummppp)) (- R))
(18)
= -0m () (0,R)) ~p.
Let = (0.RR "p. We can solve for q via the weak form
((0,R) qg u) = (p, v),
equivalently,
(g, (0,RRe) = (p, e),
(19)
where u a test function in Your
By Equation (18), we have w (0mm R)* 22 q. For any test function E Vin it holds that
(w,, v) { (0 v)
(20)
(0.. R.yo)
Note that we do not need to solve for w explicitly; we only compute (w, v) with v as the finite
element basis functions oir, ,..,, in. The cost 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 Nynd samples {mm(t) ( ) estimate of the action of operator H on any function v/ E Vin, that is,
Nows to compute the Monte Carlo
Ngnt
Hib
Σ d" ('u(m ; (0(mm);))) +)).
(21)
grad
Remark. When using the double pass randomized algorithm to obtain the first T eigenpairs in Equation (10), we need compute the action Hf on 2(r + 3) random functions in V (e.g., their nodal values are sampled from the standard Gaussian distribution), where E M+ is typically 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 dul (m((m))))))) :)) on these random functiins 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.
13

---

A.2 Proof of the equivalence of Equation (6) and Equation (7)
By the definition of C (óI - > -2 Equation (6) is equivalent to
t/ X(O7 >
(22)
We first compute (61 - ~ To do this, let p (ó1 >A) and multiply both sides of this equation by a test function v and integrate
1 p())u(x) dx S (OI - (()=))) de
8 1 (z))((z) dr +7 (Vb(x)) Ve(a)) dr,
(23)
where, in the second equality, we use integration by parts and the assumption that the function vanishes the boundary. Then we substitute with all of the finite element basis functions of ,.., dia and collect the corresponding linear equations for the nodal values of p Ng
A M
(24)
where the -entrrss of A and M are given by
A
(0)) in w + "( Vol"
M = (d)), p--)-
By Equation (22), function p satisfies X 7A)p. In the same manner we can see that
A Md
(25)
Note that both matrices M and A are symmetric positive definite and thus nonsingular. Combin- ing Equation (24) with Equation (25) yields
AA - M
(26)
or equivalently,
M (A Min
Min = XM
(27)
A.3 Proof of the equivalence of Equation (8) and Equation (10)
By Equation (8), for any test function E Vi it holds that
(Hyb, v) (AC-+ v).
(28)
In particular, for any m ~ v(m), as we let v go through all of the finite element basis functions 6 E Vin, we can show that
(d" u(m; du(m; vb)), oip)
= (V um)))" (V
(29)
(d'u(m; du (m: +),, N
Indeed, by the definition of Gateaux derivative, we have du(m;; 3) = lim
u(m + =w) - u(m)
802
N (m + EQb) (m)) @i
lim
do (x:)
N
(m, +
m No E um " m
(30)
lim
ooul
(x)
Nau
Tui 4h + ...
Ou: sb Ny )o (x:
Om
v (z)( V. u222
14

---

where (x)) = (61tt (2), Ne (z)) are the finite element basis functions of output function space V.
Then for any test function E V. it holds that
(d" u(mmp) v) (p, du(m;; v))
(by the definition of adjoint operator)
(I)) u) 60th (v) (by Equation (30))
= v" (Vmu)" M (Vm"m)
where M is the mass matrix of output function space VAA", i... M (d ) for 1 S i,j S Na Note that if we replace o by the i-th finite element basis functions o then a becomes the standard unit vector e E RRNE (with the k-th entry one and all others zero). Thus,
(d' u(m: p), ) o (Vmu) u) M (Vmu)u, Isis N Concatenating the above equations yields
(d' u(m; of)
(mmu) Mon (
(d (d'u(m;p) )
Next, we prove that
(2 o/:)
(M ub.
(31)
(XC- tln
Indeed, if we let XC- then similar to the argument in Appendix A.2, we have
(M u/ '
Note that
(w, (.) = T isis N
Thus,
(w,, (
TM -1 1m
(w, Np)))
Combining (28), (29) and (31) yields
Em (v 2 (V - 2b.
B Experimental details
B.1 Governing equations
Hyperelasticity equation. We follow the problem setup in [30]. Write X (instead of x) for material point in the domain Q and u = =(X) IR2 - 2 for the displacement of the material point. Under the influence of internal and/or external forces, the material point is mapped to spatial point x = xX)) = X + u(X): : I² + IR2 Let F Vxx = + Vxu : IR2 + _ denote the deformation gradient. For hyperelasticity material, the internal forces can be derived from strain energy density
((()) 2 (tr(C) - 3) +
A(X)
W(X C)
2 (In(J))² (-XXx((XX))) In(J).
15

---

Here, -= C=FFFF is the right Cauchy-Grev Green strain tensor, tr(C') is the trace of matrix C, J is the determinant of matrix F. and μ(X),, X(X) IR? + are the Lamé parameters which we assume to be related to Young's modulus of elasticity E(X): TR lR and Poisson ratio E R
B(X) ((X) = +1 X(X) = (1+ 2n)
VE(X)
We assume the randomness comes from the Young's modulus E(X) em((X) + Let oW
ac
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) R2 + RR? The strong form of the steady state PDE can be written as
Vx (FS) = O, u 0,
X E Ω
X E Tie:;
FS
X E T Top
bottom:
FS n =
E T right:
where left, I right top and I bottom denote the left, right, top, and bottom boundary the material domain Ω. respectively, and is the unit outward normal vector on the boundary. Our goal is to learn the operator that maps the parameter m the displacement u. For demonstration, we choose 5 = 0.4, 7 = 0.04, Poisson ratio v = 0.4, and the external force
(X)=
0.51²), 0.03(1 + 0.1X2))
In practice, solve the PDE by first formulating the energy W in the weak form
[t
W = wax
u) ds
Jaas
and then solving for u that satisfies the stationary condition, that is, the equation
dW (u; 0,
holds for any test function v in the state space. See Figure 3 for the visualization for one parameter- solution pair of the hyperelasticity equation.
Parameter
Spatial point x X u(X)
1.0
1.0
0.05
-004
0.6
X
0.03 M|
0.4
0.4
0
0.2
0.2
--2
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.6
0.8
1.0
Figure 3: Visualization of one parameter-solution pair of hyperelasticity equation. The color of the output indicates the magnitude of the displacement (which maps from domain St to R2) instead of its componentwise function 2j or ug. The skewed square shows locations of any domain point after deformation X + æ. See Figures 13 and 14 for uj and ug.
16

---

Navier-Stokes equations. Let = u(x) I2' and = p(z) E IR denote the velocity and pressure
point r f2 (0, 1) The strong form of the PDE can be written as
V Vut (u Vju Vp= 0,
!2
V u = 0,
S2,
u (1,0)T u (0,
x E I top:
E T Tippp
where top and I denote the left and whole boundary of the cavity domain Ω, respectively. Here, we assume that the randomness arises from the viscosity term e" Our goal is learn the operator the velocity For demonstration, choose
that maps the parameter thus the viscosity term dominates), 0.6, and 0.06. See Figure for the visualization for one parameter-solution pair of the
6.908 (e ss 10
equations.
Parameter
Velocity x-component
Velocity y-component
1.0
0.3
8
0.8
0.8
0.8
0.8
0.6
0.6
0.6
0.6
0.1
0.1
0.0
0.0
0.2
0.0
0.0
0.2
0.4
0.6
0.8
3.0
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.6
0.0
1.0
Figure 4: Visualization of one parameter-solution pair of Navier-Stokes equations.
B.2 Data generation
For all PDEs in this work, the class dolfin. UnitSquareMash create 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 the diagonals as 'right', while for the Navier-Stokes equation, we set the direction of the diagonals as 'crossed'. See Figure for visualization of the unit square mesh with 10 by 10 resolution.
10 by unit square mesh (right)
10 by unit square mesh (crossed)
1.0
1.0
0.8
0.8
0.4
0.4
0.2
0.2
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.6
0.8
1.0
Figure 5: Visualization of the 10 by 10 unit square mesh. Left: 'raggon Right: diago- cal='crossed'
17

---

We use the class dolfin. FunctionSpaan or dolfin. VectorFunctionSpee to create the finite element function space of imput 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 See Table 6 for details.
Table Configurations of data generation for different datasets
Daraset
Configantions
Mesh
(Nusia Nive)
Naaii
Hyperelantieity
64 x 64, right
CG (ID, deg-1)
CG (2D, deg-1)
(1500, 500)
16
16
10
Navier-Stolies
G4 64, crossed
CG (ID, deg-1)
OG (2D) deg-2), (ID
(1500, 500)
16
16
10
We generate Nuain 1500 and Nust 500 imputouutuu pairs (m(*, ())) for training and testing, respectively. We compute first = 16 KLE basis and ASM basis using double pass randomized algorithm with an oversampling parameter of In the computation of ASM basis, we use Ngra 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: (Nirain No. Nw) Derivative m labels: (Nuains Nz, No)
Recall that Muain is the number of functions used for training, N₂ is the number of nodes of the mesh, N₂ the dimension of output function, r is the number reduced basis, and d is the dimension of the domain.
B.3
Computation of derivative labels and outputs
B.3.1 Derivative labels p du(m; (b) as ground truth
Since the PDE residual R(m, (m)) = 0 holds for any E Via we have
R(m + Ex; u( + =v)) = 0, Ve > 0, 2 E v.
By the Taylor expansion and R,m, u(m)) = 0, we obtain
(0m R(m, =(m) + R0,, u(m)) ss 0, where u(m evp) - v(m). Dividing both sides of Eq. (32) by & and letting approach yields R) + R), O,
(32)
(33)
where p du(m;; v). We solve Eq. (33) for p via its weak form
(0,,. R) u) ((0, r) 0,
where is test function in Vouat
Example.
Consider the (nonlinear) diffusion-reaction equation
{
-V (e" + + 1, E Ω u(x) = 0, x E 02.
Then
R= (e" + -
Rm - VV)
R - (em + 3u²p
Thus, p satisfies the linear PDE
(e VV, Vv) + Vp, Vo) + (3u²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 ((mm)) v) and ((,, R),,, u) can be written as
import dolfin dl
18

---

R= (dl. inner exp.expppo (m) #dl. ( dl. (.)d((())) *dl. dx (u=+3-dd.Constant(( (1)1)))))))))))) *v*dl.dx)
(R..
dl. Merivatiie(R,u..)
B.3.2 Derivative outputs of neural networks
DE-DeepONet. For notation simplicity, we illustrate the case where the input reduced basis is ASM basis and the output function is real-valued. The output of the model is given by
ú(m; (z) = (T ( Bo), (xx; 0e)) + Obius:
where Viin (4/T) ... ( E Ni Xm= are the nodal values of input (ASM) reduced basis func- tions, C-1 is the inverse the covariance matrix of the Gaussian random field v where parameter m 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 shes one has
V m úi (m; (2 (z)t Vm (6( t(x; 8 (T" C
P/rest.
where i = (ID")T m. The Jacobian matrix V ()), (x;; 0z)) can be efficiently computed using, e.g., torch. func. jacrev, and further parallelized with torch. vmap. If thes is in T during model training, we can see that (Trnn C becomes a unit vector which frees us the need to compute it, otherwise in the model inference stage when e/re is (the nodal values of) random function sampled from v, we compute 7' = (yin)*C-1 and then I Teprest-
DINO.
The output of the model is given by
((m=0)(z) = hou vou So ((Tri))) C-)m,,
where (x)"(() = ( (x), (I)...., o Nan (x)) E R x No denotes the output finite element basis functions
evaluated on point Ω, Youut (z/1tl reduced basis functions, similarly for ii 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.
(pb R.N. are the nodal values of output (POD)
Thus, by the chain rule of derivative, for any test direction sausa one has
m (m; (xr (x"(x) ( V (i.
Whest:
where m (WinnTTC-m For fast evaluations given different æ, m. and thes we first compute Tien = hout (z) yrout for all points x that need to be evaluated and Trght =(D"")TC- Next, we compute info (m) using, e.g., torch. func. fta.ffunc...rcree.e.o jacrev, and finally compute JTighh thesst
DeepONet & FNO. Both models receive the full high dimensional parameter m, compute the directional derivative Vmu(m f)) Preat using the Jacobian vector product torch. jvp instead of the full Jacobian in DE-DeepONet and DINO. For the coordinates as additional inputs, we pad zeros the tensor the match the dimension of the input tensor (m, {x(1) )
B.4 Training details
For the DeepONet, we parameterize the branch net using 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 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 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 moving average parameter 0.9, to automatically adjust the loss weights { in Equation (5). Additionally, standardize the inputs and labels before training.
19

---

Table 7: Training details for DeepONet
Dataset
Hyperelasticity
Navier-Stokes
branch net
CNN
CNN
6 hidden layers
6 hidden layers
256 output dim
256 output dim
ReLU
ReLU
trunk net
ResNet
ResNet
3 hidden layers
3 hidden layers
512 output dim
512 output dim
512 hidden
512 hidden dim
ReI.U
ReL.U Kaiming Uniform
initialization
Kaiming Uniform
A (Ir, weight decay)
(10 3, 9,10-) (0)
(10 -3 9,10-9)
StepLR (gamma, step size)
(0.7, 10)
number of Fourier features
64
64
Fourier feature scale
0.5
1.0
Table 8: Training details for FNO & DE-FNO
Dataset
Hyperelasticity
Navier-Stokes
number of modes in channels out channels hidden channels numher of layers lifting channel ratio projection channel ratio activation function AdamW (Ir, weight decay) StepI.R (gamma, step size)
(32,32)
(32,32) 32)
64
64
4
4
GELU
GELU
(5 10-3 10-1) (0.9, 25)
(5 10-3, 10-1) (0.9,50)
Table 9: Training details for DINO
Dataset
Hyperelasticity
Navier-Stokes
neural network
ResNet
ResNet
3 hidden layers
hidden layers 64 output dim
64 output dim
128 hidden dim
256 hidden dim
EL.U
EL. Kaiming Normal (5 10-3 10-3,10-14 110
initialization
Kaiming Normal
AdamW (Ir. weight decay)
(5 x 10 3000)..)) (0.5, 50)
StepLR (gamma, step size)
(0.5, 50)
20

---

Table 10: Training details for DE-DeepONet
Datasct
Hyperelasticity
Navier-Stokes
branch net
ResNet
ResNet
3 hidden layers
hidden layers
128 output dim
256 output dim
hidden dim
256 hidden dim
ELU
ELU
trunk net
ResNet
ResNet
3 hidden layers
hidden layers
256 output dim
512 output dim
256 hidden dim
512 hidden dim
ReLUU
ReLU Kaiming Uniform
initialization
Kaiming Uniform
AdamW (Ir, weight decay)
(10-3, 110-$,10-1))
(10-3 3,-3,00 "100"1)
disable scheduler
True
True
number of Fourier features
64
64
Fourier feature scale o
0.5
1.0
Nbauch (= o.N2)
422 x (0.1 65
422 s (0.1 x 65")
Table 11: Number of trainable parameters in each model
parameters
Dataset
DeepONet
FNO & DE-FNO
DINO
DE-DeepONet
Hyperelasticity
4.21 M
17.88 M 17.88 M
0.04 M 0.15 M
0.27 M
Navier-Stokes
4.21 M
1.02 M
B.5
Convergence plot
Based on the training time (seconds/tteration) of each model in Table 4, we obtain the convergence
plots when the number of training samples is either limited (Nurii 16, 64) or sufficient (Nurio = 256, 1024) in Figures for the hyperelasticity equation, and in Figures 2, 10 and 11 for the Navier-Stokes equations.
Relid cam Pain 101
300%
112%
100%
55%
100%
20%
25"
405
355
20%
20%
385
31C
T6C Touriin Scoreds
35
10
Tuaring Scowed
Figure 6: Mean relative errors (+t standard deviation) over 5 trials versus model training time for the hyperelasticity equations when the number of training samples is 16.
21

---

55h
45%
100%
185
205
54%
42%
36%
185
30%
20%
20%
106
30
10
10
Figure 7: Mean relative emrors (# standard deviation) over 5 trials versus model training time for the hyperelasticity equations when the number of training samples is 64.
For ertee ener 2100 Phi 250
2045/
Adtavve ant in IEE same MU.ran
200%
PASAIN and FO 1or SIVE (NLIEI 250)
105
LO0%
100%
20%
28k
40%
45%
105
25%
04
194
10%
350
250
Flarinn
Figure 8: Mean relative errors (+ standard deviation) over trials versus model training time for the hyperelasticity equations when the number of training samples is 256.
Satais Litet LG241
108%
100%
50%
30%
109
30%
10%0
2%
32
oeceedsl
16
12
Traning tIme Decones
10
32
1
Traiing time becorest
WV
Figure 9: Mean relative errors (#: standard deviation) over 5 trials versus model training time for the hyperelasticity equations when the number of training samples is 1024.
Enteive evEr tor Inin
>alrimm any reri 25h
300%
200%
59
W
Tulees Seored
Thaeeiy See Eocreeg
127
Praniin Lime Seooeen
Figure 10: Mean relative errors (+: standard deviation) over 5 trials versus model training time for the Navier-Stoker equations when the number of training samples is 64.
22

---

1004)
6.0
59%
2004
27%
100%
19
04
10%
14
25
114
154
Flaiing Sin Records
Farii tin Jecree
Figure 11: Mean relative emrors (# standard deviation) over 5 trials versus model training time for the Navier-Stokes equations when the number of training samples is 1024.
B.6
Output reconstruction error
To measure the error induced by the projection, we define the output reconstruction emoor as follows
u lu((,m - u(m [u(m(()) )|1
where P. is the rank r linear projector. We provide the plots of the output reconstruction enror vs number reduced basis r using KLE basis and ASM basis in Figure 12. We can see that using the ASM basis results in lower output reconstruction error than the KLE basis.
Ouput reconstruction error
Ouput reconstruction error
KLE
KLE
10-1-
rel42-arT
rel-12-erd--t
10
107
10a rank
10a
107
10 rank
10"
Figure 12: Output reconstruction emro using KLE and ASM basis. Left: Hyperelasticity; Right: Navier-Stokes
B.7 Visualization of the ground truth and prediction
We present the comparisons of the ground truth of solution model prediction à using different methods, and absolute value of their difference (u(x) ū(z) in Figures 13 and 14 for the hypere- lasticity equation, and Figures 17 and 18 for the Navier-Stooe Stokes equations. In addition, we present the comparisons of the ground truth of the derivative of u with respect to m the direction of wi, denoted as du(m; w1), model prediction dü(m;; g/()) using different methods, and absolute value of their difference Idu(m;; w1) dū(m;:)) in Figures 15 and 16 for the hyperelasticity equation, and Figures 19 and 20 for the Navier-Stoken equations.
23

---

Grand kus
Fraioton
Fradotion
Pradoces
19
Prefetie
Protettn
Frotetden
Frotetter
4.00
14
05
cc0
44
20
44
2411
cooe
44
Koo
c001
4:0
408
42-
Proferde
Proderle
65
43
Figure 13: Hyperelasticity. Comparison of the predictions of #1 with (16, 64, 256, 1024) training samples using different methods. Top row -> bottom row: DeepONet, FNO, DE-FNO (Random), DINO, DE-DeepONet (KLE), DE-DeepONet (ASM).
24

---

100
100
64
34
24
04
1:09
400
Arand Tue
PeC 13
Pee
Ped 13
Figure 14: Hyperelasticity. Comparison of the predictions of u2 with (16, 64, 256, 1024) training samples using different methods. Top row -> bottom row: DeepONet, FNO, DE-FNO (Random), DINO, DE-DeepONet (KLE), DE-DeepONet (ASM).
25

---

30
955
Pow
01
81
22
cose
34
09
64
4.30
000
KMC
Figure 15: Hyperelasticity. Comparison of the predictions of directional derivative du, (m; w1) with (16, 64, 256, 1024) training samples using different methods. Top row -> bottom row: DeepONet, FNO, DE-FNO (Random), DINO, DE-DeepONet (KLE), DE-DeepONet (ASM).
26

---

Cmcond
6400
Soonnd Push
1007
Endicdos
A401
65
100
341
Figure 16: Hyperelasticity. Comparison of the predictions of directional derivative duz(m; wi) with (16, 64, 256, 1024) training samples using different methods. Top row -> bottom row: DeepONet, FNO, DE-FNO (Random), DINO, DE-DeepONet (KLE), DE-DeepONet (ASM).
27

---

Doure Tun
Recottn
Ireteres
Teeletes
65
Drund
50
01
05
M
08
M4
14
34
Figure 17: Navier-Stokes. Comparison of the predictions of veloeity-x with (16, 64, 256, 1024) training samples using diffferent methods. Top row -> bottom row: DeepONet, FNO, DE-FNO (Random), DINO, DE-DeepONet (KLE), DE-DeepONet (ASM).
28

---

Ned
Nedila
Noti
Figure 18: Navier-Stokes. Comparison of the predictions of velocity-y with (16, 64, 256, 1024) training samples using different methods. Top row bottom row: DeepONer, FNO, DE-FNO (Random), DINO, DE-DeepONet (KLE), DE-DeepONet (ASM).
29

---

33
61
123
140
103
Fablen
-20
Gruud Eve
Prddoe
40
103
6.6
10
05
Redtiies
Redetie
05
10
82
05
10
Figure 19: Navier-Stokes. Comparison of the predictions of directional derivative du, (m; wn) with (16, 64, 256, 1024) training samples using different methods. Top bottom row: DeepONet, FNO, DE-FNO (Random), DINO, DE-DeepONet (KLE), DE-DeepONet (ASM).
30

---

213
013
310
314
3.0
330
3.0
Scord FuEY
Pedictar
01
45
Figure 20: Navier-Stokes. Comparison of the predictions of directional derivative (m; w1) with (16, 64, 256, 1024) training samples using different methods. Top row -> bottom row: DeepONet, FNO, DE-FNO (Random), DINO, DE-DeepONer (KLE), DE- DE-DeepONet (ASM).
31

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: We include paper's contributions in the abstract and the second paragraph introduction and its scope in the first and third paragraphs of introduction.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification:
We discuss the limitations in the last section.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g.. independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect 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 is low images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
the paper. The authors should use their best
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]
32

---

Justification:
We provide one theorem with clearly stated assumptions and proof.
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 the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix 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: We discuss all details (including data generation and model training) that help reproduce the experiments in the appendix.
Guidelines:
The answer NA means that the paper does include experiments. If the paper includes experiments, No answer to this question will not be 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 to 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 is a specific model and empirical evaluation, it may be necessary to either make possible for others to 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 are 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) If the contribution is primarily a new algorithm, the paper should make it clear how reproduce that algorithm. (b) If the contribution is primarily
model architecture, the paper should describe the architecture clearly and fully.
(c) If 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., registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
33

---

Answer: [Yes]
Justification: We provide the code necessary to generate the data and train and test the different models considered in the paper.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips cc/ public/guides/Cdeen SoweiddeeeooeeSuimiiiiinnoneeee
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to to reproduce the results. See the NeurIPS code and data submission guidelines (https: httpsssccppubllc//ee Seeguddee/ooeeeeeeeeeeee
ssi agionPolcyy) for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, 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
Question: 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?
Answer: [Yes]
Justification:
We discuss the experimental details in the Appendix B.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, or as supplemental material.
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]
Justification: We run five trails with different random seed for each method and report the mean value and standard deviation of the corresponding relative error.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by emror bars, confi- dence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.
The factors of variability that the error bars capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall 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).
34

---

It should be clear whether the error bar is the standard deviation or the standard error of the mean.
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 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 out of 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 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 provide the machine information and time of exccution in the Tables 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 compute required for each of the individual experimental runs 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 't make it into the paper).
9. Code Of Ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics fttts:///neure ips. cc/ htcp/publcc//hhcscualliiioeeeeee ines?
Answer: [Yes]
Justification:
We have reviewed the code of ethies and checked that the research follows
them.
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
If 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 (e.g., if there is special consid- eration due to laws or regulations 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.
35

---

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 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 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 ML).
11. Safeguards
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, scraped datasets)?
Question:
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 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 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 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 dataset.
The authors should state which version of the asset used and, if possible, include URL.
The name of the license (c.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, paperswi thcode. com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of 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.
36

---

If this information is 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:
We do not release new assets at the moment.
Guidelines:
The answer NA means that the paper does not release new assets. Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification: Our research is in the field of the scientific machine learning and does not involve crowdsourcing or human subjects.
Guidelines:
The answer NA means that the paper does 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 much detail as possible should be included the main paper.
According to the NeurIPS Code Ethics, workers involved in data collection, curation, 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 an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: Our research is the field of the scientific machine learning and does not involve crowdsoureing human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor 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.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
37


Reviews:
[{"review_id": "I8gBUmqxS2", "reviewer": "Reviewer_row2", "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.", "questions": "- Could you provide a higher-level explanation of how the derivative labels are calculated? Are they derived from the underlying equation?\n- 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?\n- Why does supervision of the derivatives become less significant as more data becomes available?", "limitations": "Limiatations have been addressed", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "- The paper is engaging and well-motivated.\n- The methodology presented appears novel and is supported by solid theoretical underpinnings. It demonstrates consistent enhancements across various metrics, particularly with the hyperelasticity dataset.", "weaknesses": "- 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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your recognition of our work. We'd like to address the questions:\n\n1. The compution of derivative labels $p=du(m;\\psi)$ is equivalent to solving Eq.(13). The reviewer can find more details in the former part of the proof Theorem 1. The discretized derivative labels $p(x_i)$ are obtained by evaluating p on the grid points $x_i$. Note that these labels are a little different from the derivative labels used in DINO. In the revised version, we will add more details about the computation of our labels and DINO's labels.\n\n2. Although we currently do not have extensive experiments in this direction of improvement, we provide a small test. We train DE-DeepONet using dm loss only at the later 10% epoch of total epochs (denoted as DE-DeepONet later) on Navier--Stokes equations. The following are comparisions of the test relative error when using 16 samples to train DE-DeepONet/DeepONet for 32768 iterations.\n\nmethod            | rel-L2-err   | rel-H1-err | rel-Fro-err |\n------------------|   -----------| -----------| ------------|\nDE-DeepONet (ASM) |   9.10 %     |   13.01 %   |     39.76 %  |\nDE-DeepONet later (ASM) |   9.49 % |        13.45 %     |     40.87 %  |\nDE-DeepONet (KLE) |   13.75 %        |     18.03 %  |   78.70 %    |\nDE-DeepONet later (KLE) |   19.11 % |       25.15 %      |       136.30 %   |\nDeepONet                 |  27.88 %    |    26.70 %         |     130.79 %   |\n\nWe can see that DE-DeepONet later (ASM) has almost the same accuracy as DE-DeepONet (ASM). However, when replacing the ASM basis with the KLE basis, the significant increase in relative error makes this technique less appealing -- it may be necessary to train with dm labels for more epochs.\n\n3. Intuitively, derivatives can be viewed as sensitivity information for the neighborhood of each input-output pair. In other words, if we slightly perturb the input, we can predict how the output will change accordingly. It's akin to enhancing the training dataset from (input point, output point) to (input neighborhood, output neighborhood), where each point represents a function. When more (input point, output point) pairs are available, the neighborhood information can be revealed by the newly added nearby points, making the derivatives less useful.", "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": "GbhbvHHW2v", "reviewer": "Reviewer_41ML", "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.", "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? \n\n- Can the authors describe how the Frechet derivative is computed (the derivative with respect to the reduced input) to collect the ground truth data? \n\n- Is the reason for using the relative error for measuring the loss for preventing numerical issues caused by different numerical scales?", "limitations": "- A thorough empirical investigation on in which scenarios this proposed method outperforms baseline (and also the opposite scenarios) would be needed.", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "- The paper is written clearly, explaining the relevant preliminaries in a concise way and delivering the core contributions of the proposed method.\n\n- The proposed method leverages well-established classical methods (such as KLE and ASM) for building low-dimensional surrogate and incorporate them into advanced neural operator settings (DeepONet). \n\n- Although not very extensive, the paper present a set of numerical results from experimentation on important benchmark problems.", "weaknesses": "- 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. \n\n- 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. \n\n- 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. \n\n- 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. \n\n- 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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate your thoughtful comments. Regarding the weaknesses: \n\n1. We updated the experimental results as shown in the general response. Our method have much lower test relative errors compared to FNO when the training samples are scarce (N_train = 16 or 64). The main reason our method does not perform well with a large number of traning samples (N_train = 256 or 1024) is due to input dimension reduction, where the reduction error becomes the dominant component of the total error. This can be seen from the fact that relative L2 error (when N_train is large) is quite close to the output construction error at rank = 16 shown in Figure 6 in the appendix. Moreover, we show that the additional dm supervision can also be applied to FNO to enhance its prediction of solution and functional derivatives. \n\n2. The gains and limitations of our method indeed inherits from the linear reduction techniques such as KLE and ASM. We will include a more detailed discussion about this in the revised version. \n\n3. We show the plots of test errors vs training time in the pdf file in the general response. \n\n4. Thanks for your suggestion. We provide one table showing the inference time of each method and one table showing the data generation with different number of training samples. \n\n5. Figure 2 may be a bit misleading. The color of the output indicates the magnitude to the displacement $u$ (which maps from domain $\\Omega$ to $\\mathbb{R}^2$), rather than any the component of $u$. And the skewed square shows locations of the any domain point after deformation. This is one of the commonly used visualization methods in studying elasticity problems. To see the componentwise function $u_1$ and $u_2$, please refer to Figure 7 and 8 in the appendix. \n\nRegarding the questions: \n\n1. In Eq(4), the branch net takes the reduced parameters as inputs. If the branch net takes original input, we believe it is not possible to compute Gateaux derivatives (which we actually use) only using reduced inputs by automatic differentiation (AD). But it is possible to compute the Gateaux derivatives of original inputs using AD, and match the these directives outputs with the same labels used in this work. This approach is similar to using dm labels in FNO.\n\n2. The high level of computation of Gateaux derivative functions $p=du(m;\\psi)$ shown in Eq. (13) can be found in the former part of the proof of Theorem 1. The discretized derivative labels are obtained by evaluating function $p$ on the grid points $x_i$. Note that these labels are a little different from the derivative labels used in DINO. In the revised version, we will add more details about the computation of our labels and DINO's labels.\n\n3. There are two main reasons we use relative error for measuring the loss. First, we find that it aligns the magnitudes of different loss terms (i.e., evaluation loss and derivative losses), with the loss weights lambda_i also being of the same mangnitude. Intuitively, this allows the network to learn multiple objectives more effectively. Our experiments show that in our problem setups, relative L2 error performs better than MSE. Second, our evaluation metrics use relative error instead of absolute error (which is fairer for comparision of ground truth and prediction), it is natural to train the neural network to optimize according to the same metric (with minimal additional computational cost compared to MSE).", "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": "wmxOe1KTQj", "reviewer": "Reviewer_WTPz", "summary": "This manuscript introduces a derivative-enhanced deep operator network that utilizes derivative information to improve prediction accuracy.", "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]?\n\n[1] DEEP MULTI-SCALE VIDEO PREDICTION BEYOND MEAN SQUARE ERROR. Michael Mathieu, Camille Couprie & Yann LeCun.\n\nQ2) What is the impact of proposed losses when combined with other operators like FNO? Do they enhance the performance of different operators?\n\nQ3) Could approximating the Frobenius norm along a random direction be analogous to score-based diffusion models?\n\nQ4) 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?\n\nQ5) How does incorporating this noisy gradient version aid in training?", "limitations": "Yes, the authors adequately addressed the limitations.", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 3, "strengths": "S1) This paper presents a new method for improving the accuracy of approximating the output function for DeepONet, along with its directional derivative relative to the input function and spatial derivatives.\n\nS2) The suggested enhancements improve the performance of the basic DeepONet and perform well in scenarios with limited training samples.", "weaknesses": "W1) The baselines appear limited as recent benchmarks like GNOT, LSM, ONO, Transsolver, etc., have ye\nt to be compared. Numerous variations of FNO have been introduced, such as FFNO, UNO, CoNO, GFNO, UFNO, etc.\n\nW2) The benchmarks dataset appears to be restricted to datasets for Navier-Stokes equations and hyperelasticity, which are not open-source.\n\nW3) ASM is costly in computation, and derivative information is necessary when the PDE form is unknown. Additionally, it needs to scale better with dimensions.\n\nW4) The proposed method builds upon the DINO paper by introducing additional informed losses.\n\nW5) The new losses were found effective only for DeepOnet and lacked experimental evidence compared to FNO.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thoughtful comments. Regarding the weaknesses: \n\n1. Most of the benchmarks you mention are improved by changing the architecture (or essentially the parameterization method) rather than focusing on adding regularization terms to the loss function, which is the focus of our paper. Also, these benchmarks primarily focus on output solution accuracy. While it is possible to modify the existing code to support dm prediction and dm test error evaluation, it would require a significant amount of effort and time. Thanks for your understanding. \n\n2. We have not found any public datasets containing dm labels, and these labels cannot be directly computed from input-output pairs; they need to be computed by solving a linear PDE induced by the original PDE. Therefore, we have to compute dm labels by ourselves. We use FEniCS to automate the generation of dm labels, thanks to its supporting for Gateaux derivative. The code will be made available upon publication.\n\n3. Computing ASM basis is approximately equivalent to solving O(N_{grad}) linear PDEs (regardless the original PDEs is nonlinear or not) together with a generlized eignvalue problem, where N_{grad} is the number of input-output pairs used for the Monte Carlo estimation. As shown in Fig 3 \\& 5 in [1], at least in some cases, a small number of samples is sufficient to compute the basis accurately. In our problem setups, we find that setting N_{grad} to 16 already provides a satifactory approximation. Also, we use the double pass randomized algorithm when solving the generalized eigenvalue problem to improve scalability. \n\n[1] Zahm, Oliver, et al. \"Gradient-based dimension reduction of multivariate vector-valued functions\"\n\n4. The arguments in the general response provide a more detailed explanation of how DE-DeepONet differs from DINO.\n\n5. We add experiments to show that adding dm regularization also significantly improve FNO's accuracy. We believe this idea can be applied to other variants that have the same two ends (i.e., input and output) of networks. \n\nRegarding the questions: \n\n1. For the approximation for spatial derivative, [2] use finite difference to compute both spatial derivative outputs and labels. We use PyTorch's AD to compute outputs and FEniCS's builtin gradient function to compute labels. \n\n[2] Mathieu, Michael, et al. \"Deep multi-scale video prediction beyond mean square error\"\n\n2. The results in general response show that dm loss greatly enhance FNO's prediction accuracy.\n\n3. We observe a certain similarity between our dm loss formulation and the loss function in score-based diffusion models. However, a key difference is that in our work, we use the same neural network to approximate both the solution and its derivatives with respect to input function, whereas in score-based diffusion models, the neural network is trained solely to approximate the derivatives. Also, we estimate the derivative of $u(x_i)$ (where x_i are grid points), which differs from estimate the derivative of log of conditional probability in score-bsaed diffusion models.\n\n5. 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 traning samples) whereas in FNO we use full inputs.\n\n6. We do not add noise to any gradients/derivatives involved (including dm, dx, and the gradient of trainable parameters). We are not sure that we fully understand your question. Would you mind explaining your question more clearly? Thank you.", "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": "V9O8t9r7yv", "reviewer": "Reviewer_nyqb", "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.", "questions": "-\tCould 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. \n\n-\tHow 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. \n\n-\tIs 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?\n\n-\tThe 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?\n\n-\tThe 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. \n\n[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.\n\n[3] Wang, Sifan, et al. \"Learning the solution operator of parametric partial differential equations with physics-informed DeepONets.\" Science advances 7.40 (2021): eabi8605.\n\n-\tHow 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. \n\n[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": "-\tAdequately described in manuscript.", "rating": 7, "confidence": 5, "soundness": 3, "presentation": 4, "contribution": 2, "strengths": "- The paper is very well written, well structured, and with the right level of detail and understanding to describe the setup in the preliminary section without lingering too long. The authors clearly have a deep understanding of all components of this method and convey the steps well in the space given. \n\n- The appendix is thorough, high quality, and provides relevant details to better comprehend and reproduce the results in the main text. I found the visualizations in section B.6 particularly compelling for the low data case with vanilla FNO and DeepONet failing compared to the gradient enhanced models.\n \n- The improvement in dm prediction for control is an important but overlooked metric for neural operators which primary focus on the output solution accuracy despite the gradient information being needed in real-world problems. The authors may want to consider the existing real-world example of DeepONets for optimization in [1], shown in Figure 7 to require gradient information of the objective w.r.t to the input parameters to optimize aerodynamic shapes which could benefit from this methodology which is a strength over accuracy improvement alone. \n\n[1] Shukla, Khemraj, et al. \"Deep neural operators as accurate surrogates for shape optimization.\" Engineering Applications of Artificial Intelligence 129 (2024): 107615.", "weaknesses": "-\tThe 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. \n\n-\tI 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. \n\n-\tFollowing 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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your recognition of our work and for your constructive feedback. They are valuable to us.\n\n1. We address the concerns on computational cost and dfference between our method and DINO in the general response. Additionally, we provide the convergence plots in the pdf file.  We hope the updated results can better help evaluate our work. \n\n2. According to ML community, iterations and epochs are different concepts, even though sometimes they can be equal. One epoch is a complete pass of training data while an iteration is a single update to the model's parameters using a batch of data (typically not the whole training data). Therefore, one epoch equals to ceil(N_{train}/batch size) iterations. In our work, we set batch size 8. When N_train = 1024, epochs = 1000, the model's parameters update (1024/8)*1000=128,000 times, which lies in the range of 50,000 - 500,000 the original DeepONet paper used. To make a more fair comparision, we train all models using different number training samples for the same number of iterations (32768) instead of epochs (1000), and evaluate them on test dataset at milestones: 128, 512, 2048, 8192, 32768 and compute the three test relative errors.\n\n3. The raw inputs are the same for all models, while in vanilla DeepONet and FNO we feed the network with raw inputs (no dimensionality reduction), and in DE-DeepONet and DINO, we feed the network with post-processed inputs (that is, the reduced inputs by projecting raw inputs into low dimensional linear subspace spanned by KLE or ASM basis).\n\n4. Thanks for the suggestion. Actually we use CNN for the branch (which receives high dimensional vectors m) and ResNet for the trunk (which receives two dimensional spatical coordiates x). We find that using CNN for the branch yields lower rel-L2-err than using ResNet. Some works also consider this choice (see e.g., Table C.1 in [1]). Since in our problem setup, discretizated parameter m on 2D domain dcan be viewed as a image, it is quite natural to consider using neural network architectures that are originally designed to solve image classification tasks. \n\n[1] Lu, Lu, et al. \"A comprehensive and fair comparison of two neural operators (with practical extensions) based on FAIR data\"\n\n5. Thanks for reminding us of the closely related work PI-DeepONet and gPINN. We notice that our dm loss can be directly incorporated into PI-DeepONet's branch net without influencing its trunk net. This could possibly further reduce the requirement of training data. Similar to gPINN, it is possible to enforce the outputs obey not only the residual equal to zero but also the directional derivative of the residual w.r.t. m equal to zero. But this method can be potentially more difficult to implement in practice since the Gateaux derivative of residual is more complex than residual. \n\n6. We use PyTorch to do AD. We thank the reviewer for mentioning Jax, which we will consider in our later work.", "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 are sincerely grateful to the reviewers for dedicating their time and effort to reviewing our work and providing helpful feedback. 

Compared to DINO, although the DeepONet architecture (and its formulation of dm loss) requires longer training time, it offers the following advantages

- Much shorter inference time of $du(m;\psi)(x)$. The additional trunk net (which receives spatial coordinates) allows us to quickly query the sensitivity of output function $u$ at any point $x$ when input function $m$ is perturbed in any direction $\psi$. 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. We provide more details in the revised version. 

- 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, no matter how we modify the underlying ResNet/MLP and/or loss function. 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 [1] and replacing the  inner product of the outputs of two networks by more flexible operations [2] [3] (to enhance expressive power). The dm loss formulation of our work differs nontrivally by DINO's, but it is broadly suitable any network architecture that has multiple subnetworks, where at least one of them receives high-dimensional inputs.

[1] Wang, Sifan, et al. "Learning the solution operator of parametric partial differential equations with physics-informed DeepONets."

[2] Pan, Shaowu, et al. "Neural Implicit Flow: a mesh-agnostic dimensionality reduction paradigm of spatio-temporal data"

[3] Hao, Zhongkai, et al.  "GNOT: A General Neural Operator Transformer for Operator Learning" 

We made an update to the experimental results. The modifications include 

1. Adding three models -- FNO trained with dm loss along ASM, KLE and Random (sampled from same Gaussian random field of parameter m) directions -- into comparision.
2. Removing dx loss in training. Our experiments show that dx labels typically do not help improve model's prediction accuracy if dm labels are used in training.
3. Fixing iterations (32,768) instead of epochs (1000) for model training, no matter which model or how many training samples are used. We compute the test relative errors at iteration milestones: 128, 512, 2048, 8192, 32768.  
4. Adjusting the hyperparameters of DE-DeepONet and DINO. We double their widths and halving their depths. We disable the learning rate scheduler step_LR. For the DE-DeepONet, we disable dx loss and add Fourier features embedding (Gaussian RFF mapping) [4] into trunk net.

[4] Tancik, Matthew, et al. "Fourier features let networks learn high frequency functions in low dimensional domains"

In the file **rebuttal_rel_err.pdf**, we provide eight new plots of relative test errors on 500 test samples on Navier--Stokes equations in L2 norm and Fro norm (H1 norm is similar to L2 norm, we omit it here due to space limitation) versus different number of iterations and total training time when using 16 or 256 training samples (representative of limited data and sufficient data scenarios). We can see that when N_train=16, DE-DeepONet (ASM) and DINO have much lower error compared to FNO, FNO+dm and DeepONet. However, given 256 training samples and sufficient training time, FNO+dm achieves the lowest test error. The choice of perturbation directions (ASM, KLE, and Random) does not make a significant difference. Furthermore, FNO+dm consistently outperforms the vanilla FNO, although it requires much longer training time. 

We also provide the total wall clock time of each model inferencing on **500** test samples of **both solution and dm in 128 random directions** on a single GPU and a single CPU (which is needed for post processing of outputs in DE-DeepONet and DINO). Note that the time for predicting dm dominates the total inference time. 

| Model | Inference time (seconds) | 
|-|-|
|FNO|53| 
|DeepONet|4|  
|DE-DeepONet|47|  
|DINO|415| 
|Numerical solver (16 CPUs)|1103|

The reason why DINO has much more inference time compared to DE-DeepONet is due to the post processing the reduced dm requires very large matrices multiplications (one of the dimensions is the degree of freedom of high fedility solution, which in our case is 66050). And this is even an overoptimistic estimate since we exclude the time for computing a repeated used large matrix, i.e., the output finite element basis functions evaluating on gird points and its matrix multiplication with the nodal values of output reduced basis. 

Finally, we show the total wall clock time of data generation of our DE-DeepONet (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.

| N_train | Data generation (s) | 
|-|-|
|16|17| 
|64| 38|  
|256|125|  
|1024|470| 

We hope the above figures and tables provide readers with a more clear view of the offline and online cost vs accuracy.

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