DeepONet for Solving PDEs: Generalization Analysis in Sobolev Training

The rest of the paper is organized as follows. In Section 2, we introduce the notations used throughout the paper and set up the operator and loss functions from the problem statement. In Section 3, we estimate the approximation rate. In the first subsection, we reduce operator learning to function learning, showing that the trunk network does not benefit from a deep structure. In the second subsection, we reduce functional approximation to function approximation and demonstrate that the encoding of the input space for the operator can still rely on traditional sampling methods, without the need to incorporate derivative information, even in Sobolev training. In the last subsection, we approximate the function and show that the branch network can benefit from a deep structure. In Section 4, we provide the generalization error analysis of DeepONet in Sobolev training. All omitted proofs from the main text are presented in the appendix.

Function Learning: Unlike operator learning, function learning focuses on learning a map from a finite-dimensional space to a finite-dimensional space. Many papers have studied the theory of function approximation in Sobolev spaces, such as [27, 30, 31, 32, 33, 42, 40, 39, 28, 44, 45], among others. In [27, 39, 45], the approximation is considered under the assumption of continuous approximators, as mentioned in [8]. However, for deep neural networks, the approximator is often discontinuous, as discussed in [45], which still achieves the optimal approximation rate with respect to the number of parameters. This is one of the benefits of deep structures, and whether this deep structure provides a significant advantage in approximation is a key point of discussion in this work.

Functional Learning: Functional learning can be viewed as a special case of operator learning, where the objective is to map a function to a constant. Several papers have provided theoretical analyses of using neural networks to approximate functionals. In [6], the universal approximation of neural networks for approximating continuous functionals is established. In [35], the authors provide the approximation rate for functionals, while in [41], they address the approximation of functionals without the curse of dimensionality using Barron space methods. In [43], the approximation rate for functionals on spheres is analyzed.

Operator Learning: The operator learning aspect is closely related to this work. In [18, 7], they perform an error analysis of DeepONet; however, their analysis is primarily focused on solving PDEs in the $L^{2}$-norm, and they do not consider network design, i.e., which parts of the network should be deeper and which should be shallow. In [17], the focus is on addressing the curse of dimensionality and providing lower bounds on the number of parameters required to achieve the target rate. In their results, there is a parameter $\gamma$ that obscures the structural information of the network, so the benefit of network architecture in operator learning remains unclear. In [22], they summarize the general theory of operator learning, but the existence of the benefits of deep neural networks is still not explicitly visible in their results. Other works, such as those focusing on Fourier Neural Operators (FNO) [16], explore different types of operators, which differ from the focus of this paper.

Let us summarize all basic notations used in the paper as follows:

1. Matrices are denoted by bold uppercase letters. For example, $\bm{A}\in\mathbb{R}^{m\times n}$ is a real matrix of size $m\times n$ and $\bm{A}^{\intercal}$ denotes the transpose of $\bm{A}$.

2. For a $d$-dimensional multi-index $\bm{\alpha}=[\alpha_{1},\alpha_{2},\cdots\alpha_{d}]\in\mathbb{N}^{d}$, we denote several related notations as follows: $(a)~{}|\bm{\alpha}|=\left|\alpha_{1}\right|+\left|\alpha_{2}\right|+\cdots+\left|\alpha_{d}\right|$; $(b)~{}\bm{x}^{\alpha}=x_{1}^{\alpha_{1}}x_{2}^{\alpha_{2}}\cdots x_{d}^{\alpha_{d}},~{}\bm{x}=\left[x_{1},x_{2},\cdots,x_{d}\right]^{\intercal}$; $(c)~{}\bm{\alpha}!=\alpha_{1}!\alpha_{2}!\cdots\alpha_{d}!.$

3. Let $B_{r,|\cdot|}(\bm{x})\subset\mathbb{R}^{d}$ be the closed ball with a center $\bm{x}\in\mathbb{R}^{d}$ and a radius $r$ measured by the Euclidean distance. Similarly, $B_{r,\|\cdot\|_{\ell_{\infty}}}(\bm{x})\subset\mathbb{R}^{d}$ be the closed ball with a center $\bm{x}\in\mathbb{R}^{d}$ and a radius $r$ measured by the $\ell_{\infty}$-norm.

4. Assume $\bm{n}\in\mathbb{N}_{+}^{n}$, then $f(\bm{n})=\bm{O}(g(\bm{n}))$ means that there exists positive $C$ independent of $\bm{n},f,g$ such that $f(\bm{n})\leq Cg(\bm{n})$ when all entries of $\bm{n}$ go to $+\infty$.

5. Define $\sigma_{1}(x):=\sigma(x)=\max\{0,x\}$ and $\sigma_{2}:=\sigma^{2}(x)$. We call the neural networks with activation function $\sigma_{t}$ with $t\leq i$ as $\sigma_{i}$ neural networks ($\sigma_{i}$-NNs). With the abuse of notations, we define $\sigma_{i}:\mathbb{R}^{d}\to\mathbb{R}^{d}$ as $\sigma_{i}(\bm{x})=\left[\begin{array}[]{c}\sigma_{i}(x_{1})\\ \vdots\\ \sigma_{i}(x_{d})\end{array}\right]$ for any $\bm{x}=\left[x_{1},\cdots,x_{d}\right]^{T}\in\mathbb{R}^{d}$.

6. Define $L,W\in\mathbb{N}_{+}$, $W_{0}=d$ and $W_{L+1}=1$, $W_{i}\in\mathbb{N}_{+}$ for $i=1,2,\ldots,L$, then a $\sigma_{i}$-NN $\phi$ with the width $W$ and depth $L$ can be described as follows:

$$\bm{x}=\tilde{\bm{h}}_{0}\stackrel{{\scriptstyle W_{1},b_{1}}}{{\longrightarrow}}\bm{h}_{1}\stackrel{{\scriptstyle\sigma_{i}}}{{\longrightarrow}}\tilde{\bm{h}}_{1}\ldots\tilde{\bm{h}}_{L}\stackrel{{\scriptstyle W_{L+1},b_{L+1}}}{{\longrightarrow}}\phi(\bm{x})=\bm{h}_{L+1},$$

where $\bm{W}_{i}\in\mathbb{R}^{W_{i}\times W_{i-1}}$ and $\bm{b}_{i}\in\mathbb{R}^{W_{i}}$ are the weight matrix and the bias vector in the $i$-th linear transform in $\phi$, respectively, i.e., $\bm{h}_{i}:=\bm{W}_{i}\tilde{\bm{h}}_{i-1}+\bm{b}_{i},~{}\text{ for }i=1,\ldots,L+1$ and $\tilde{\bm{h}}_{i}=\sigma_{i}\left(\bm{h}_{i}\right),\text{ for }i=1,\ldots,L.$ In this paper, an DNN with the width $W$ and depth $L$, means (a) The maximum width of this DNN for all hidden layers less than or equal to $W$. (b) The number of hidden layers of this DNN less than or equal to $L$.

7. Denote $\Omega$ as $[0,1]^{d}$, $D$ as the weak derivative of a single variable function and $D^{\bm{\alpha}}=D^{\alpha_{1}}_{1}D^{\alpha_{2}}_{2}\ldots D^{\alpha_{d}}_{d}$ as the partial derivative where $\bm{\alpha}=[\alpha_{1},\alpha_{2},\ldots,\alpha_{d}]^{T}$ and $D_{i}$ is the derivative in the $i$-th variable. Let $n\in\mathbb{N}$ and $1\leq p\leq\infty$. Then we define Sobolev spaces

$$W^{n,p}(\Omega):=\left\{f\in L^{p}(\Omega):D^{\bm{\alpha}}f\in L^{p}(\Omega)\text{ with }|\bm{\alpha}|\leq n\right\}$$

with a norm

$$\|f\|_{W^{n,p}(\Omega)}:=\left(\sum_{0\leq|\alpha|\leq n}\left\|D^{\bm{\alpha}}f\right\|_{L^{p}(\Omega)}^{p}\right)^{1/p},~{}p<\infty,$$

and $\|f\|_{W^{n,\infty}(\Omega)}:=\max_{0\leq|\alpha|\leq n}\left\|D^{\bm{\alpha}}f\right\|_{L^{\infty}(\Omega)}$. For simplicity, $\|f\|_{W^{n,2}(\Omega)}=\|f\|_{H^{n}(\Omega)}$. Furthermore, for $\bm{f}=(f_{1},f_{2},\ldots,f_{d})$, $\bm{f}\in W^{1,\infty}(\Omega,\mathbb{R}^{d})$ if and only if $f_{i}\in W^{1,\infty}(\Omega)$ for each $i=1,2,\ldots,d$ and $\|\bm{f}\|_{W^{1,\infty}(\Omega,\mathbb{R}^{d})}:=\max_{i=1,\ldots,d}\{\|f_{i}\|_{W^{1,\infty}(\Omega)}\}$.

In this paper, we consider applying DeepONet to solve partial differential equations (PDEs). Without loss of generality, we focus on the following PDEs (3), though our results can be easily generalized to more complex cases. In this paper, we assume that $\mathcal{L}$ satisfies the following conditions:

The above assumption holds in many cases. For example, based on [12], when $\mathcal{L}$ is a linear elliptic operator with smooth coefficients and $\Omega$ is a polygon (which is valid for $[0,1]^{d}$), condition (i) in Assumption 1 is satisfied. Regarding condition (ii) in Assumption 1, this assumption reflects that the operator we aim to learn is Lipschitz in $W^{n}(\Omega)$, which is a regularity assumption. If the domain is smooth and $\mathcal{L}$ is linear and has smooth coefficients, then $\|u\|_{W^{n+2,\infty}(\Omega)}\leq L\|f\|_{W^{n,\infty}(\Omega)}$ [11], which is a stronger condition than the one we assume here.

Based on above assumption, the equation above has a unique solution, establishing a mapping between the source term $f\in\mathcal{X}$ and the solution $u\in\mathcal{Y}$, denoted as $\mathcal{G}_{*}$. In [24], training such an operator is treated as a black-box approach, i.e., a data-driven method. However, in many cases, there is insufficient data for training, or there may be no data at all. Thus, it becomes crucial to incorporate information from the PDEs into the training process, as demonstrated in [37, 21, 13]. This technique is referred to as Sobolev training.

In other words, we aim to use $\mathcal{G}(~{}\cdot~{};\bm{\theta})$ (see Eq. (2)) to approximate the mapping between $f$ and $u$. The chosen loss function is:

	$\displaystyle L_{D}(\bm{\theta})$	$\displaystyle:=L_{\text{in}}(\bm{\theta})+\lambda L_{\text{bound}}(\bm{\theta}),$
	$\displaystyle L_{\text{in}}(\bm{\theta})$	$\displaystyle:=\int_{\Omega}\int_{\mathcal{X}}|\mathcal{L}\mathcal{G}(f;\bm{\theta})(\bm{y})-f(\bm{y})|^{2}\,\mathrm{d}\bm{y}\,\mathrm{d}\mu_{\mathcal{X}},$
	$\displaystyle L_{\text{bound}}(\bm{\theta})$	$\displaystyle:=\int_{\mathcal{X}}\int_{\partial\Omega}|\mathcal{G}(f;\bm{\theta})(\bm{y})|^{2}\,\mathrm{d}\bm{y}\,\mathrm{d}\mu_{\mathcal{X}},$		(4)

where $\mu_{\mathcal{X}}$ is a measure on $\mathcal{X}$, and $\mathcal{X}$ is the domain of the operator, i.e., $f\in\mathcal{X}$. The constant $\lambda$ balances the boundary and interior terms. The error analysis of $L_{\text{bound}}(\bm{\theta})$ is the classical $L^{2}$ estimation. Since we focus on the Sobolev training part in this paper, we set $\lambda=0$ for simplicity. Due to (i) in Assumption 1, we obtain:

$\displaystyle L_{D}(\bm{\theta})\leq C\int_{\mathcal{X}}\|\mathcal{G}(f;\bm{\theta})(\bm{y})-\mathcal{G}_{*}(f)(\bm{y})\|^{2}_{H^{2}(\Omega)}\,\mathrm{d}\mu_{\mathcal{X}}.$

(5)

Furthermore, if we consider that $\mu_{\mathcal{X}}$ is a finite measure on $\mathcal{X}$ and, for simplicity, assume $\int_{\mathcal{X}}\mathrm{d}\mu_{\mathcal{X}}=1$, we have:

$\displaystyle L_{D}(\bm{\theta})\leq C\sup_{f\in\mathcal{X}}\|\mathcal{G}(f;\bm{\theta})(\bm{y})-\mathcal{G}_{*}(f)(\bm{y})\|^{2}_{H^{2}(\Omega)}.$

(6)

This error is called the approximation error. The other part of the error comes from the sample, which we refer to as the generalization error. To define this error, we first introduce the following loss, which is the discrete version of $L_{D}(\bm{\theta})$ for $\lambda=0$:

$$L_{S}(\bm{\theta}):=\frac{1}{MP}\sum_{i=1}^{M}\sum_{j=1}^{P}|\mathcal{L}\mathcal{G}(f_{i};\bm{\theta})(\bm{y}_{j})-f_{i}(\bm{y}_{j})|^{2},$$

(7)

where $\bm{y}_{j}$ is an i.i.d. uniform sample in $\Omega$, and $f_{i}$ is an i.i.d. sample based on $\mu_{\mathcal{X}}$.

Then, the generalization error is given by:

$\displaystyle\mathbb{E}L_{D}(\bm{\theta}_{S})\leq L_{D}(\bm{\theta}_{D})+\mathbb{E}\left[L_{D}(\bm{\theta}_{S})-L_{S}(\bm{\theta}_{S})\right],$

(8)

where $\bm{\theta}_{D}=\arg\min L_{D}(\bm{\theta})$ and $\bm{\theta}_{S}=\arg\min L_{S}(\bm{\theta})$. The expectation symbol is due to the randomness in sampling $\bm{y}$ and $f$. The generalization error is represented by the last term in Eq. (22).

The sketch of the proof of Proposition 1 is as follows. First, we divide $\Omega$ into $p$ parts and apply the Bramble–Hilbert Lemma [3, Lemma 4.3.8] to locally approximate $v$ on each part using polynomials. The coefficients of these polynomials can be regarded as continuous linear functionals in $W^{n,\infty}(\Omega)$. Next, we use a partition of unity to combine the local polynomial approximations. Finally, we represent both the polynomials and the partition of unity using neural networks.

Based on Proposition 1, we can reduce the operator learning problem to functional learning as follows:

		$\displaystyle\left\|\mathcal{G}_{*}(f)-\sum_{k=1}^{p}c_{k}(\mathcal{G}_{*}(f))\mathcal{T}(\bm{y};\bm{\theta}_{2,k})\right\|_{H^{2}(\Omega)}\leq Cp^{-\frac{n-2}{d}}\|\mathcal{G}_{*}(f)\|_{H^{n}(\Omega)}\leq Cp^{-\frac{n-2}{d}}\|\mathcal{G}_{*}(f)\|_{W^{n,\infty}(\Omega)}$
	$\displaystyle\leq$	$\displaystyle Cp^{-\frac{n-2}{d}}\|f\|_{W^{n,\infty}(\Omega)},$		(11)

where the last inequality follows from Assumption 1. Therefore, without considering the approximation of $c_{k}(\mathcal{G}_{*}(f))$, the approximation error can be bounded by $\sup_{f\in\mathcal{X}}Cp^{-\frac{n-2}{d}}\|f\|_{W^{n,\infty}(\Omega)}$. The next step is to consider the approximation of the functional $c_{k}(\mathcal{G}_{*}(f))$.

The approximation order achieved in Proposition 1 corresponds to the optimal approximation rate for neural networks in the continuous approximation case. If we consider deep neural networks, we can achieve an approximation rate of $p^{-\frac{2(n-2)}{d}}$ using the neural network $\mathcal{T}(\bm{y};\bm{\theta}_{2}(v))$. As discussed in [40], most of the parameters depend on $v$, leading to significant interaction between the trunk and branch networks, which makes the learning process more challenging. Furthermore, this approximator can be a discontinuous one, as highlighted in [45, 39], where the functionals $c_{k}$ may no longer be continuous, thereby invalidating subsequent proofs. Therefore, while it is possible to construct deep neural networks for trunk networks, we do not gain the benefits associated with deep networks, and the approximation process becomes more difficult to learn.

In this part, we reduce each functional $c_{k}(\mathcal{G}_{*}(f))$ to a function learning problem by applying a suitable projection $\mathcal{D}$ in Eq. (2) to map $u$ into a finite-dimensional space. To ensure that the error in this step remains small, we require $f-\mathcal{P}\circ\mathcal{D}f$ to be small for a continuous reconstruction operator $\mathcal{P}$, which will be defined later, as shown below:

		$\displaystyle|c_{k}(\mathcal{G}_{*}(f))-c_{k}(\mathcal{G}_{*}(\mathcal{P}\circ\mathcal{D}f))|\leq C\|\mathcal{G}_{*}(f)-\mathcal{G}_{*}(\mathcal{P}\circ\mathcal{D}f)\|_{W^{n,\infty}(\Omega)}$
	$\displaystyle\leq$	$\displaystyle C\|f-\mathcal{P}\circ\mathcal{D}f\|_{W^{n,\infty}(\Omega)}.$		(12)

The first inequality holds because $c_{k}$ is a bounded linear functional in $W^{n-1,\infty}(\Omega)$, and thus also a bounded linear functional in $W^{n,\infty}(\Omega)$. The second inequality follows from (ii) in Assumption 1.

For the traditional DeepONet [24], $\mathcal{D}$ is a sampling method, i.e.,

$$\mathcal{D}(f)=(f(\bm{x}_{1}),f(\bm{x}_{2}),\cdots,f(\bm{x}_{m})),$$

and then a proper basis is found such that $\mathcal{P}\circ\mathcal{D}(f)=\sum_{j=1}^{m}\alpha_{j}(\mathcal{D}(f))\phi_{j}(\bm{x})$. To achieve this reduction in $W^{n,\infty}(\Omega)$, one approach is to apply pseudo-spectral projection, as demonstrated in [26, 4] and used in the approximation analysis of Fourier Neural Operators (FNO) in [16].

Based on Proposition 2, we can reduce functional learning to function learning. Note that, even though we are considering Sobolev training for operator learning, we do not need derivative information in the encoding process. Specifically, the projection $\mathcal{D}$ defined as $\mathcal{D}(f)=(f(\bm{x}_{\bm{\nu}}))_{\bm{\nu}\in\{0,1,\ldots,2N\}^{d}}$ is sufficient for tasks like solving PDEs. Therefore, when designing operator neural networks for solving PDEs, it is not necessary to include derivative information of the input space in the domain. However, whether adding derivative information in the encoding helps improve the effectiveness of Sobolev training remains an open question, which we will address in future work.

In this part, we will use a neural network to approximate each functional $c_{k}(\mathcal{G}_{*}\circ\mathcal{P}\circ\mathcal{D}(f))$, i.e., we aim to make

$$\sup_{f\in\mathcal{X}}\left|c_{k}(\mathcal{G}_{*}\circ\mathcal{P}\circ\mathcal{D}(f))-\mathcal{B}(\mathcal{D}(f);\bm{\theta}_{1,k})\right|$$

(14)

as small as possible. For simplicity, we consider all functions in $\mathcal{X}$ with $\|f\|_{L^{\infty}(\Omega)}\leq M$ for some $M>0$. We then need to construct a neural network in $[-M,M]^{m}$ such that

$$\sup_{\bm{z}\in[-M,M]^{m}}\left|c_{k}(\mathcal{G}_{*}\circ\mathcal{P}(\bm{z}))-\mathcal{B}(\bm{z};\bm{\theta}_{1,k})\right|$$

(15)

is small. The idea for achieving this approximation is to first analyze the regularity of $c_{k}(\mathcal{G}_{*}\circ\mathcal{P}(\bm{z}))$ with respect to $\bm{z}$. Once the regularity is established, we can construct a neural network that approximates it with a rate that depends on the regularity.

Now, we can establish a neural network to approximate $c_{k}(\mathcal{G}_{*}\circ\mathcal{P}(\bm{z}))$. For shallow neural networks, the result from [27] is as follows:

For deep neural networks, the result from [32] is as follows:

Therefore, for shallow neural networks, based on Lemma 2, we know that using $(q)$ parameters can only achieve an approximation rate of $\bm{O}(q^{-\frac{1}{m}})$ for each $c_{k}(\mathcal{G}_{*}\circ\mathcal{P}(\bm{z}))$. However, for deep neural networks, the approximation rate can be $\bm{O}(q^{-\frac{\lambda}{m}})$, where $\lambda\in[1,2]$. Here, $\lambda=1$ corresponds to the shallow or very wide case (i.e., the depth is logarithmic in the width), and $\lambda=2$ corresponds to the very deep case (i.e., the width is logarithmic in the depth). This is because, in Lemma 3, a $\sigma_{1}$-NN with width $C_{1}W$ and depth $C_{2}L$ has a number of parameters $\bm{O}(W^{2}L)$, and the approximation rate is $\bm{O}\left(\left(W^{2}L^{2}\log^{3}(W+2)\right)^{-1/m}\right)$. When the depth $L$ is large and the width $W$ is fixed, the rate corresponds to $\lambda=2$. The case $\lambda\in(1,2)$ represents the middle ground, where the network is neither very deep nor very wide.

In [45, 39], it is mentioned that an approximation rate better than $\bm{O}(q^{-\frac{1}{m}})$ may cause the neural network approximator to become discontinuous. However, unlike the approximation of the trunk network, we do not require the approximation of the branch network to be continuous. Thus, we can establish a deep neural network for the approximation of the branch network and obtain the benefits of a deep structure, as indicated by the improved approximation rates.

The result in Proposition 3 exhibits a curse of dimensionality when $m$ is large. One potential approach to mitigate this issue is to incorporate a proper measure $\mu_{\mathcal{X}}$ instead of using $\sup_{f\in\mathcal{X}}$ in the error analysis. However, since this paper focuses on establishing the framework for generalization error in Sobolev training for operator learning and investigating whether deep structures offer benefits, addressing the curse of dimensionality in Sobolev training is left for future work. Solving this challenge would likely require $\mathcal{X}$ and $\mathcal{Y}$ to exhibit smoother structures.