A hybrid iterative method based on MIONet for PDEs: Theory and numerical examples¹¹1The research is supported by NSFC project 12288101.

Partial differential equations (PDEs) play a crucial role in describing physical reality and modeling engineering problems. As the vast majority of PDEs have no computable analytical solutions, seeking numerical solutions is an important issue. There have been many classical numerical methods developed for solving PDEs, such as the finite element method (FEM) [3], the finite difference method (FDM) [32], the finite volume method [6], and the spectral method [29]. With these methods, linear PDEs are usually discretized as linear algebraic systems, and subsequently be solved by direct methods (e.g. Gaussian elimination), or iterative methods such as Richardson, (damped) Jacobi, Gauss-Seidel, SOR, SSOR [33]. When the scale of the problem is very large, iterative methods show significant advantages like low memory requirements. However, a large scale system results in a long solving time, which is unacceptable. To accelerate the iteration, multilevel iterative methods are developed [8, 34, 38], a multigrid solver will greatly reduce the steps of iteration and lead to fast convergence. The shortcoming of the multigrid method is that it is mesh-depended, and requires users to generate several levels of meshes, which may be a very heavy task. In order to overcome this drawback, the algebraic multigrid (AMG) methods [2, 22, 28, 39] are proposed, which deal with the algebraic system directly without geometric meshes. To further promote the development of fast and low usage cost iterative methods, in this work we propose a hybrid iterative method by introducing neural network-based tools.

The rise of neural networks (NNs) brings a huge impact on searching numerical solutions of PDEs. It has been shown that NNs have powerful expressivity [4, 30], thus NNs are considered as universal approximators for the solutions of PDEs. Mainstream neural networks-based solving methods for PDEs are optimizing the parameters of such NNs to minimize the corresponding loss functions constructed by strong or variational forms of PDEs, e.g. the PINNs [14, 21, 26] , the deep Ritz method [5], and the deep Galerkin method [31]. These methods show the possibility for solving high-dimensional problems, which is indeed a challenge for traditional methods due to the curse of dimensionality. As meshless methods, they are easily implemented and friendly for coders. The drawbacks of such methods are also obvious. On the one hand, these methods are usually hard to obtain solutions as accurate as the traditional methods like FEM, on the other hand, one needs to restart a training process once some parameters of the given PDE problem (e.g. the right-hand side term of Poisson equation) are changed, which is costly to deal with parametric PDE problems. In recent years, a new research field called neural operator, is emerging to achieve the fast solving for PDEs based on neural networks. The pioneer works are the DeepONet [19, 20] and the FNO [17, 18], they directly learn the solution operators of PDEs, for example, the mapping from the right-hand side of Poisson equation to its solution. These methods present a new landscape for understanding physical reality, one may study the causality via a data-driven manner besides the PDE model. Recently numerous neural operators are proposed to achieve higher performance [7, 9, 12, 13, 25, 27, 37], as this field is developing rapidly. There are also theoretical works for the above neural operators including error analysis [15, 16, 23].

The most attractive neural operator for this work is the multiple-input operator network (MIONet) [13]. MIONet firstly faces up to the problems with multiple inputs, and generalizes the architecture and the theory of DeepONet by the tool of tensor decomposition. Such architecture also induces the tensor-based models for solving eigenvalue problems [10, 36]. An important characteristic of DeepONet/MIONet series architectures is that they encode the output functions as neural networks, rather than discrete vectors. This feature leads to efficient interpolation for the predicted solutions, and also allows the computing of derivatives of output functions during training, which prompts the work of physics-informed DeepONets [35] that learns the solution operator without any input-output pairs of data, by training a PINN-style loss. Another interesting work based on DeepONet is the HINTS [41], which in fact inspires this work. Almost all of the works of neural operators use a pure machine learning-style to solve the PDE problems, and they break away from traditional PDE numerical methods. The neural operators are developed for fast solving, but usually have limited prediction accuracy, while the traditional numerical methods can provide solutions as accurate as needed at the cost of long solving time. Surprisingly, HINTS combines the two methods, i.e. the traditional numerical solver and the neural operator. It has been shown a spectral bias/frequency principle of the training behavior of NNs in [24, 40], i.e. the NNs often fit functions from low to high frequency during the training. Based on this feature, [11] employs the deep learning solutions (e.g. PINNs) to initialize the iterative methods for PDEs. HINTS further utilizes the spectral bias, where a pre-trained DeepONet is periodically used to eliminate the low-frequency components of the error during the iteration process, thus achieve acceleration since the smoothers (e.g. Gauss-Seidel) are difficult to deal with the low-frequency errors. This idea is similar to the mulltigrid method. A significant difference is that HINTS does not require the users to generate meshes, and the cost of this method has been transferred to the developers who are responsible for the pre-trained DeepONets.

Contributions. As HINTS studies the hybrid method empirically, the lack of theoretical guidance leads to the imperfection of this method. In this work, we propose a hybrid iterative method based on MIONet. We theoretically show that MIONet is a more reasonable framework to complete the hybrid iterative method. An ideal operator regression model for the hybrid method is expected to satisfy the following three key points:

• 1.

  Supporting multiple inputs. Taking the Poisson equation

  $$\begin{cases}-\nabla\cdot(k\nabla u)=f&\quad\mbox{in}\;\Omega,\\ u=0&\quad\mbox{on}\;\partial\Omega,\end{cases}$$

  (1)

  as an example, its solution operator $u=\mathcal{G}(k,f)$ has two inputs, hence the chosen operator regression model should efficiently deal with such multiple-input case.

• 2.

  Separating linearity and nonlinearity. Since the solution operator $\mathcal{G}(k,f)$ of (1) is linear w.r.t. $f$ but nonlinear w.r.t. $k$, the operator regression model $\mathcal{M}(k,f)$ is also expected to be linear w.r.t. $f$ and nonlinear w.r.t. $k$. The linearity on $f$ will guarantee the convergence of the hybrid iterative method which will be discussed later.

• 3.

  Supporting efficient interpolation. After obtaining the predicted solution by operator regression model $\mathcal{M}$, we have to compute the values of $\mathcal{M}(k,f)(\cdot)$ at the grid points which are independent of the model $\mathcal{M}$, so that the model $\mathcal{M}$ should support efficient interpolation.

Fortunately, we find that MIONet satisfies all the three points perfectly, while other models do not (note that MIONet is a high-level framework, and one may design a targeted architecture for each branch net on a specific problem, rather than directly use FNNs.). With the proposed MIONet-based hybrid iterative method, we systematically analyze its theoretical properties, including the convergence condition, the spectral behavior, as well as the convergence rate, in terms of the errors of the discretization and the model inference. We show the theoretical results for the frequently-used smoothers, i.e. Richardson (damped Jacobi) and Gauss-Seidel. We give an upper bound of the convergence rate of the hybrid method w.r.t. the model correction period, which indicates a minimum point to make the hybrid iteration converge fastest. Several numerical examples including the Richardson (Gauss-Seidel) iteration for the 1-d (2-d) Poisson equation are shown to verify our theoretical results.

This paper is organized as follows. We first briefly introduce the operator regression model MIONet as the necessary tool in Section 2. In Section 3, we present the algorithm of the proposed hybrid iterative method in detail. We then comprehensively analyze the theoretical properties of this method in Section 4. The numerical experiments are shown in Section 5. Finally, Section 6 summarizes this work.

Let us simply introduce the operator regression model MIONet, which is a necessary tool for the later discussion of the hybrid iterative method. All the content of this section can be found in the MIONet paper [13], and the readers may refer to that paper for more details if interested. The readers familiar with MIONet may skip this section.

Let $X_{1},X_{2},...,X_{n}$ and $Y$ be $n+1$ Banach spaces, and $K_{i}\subset X_{i}$ ($i=1,...,n$) is a compact set. Now we aim to learn a continuous operator

where $v_{i}\in K_{i}$ and $u=\mathcal{G}(v_{1},...,v_{n})$. Such $\mathcal{G}$ form the space $C(K_{1}\times\cdots\times K_{n},Y)$ which is studied below. Note that $Y$ could be any Banach space. We usually consider $Y$ as a function space, such as $Y=C(K_{0}),L^{2}(K_{0}),H^{1}(K_{0})$ for a compact domain $K_{0}$.

To deal with the inputs from infinite-dimensional Banach spaces, it is necessary to firstly project them onto finite-dimensional spaces. Here we utilize the tools of Schauder basis and canonical projections.

We show two useful examples of Schauder basis as follows.

We denote the coordinate functional of $e_{i}$ by $e_{i}^{*}$, and thus

Then for a positive integer $n$, the canonical projection $P_{n}$ is defined as

We have the following property for $P_{n}$, according to which, we can represent points in an infinite-dimensional Banach space by finite coordinates within a sufficiently small projection error.

For convenience, we decompose the $P_{n}$ as

where $\varphi_{n}:X\to\mathbb{R}^{n}$ and $\psi_{n}:\mathbb{R}^{n}\to X$ are defined as

The $\varphi_{n}(x)$ are essentially the truncated coordinates for $x$. Moreover, sometimes we can further replace $\{e_{1},...,e_{n}\}$ with an equivalent basis for the decomposition of $P_{n}$, i.e.,

with a nonsingular matrix $Q\in\mathbb{R}^{n\times n}$. For example, when applying the Faber-Schauder basis (Example 1), instead of using the coordinates based on the sequence $\{e_{i}\}_{i=1}^{\infty}$, we adopt the function values evaluated at certain grid points as the coordinates, which is the same as the linear element basis in the finite element method.

By considering the canonical projection property and the injective tensor product decomposition

we can easily obtain the following approximation result.

The MIONet is immediately constructed based on this theorem, where we just replace the $g_{j}^{i}$ with neural networks. Assume that now $Y=C(K_{0})$ ($L^{2}(K_{0})$, $H^{1}(K_{0})$ etc. are also available) for a compact domain $K_{0}$.

MIONet. We construct the MIONet according to Eq. (3). Specifically, $\mathbf{g}_{i}=(g_{1}^{i},...,g_{p}^{i})^{T}\in C(\mathbb{R}^{q_{i}},\mathbb{R}^{p})$ and $\mathbf{f}=(u_{1},...,u_{p})^{T}\in C(K_{0},\mathbb{R}^{p})$ are approximated by neural networks $\tilde{\mathbf{g}}_{i}$ (called branch net $i$) and $\tilde{\mathbf{f}}$ (called trunk net). Then the network (Figure 1) is

where $\varphi_{q_{i}}^{i}$ truncates the previous $q_{i}$ coordinates w.r.t. the corresponding Schauder basis as defined in Eq. (3), $\odot$ is the Hadamard product (element-wise product), e.g.

and $\mathcal{S}$ is the summation of all the components of a vector, $b\in\mathbb{R}$ is a trainable bias.

The universal approximation theorem for MIONet can be easily obtained via Theorem 1 and the universal approximation theorems for the applied neural networks (i.e. $\tilde{\mathbf{g}}_{i}$ and $\tilde{\mathbf{f}}$) such as FNNs.

Up to now, we have already shown the architecture and properties of MIONet, which will be applied to the iterative methods later.

We present a new MIONet-based hybrid iterative method for general linear equations:

where $\mathcal{L}_{\alpha}$ is a linear differential operator depends on $\alpha$. $\alpha=\{\alpha_{1},\alpha_{2},...\}$ is a finite set consisting of several parameters which can be scalars, vectors, functions, etc.. Other boundary conditions are also available.

Without loss of generality we take the well-known Poisson equation as an example. Consider

where $\Omega\subset\mathbb{R}^{d}$ is a bounded open set. Here $f\in L^{2}(\Omega),k\in L^{\infty}(\Omega)$, and $k\geq C$ a.e. for a constant $C>0$, so that it is a classical second-order elliptic equation. The variational formula is

Now assume that $\Omega$ is a bounded polyhedral domain, $\mathcal{T}_{h}$ is a quasi-uniform and shape regular triangulation of $\Omega$. We apply the simple linear element space

where $P_{1}(K)$ denotes the function space consisting of all the linear polynomials over $K$. The finite element method is

Let $\{\phi_{i}\}_{i=1}^{n_{h}}$ be the nodal basis functions corresponding to all the interior nodes of $\mathcal{T}_{h}$. Then we obtain a discrete linear system for Eq. (10) as

where $A_{h}=(\int_{\Omega}k\nabla\phi_{i}\cdot\nabla\phi_{j}dx)_{n_{h}\times n_{h}}$, $b_{h}=(\int_{\Omega}f\phi_{i}dx)_{n_{h}\times 1}$. Denote $\mu_{h}=(\alpha_{1},...,\alpha_{n_{h}})^{T}$, then $u_{h}:=\sum_{i=1}^{n_{h}}\alpha_{i}\phi_{i}$ is the numerical solution obtained by the above system. What we next to do is solving Eq. (11) via a suitable iterative method.

For convenience, we temporarily omit the “$h$” in above notations, but we keep in mind that they indeed depend on $h$. For the discrete linear system

we consider the iterative form

where $B$ is an approximation of $A^{-1}$. The frequently-used choices of $B$ are

where $A=D+L+U$, $D$ is the diagonal of $A$, $L$ and $U$ are the strictly lower and upper triangular parts of $A$, respectively. The $\omega$ is usually considered as the relaxation parameter to be chosen. By performing a suitable iterative method, we will get a numerical solution $\mu$ for system (12).

Briefly, our hybrid iterative method is just inserting some inference steps via MIONet into the iterative process. Now we return to the operator regression model MIONet. For the problem (7), what we are concerned with is in fact the solution operator

Now assume that we have a dataset as $\mathcal{T}=\{(k_{i},f_{i}),u_{i}\}_{i=1}^{N}$ satisfying $\mathcal{G}(k_{i},f_{i})=u_{i}$. We pre-train a MIONet, denoted by $\mathcal{M}$, on this dataset with loss function

where $Y$ could be $C(\overline{\Omega})$, $L^{2}(\Omega)$, $H^{1}(\Omega)$, etc.. Note that here we set $\mathcal{M}$ linear with respect to $f$ due to Theorem 1. An illustration is presented in Figure 2.

Denote $r^{(m)}:=b-A\mu^{(m)}$, which is the residual vector of $m$-th step. We expect to find a residual function related to the residual vector $r^{(m)}$. Define the residual function $\bar{r}^{(m)}:=\sum_{i=1}^{n}\beta_{i}\phi_{i}$, where $\beta=(\beta_{1},...,\beta_{n})^{T}$ satisfies

Specially, we let $\bar{r}^{(0)}=f$. Note that $\bar{r}^{(m)}$ is linear with respect to $r^{(m)}$, we denote their relationship by a linear operator $\mathcal{H}$ as

We further define the interpolation operator as

where $x_{i}\in\Omega$ is the corresponding nodal coordinate of $\phi_{i}$, i.e. $\phi_{i}(x_{j})=\delta_{i}^{j}$ (Dirac delta function).

Then the hybrid iterative method can be written as

where $M$ is a positive integer to be chosen, we regard it as the correction period. The algorithm is summarized in Algorithm 1.

Approximation and padding for residual function. There is a new linear system (19) involved. We may choose an approximation for $\beta$, for example:

With the solved $\beta$, we obtain the $\bar{r}^{(m)}=\mathcal{H}r^{(m)}$ expressed as a continuous piece-wise linear function w.r.t. $\{\phi_{i}\}$, which does not contain the boundary nodes. Hence the current $\bar{r}^{(m)}$ equals to 0 on $\partial\Omega$, which is regarded as a “zero padding”. Here we may use other better padding for $\bar{r}^{(m)}$, for example, add boundary nodes whose values equal to the values of their corresponding nearest interior nodes (we can also apply other padding strategy like reflection padding, etc.). For convenience, we still use the notation $\bar{r}^{(m)}$ and $\mathcal{H}$, i.e.

Note that the above approximation brings an acceptable error for this hybrid method. Since $\bar{r}^{(0)}=f$, $\mu^{(1)}$ is exactly the prediction by MIONet. We may also think of $\mu^{(1)}$ as an initial guess provided by MIONet for the iterative method. When feeding $\bar{r}^{(m)}$ to $\mathcal{M}$, we need to evaluate the linear element function $\bar{r}^{(m)}$ at some sampling points (usually with smaller scale compared to the number of nodes $n$) for $\mathcal{M}$.

Efficient interpolation. The interpolation $\mathbb{I}$ is easy to compute, since the prediction

is expressed by a neural network (encoded as a trunk net) rather than a discrete vector, and we only need to directly feed the interpolated points into the trunk net of the ONet. This is indeed an important advantage of the ONets series architectures.

Hybrid iteration with multigrid method. Since there is no requirement for the basic iterative method in this algorithm, we can also utilize other advanced iterative methods such as the multigrid method here. We will show more details in the experiment part.

Inhomogeneous boundary condition. Here we further investigate the inhomogeneous boundary condition. Consider

where $T$ is the trace operator. In such a case, we aim to learn an operator mapping from $(k,f,g)$ to the solution $u$, i.e.,

Of course we can train an MIONet $\mathcal{M}(k,f,g)$ to learn this operator. However, this strategy losses the linearity w.r.t. the residual vector/function during the iteration, and we strive to preserve this linearity to guarantee the convergency of the hybrid method. Denote by $E$ the trace extension operator for (27), then $T(Eg)=g$ for $g\in H^{1/2}(\partial\Omega)$. We derive that

where $u=u_{0}+Eg$. Therefore the model can be modified to

where $\mathcal{M}$ is a standard MIONet which is linear w.r.t. the second input as before. Although the trace extension operator $E$ is unique defined, we may choose a suitable constructed extension as needed, for example the continuous piece-wise linear function related to a grid, which takes 0 at the interior nodes and matches $g$ at the boundary nodes. Note that here we regard $f+\nabla\cdot(k\nabla Eg)$ as an element in $H^{-1}(\Omega)$.

There is in fact another way which is much easier to achieve. It is not difficult to find that $\mathcal{G}(k,f,g)$ is linear w.r.t. $(f,g)\in L^{2}(\Omega)\times H^{1/2}(\partial\Omega)$, i.e.,

then we naturally apply the model

where $(f,g)$ is a simple cartesian product of $f$ and $g$, and it will be fed into the linear branch net of $\mathcal{M}$. In such a case, the discrete algebraic system of (27) is augmented with several rows and columns taking into account the boundary nodes, where the boundary block is the identity. The unknown and residual vectors contain not only the coefficients of the interior nodes but also the boundary nodes. We apply this strategy in our experiment.

Necessary convergence condition for operator regression models. There is a necessary condition for the convergence of the hybrid iterative method. Assume that the hybrid method converges with a model $\mathcal{M}$, then

and we have

As the interpolation operator $\mathbb{I}$ is related to the mesh $\mathcal{T}_{h}$, and the above equation holds for arbitrary $\mathbb{I}_{\mathcal{T}_{h}}$, we know

Obviously, Eq. (35) holds as long as $\mathcal{M}$ is linear w.r.t. the second input.

Consider the one-dimensional two-point boundary value problem of system (7) where $k=1$, $\Omega=(0,1)$, i.e.,

We choose a uniform grid with size $h=\frac{1}{n+1}$ and apply the linear element bases functions $\{\phi_{i}\}_{i=1}^{n}$ whose nodal points are interior, then the above system leads to a discrete system as

where

and $b_{i}=\int_{0}^{1}f\phi_{i}dx$. We can easily write out the eigenvalues and the eigenvectors of $A$ as

Denote $\Xi=(\xi^{1},...,\xi^{n})$, then $\Xi$ is symmetric and orthogonal, i.e. $\Xi^{T}=\Xi$ and $\Xi^{-1}=\Xi$.