Solving multiscale elliptic problems by sparse radial basis function
  neural networks
Machine learning has been successfully applied to various fields of
scientific computing in recent years. In this work, we propose a sparse radial
basis function neural network method to solve elliptic partial differential
equations (PDEs) with multiscale coefficients. Inspired by the deep mixed
residual method, we rewrite the second-order problem into a first-order system
and employ multiple radial basis function neural networks (RBFNNs) to
approximate unknown functions in the system. To aviod the overfitting due to
the simplicity of RBFNN, an additional regularization is introduced in the loss
function. Thus the loss function contains two parts: the $L_2$ loss for the
residual of the first-order system and boundary conditions, and the $\ell_1$
regularization term for the weights of radial basis functions (RBFs). An
algorithm for optimizing the specific loss function is introduced to accelerate
the training process. The accuracy and effectiveness of the proposed method are
demonstrated through a collection of multiscale problems with scale separation,
discontinuity and multiple scales from one to three dimensions. Notably, the
$\ell_1$ regularization can achieve the goal of representing the solution by
fewer RBFs. As a consequence, the total number of RBFs scales like
$\mathcal{O}(\varepsilon^{-n\tau})$, where $\varepsilon$ is the smallest scale,
$n$ is the dimensionality, and $\tau$ is typically smaller than $1$. It is
worth mentioning that the proposed method not only has the numerical
convergence and thus provides a reliable numerical solution in three dimensions
when a classical method is typically not affordable, but also outperforms most
other available machine learning methods in terms of accuracy and robustness.
Zhiwen Wang, Minxin Chen, Jingrun Chen
# Solving multiscale elliptic problems by sparse radial basis function neural networks

###### Abstract

Machine learning has been successfully applied to various fields of scientific computing in recent years. For problems with multiscale features such as flows in porous media and mechanical properties of composite materials, however, it remains difficult for machine learning methods to resolve such multiscale features, especially when the small scale information is present. In this work, we propose a sparse radial basis function neural network method to solve elliptic partial differential equations (PDEs) with multiscale coefficients. Inspired by the deep mixed residual method, we rewrite the second-order problem into a first-order system and employ multiple radial basis function neural networks (RBFNNs) to approximate unknown functions in the system. To avoid the overfitting due to the simplicity of RBFNN, an additional regularization is introduced in the loss function. Thus the loss function contains two parts: the \(L_{2}\) loss for the residual of the first-order system and boundary conditions, and the \(\ell_{1}\) regularization term for the weights of radial basis functions (RBFs). An algorithm for optimizing the specific loss function is introduced to accelerate the training process. The accuracy and effectiveness of the proposed method are demonstrated through a collection of multiscale problems with scale separation, discontinuity and multiple scales from one to three dimensions. Notably, the \(\ell_{1}\) regularization can achieve the goal of representing the solution by fewer RBFs. As a consequence, the total number of RBFs scales like \(\mathcal{O}(\varepsilon^{-nT})\), where \(\varepsilon\) is the smallest scale, \(n\) is the dimensionality, and \(\tau\) is typically smaller than 1. It is worth mentioning that the proposed method not only has the numerical convergence and thus provides a reliable numerical solution in three dimensions when a classical method is typically not affordable, but also outperforms most other available machine learning methods in terms of accuracy and robustness.

Keywords:Multiscale elliptic problem, Radial basis function neural network, multi-scale, \(\ell_{1}\) regularization. +
Footnote †: journal:

0224007006@stu.suda.edu.cn

chenminxin@suda.edu.cn

20224007006@stu.suda.edu.cn

(Z. Wang), chenminxin@suda.edu.cn

(M. Chen), jingrunchen@ustc.edu.cn
Introduction

Multiscale phenomena are common in nature and mathematical modeling typically involves a small parameter \(\varepsilon\) which characterizes the ratio between the smallest characteristic length over the size of interest [6]. In some scenarios, the mechanical response of composite materials and flow property of porous media can be modeled by the elliptic problems with multiscale coefficients

\[\begin{cases}-\operatorname{div}(a^{\varepsilon}(\mathbf{x})\nabla u^{ \varepsilon}(\mathbf{x}))=f(\mathbf{x})&\mathbf{x}\in\Omega,\\ u^{\varepsilon}(\mathbf{x})=g(\mathbf{x})&\mathbf{x}\in\partial\Omega,\end{cases} \tag{1.1}\]

where \(\Omega\subset\mathbb{R}^{n}\),\(n=1,2,3\) is a bounded domain, \(\varepsilon\ll 1\), and \(f:\Omega\to\mathbb{R}\) is the source term.

For small \(\varepsilon\), classical numerical methods cannot solve (1.1) in an effective manner, where the number of unknowns scales like \(\mathcal{O}(\varepsilon^{-n\tau})\) and \(\tau>1\). Meanwhile, by the homogenization theory [3, 26], as \(\varepsilon\to 0\), an effective model can be obtained from (1.1) based on the periodicity assumption of \(a^{\varepsilon}(\mathbf{x})\). Over the past a couple of decades, a large number of multiscale methods have been developed for multiscale elliptic problems in the literature, such as the multiscale finite element method (MsFEM) [9, 11], wavelet homogenization techniques [5], the heterogeneous multiscale method (HMM) [7, 8], the metric-based upscaling [19], and the local orthogonal decomposition (LOD) [16]. Typically, these methods solve multiple local problems in advance to extract the macroscopic information and then solve one global problem. The total computational cost is at least \(\mathcal{O}(\varepsilon^{-n})\).

In recent years, a series of deep neural network methods for solving PDEs are derived. Deep Ritz method (DRM) [27] employs the variational formulation as the loss function. Deep Galerkin method (DGM) [24] and Physics-informed neural networks (PINNs) [22] utilize the residual of the PDE as the loss function. Weak Adversarial Network (WAN) [28] employs the weak formulation to construct the loss function. Multilayer perception (MLP) and residual neural network (ResNet) are commonly used to approximate the solution. For multiscale problems, DNN-based methods are also developed. Multiscale deep neural network (MscaleDNN) [14] approximates the solution by converting the original data to a low frequency space and [13] improves the MscaleDNN algorithm by a smooth and localized activation function.

Radial basis function neural network (RBFNN) represents an attractive alternative to other neural network models [18]. One reason is that it forms a unified connection between function approximation, regularization, noisy interpolation, classification, and density estimation. It is also the case that training a RBFNN is faster than training MLP networks [18]. RBFNN is a simple three-layer feedforward neural network. The first layer is the input layer, representing the input features. The second layer is the hidden layer, composed of radial basis activation functions. The third layer is the output layer. Fig. 1 shows the network architecture of RBFNN in four dimensions with five neurons. Like most neural networks, RBFNN also has the universal approximation property [20]. Meanwhile, compared with deep neural networks, RBFNN avoids the tedious
and lengthy calculation of back propagation between the input layer and the output layer, significantly accelerating the training process of the network.

Often, optimizing the empirical risk leads to the overfitting problem owing to the simplicity of RBFNN, and causes a poor generalization capability. To tackle this issue, \(\ell_{1}\) regularization techniques [21, 23], become common components in many active fields, such as sparse representations of images [2], molecular surface [10] and point cloud surface reconstruction [25]. The basic idea is to add a regularization term to the loss function to penalize over-complicated solutions. For \(\ell_{1}\) regularization, a weighted \(\ell_{1}\) norm of the parameter vector is added to the loss function, which penalizes the sum of the absolute values of the parameters.

In this work, we propose a sparse radial basis function neural network (SRBFNN) method for solving multiscale elliptic equations, where RBFNN is used to approximate the multiscale solution and the \(\ell_{1}\) regularization term is added to the loss function to guarantee the sparsity of RBF representation. In the modeling stage, by rewriting the second-order equation into a first-order system as shown in the deep mixed residual method [15], we obtain an augmented problem in the sense that both the solution and the auxiliary variables related to first-order derivatives of the solution are unknown functions to be approximated. An algorithm for optimizing the specific loss function is introduced to accelerate the training process. The accuracy and effectiveness of the proposed method are demonstrated through a collection of multiscale problems with scale separation, discontinuity and multiple scales from one to three dimensions. Notably, the \(\ell_{1}\) regularization can achieve the goal of representing the solution by fewer RBFs. As a consequence, the total number of RBFs scales like \(\mathcal{O}(\varepsilon^{-n\tau})\), where \(\varepsilon\) is the smallest scale, \(n\) is the dimensionality, and \(\tau\) is typically smaller than 1.

The paper is organized as follows. Section 2 introduces the SRBFNN. An algorithm is proposed to optimize the loss function with the \(\ell_{1}\) regularization during the training pro

Figure 1: The network architecture of radial basis function neural network.
cess in Section 3. Section 4 provides numerical experiments from one to three dimensions to demonstrate the accuracy and effectiveness of the proposed method. Conclusions and discussions are drawn in Section 5.

## 2 Sparse radial basis function neural network

### Radial basis function

The function \(\Phi:\mathbb{R}^{n}\rightarrow\mathbb{R}\) is called a radial basis function, which can be passed through a one-dimensional (1D) function \(\phi:[0,\infty)\rightarrow\mathbb{R}\), this means

\[\Phi(\mathbf{x})=\phi(\|\mathbf{x}\|)=\phi(r),\quad r:=\|\mathbf{x}\|,\]

where \(\|\cdot\|\) represents the \(L_{2}\) norm in \(\mathbb{R}^{n}\).

Some common types of RBFs:

* Gaussian function \[\Phi(r)=e^{-d^{2}r^{2}};\]
* Multiquadric function (MQ) \[\Phi(r)=\sqrt{1+d^{2}r^{2}};\]
* Inverse multiquadric function (IMQ) \[\Phi(r)=1/\sqrt{1+d^{2}r^{2}},\]

where \(d\) is the shape parameter that plays an important role in the approximation accuracy. As shown in [4], a small shape parameter results in an ill-conditioned system, while the corresponding RBF method exhibits the high approximation accuracy. In contrast, a large shape parameter results in an well-conditioned system, but the approximation accuracy is low.

In this work, we use the ellipse Gaussian RBFs

\[\Phi(\mathbf{x};\mathbf{c},\mathbf{D})=\mathrm{e}^{-\|\mathbf{D}(\mathbf{x}- \mathbf{c})\|}, \tag{2.1}\]

where \(\mathbf{D}=\mathrm{diag}(d_{1},...,d_{n})\) is the shape parameter and \(\mathbf{c}=(c_{1},...,c_{n})^{T}\) is the center. The ellipse Gaussian RBF has the good approximation property in a specific neighborhood of the center with different radius.
### Sparse radial basis function neural network

The RBFNN can be expressed as

\[\mathbf{NN}(\mathbf{x};\mathbf{w},\mathbf{c},\mathbf{D})=\sum_{i=1}^{N}w_{i}\Phi_{ i}(\mathbf{x};\mathbf{c}_{i},\mathbf{D}_{i}), \tag{2.2}\]

where \(\Phi_{i}(\mathbf{x};\mathbf{c}_{i},\mathbf{D}_{i}),i=1,...,N\) are the ellipse Gaussian RBFs and \(\mathbf{w}=(w_{1},..,w_{N})^{T}\) is the corresponding weight. For the ease of exposition, we denote the trainable parameters \(\mathbf{w},\mathbf{c},\mathbf{D}\) by \(\boldsymbol{\theta}\), i.e.,

\[\boldsymbol{\theta}=\{\mathbf{w}\in\mathbb{R}^{N},\mathbf{c}\in\mathbb{R}^{N \times n},\mathbf{D}\in\mathbb{R}^{N\times n}|N\text{ is the number of basis, }n\text{ is the dimension of the input }\mathbf{x}\}.\]

The RBFNN whose loss function contains the \(\ell_{1}\) regularization term for the weight vector \(\mathbf{w}\) is called the sparse RBFNN (SRBFNN).

### Loss function

We consider multiscale elliptic equations of the form (1.1). The total loss function for SRBFNN is defined as

\[Loss(\mathbf{x};\boldsymbol{\theta})=L_{s}(\mathbf{x};\boldsymbol{\theta})+ \lambda_{3}\|\mathbf{w}^{u}\|_{1}, \tag{2.3}\]

where \(L_{s}\) is the \(L_{2}\) loss function for the first-order system, \(\lambda_{3}\) is the penalty parameter, and \(\mathbf{w}^{u}\) is the weight of the RBFNN that approximates the solution. \(\|\cdot\|_{1}\) denotes the \(\ell_{1}\) norm. Explicit forms of \(L_{s}\) from one to three dimensions will be specified in (2.6), (2.9), and (2.12).

#### 2.3.1 \(L_{2}\) loss function in 1D

Inspired by the MIM method [15], we rewrite (1.1) in 1D, using an auxiliary variable \(p\), into a first-order system:

\[\begin{cases}p-a^{\varepsilon}u_{x}=0&x\in\Omega,\\ p_{x}+f=0&x\in\Omega,\\ u=g&x\in\partial\Omega.\end{cases} \tag{2.4}\]

Then we use two RBFNNs (\(\mathbf{NN}_{1}\), \(\mathbf{NN}_{2}\)) to approximate \(u\) and \(p\) as

\[\widetilde{u}=\mathbf{NN}_{1}(x;\boldsymbol{\theta}^{u}),\quad\widetilde{p}= \mathbf{NN}_{2}(x;\boldsymbol{\theta}^{p}), \tag{2.5}\]

where \(\boldsymbol{\theta}^{u}\) and \(\boldsymbol{\theta}^{p}\) are the neural network parameters used to approximate \(u\) and \(p\), respectively. Then the \(L_{2}\) loss function for the system (2.4) is defined by

\[L_{s}(\mathbf{x};\boldsymbol{\theta}^{u},\boldsymbol{\theta}^{p})=\|a^{ \varepsilon}\widetilde{u}_{x}-\widetilde{p}\|_{2,\Omega}^{2}+\lambda_{1}\| \widetilde{p}_{x}+f\|_{2,\Omega}^{2}+\lambda_{2}\|\widetilde{u}-g\|_{2,\partial \Omega}^{2}, \tag{2.6}\]

where \(\lambda_{1}\) and \(\lambda_{2}\) are penalty parameters given in advance. The three terms in (2.6) measure how well the approximate solution and the auxiliary variable satisfy the first-order system (2.4) and the boundary condition, respectively. Due to the simple form of the RBFNN defined in (2.2), the corresponding derivatives \(\widetilde{u}_{x}\), \(\widetilde{p}_{x}\) can be easily calculated.
#### 2.3.2 \(L_{2}\) loss function in 2D

Analogously, we rewrite (1.1) in 2D, using auxiliary variables \(p,q\), as

\[\left\{\begin{aligned} p-a^{\varepsilon}\cdot\frac{\partial u}{ \partial x}&=0&\mathbf{x}\in\Omega,\\ q-a^{\varepsilon}\cdot\frac{\partial u}{\partial y}& =0&\mathbf{x}\in\Omega,\\ \frac{\partial p}{\partial x}+\frac{\partial q}{\partial y}+f& =0&\mathbf{x}\in\Omega,\\ u&=g&\mathbf{x}\in\partial\Omega.\end{aligned}\right. \tag{2.7}\]

Then we use three RBFNNs (\(\mathbf{NN}_{1}\), \(\mathbf{NN}_{2}\), \(\mathbf{NN}_{3}\)) to approximate \(u,p,q\) as

\[\widetilde{u}=\mathbf{NN}_{1}(\mathbf{x};\boldsymbol{\theta}^{u}),\quad \widetilde{p}=\mathbf{NN}_{2}(\mathbf{x};\boldsymbol{\theta}^{p}),\quad \widetilde{q}=\mathbf{NN}_{3}(\mathbf{x};\boldsymbol{\theta}^{q}). \tag{2.8}\]

The \(L_{2}\) loss function for the system (2.7) is defined as

\[L_{s}(\mathbf{x};\boldsymbol{\theta}^{u},\boldsymbol{\theta}^{p},\boldsymbol{ \theta}^{q})=\|a^{\varepsilon}\cdot\frac{\partial\widetilde{u}}{\partial x} -\widetilde{p}\|_{2,\Omega}^{2}+\|a^{\varepsilon}\cdot\frac{\partial \widetilde{u}}{\partial y}-\widetilde{q}\|_{2,\Omega}^{2}+\lambda_{1}\|\frac{ \partial\widetilde{p}}{\partial x}+\frac{\partial\widetilde{q}}{\partial y}+ f\|_{2,\Omega}^{2}+\lambda_{2}\|\widetilde{u}-g\|_{2,\partial\Omega}^{2}. \tag{2.9}\]

#### 2.3.3 \(L_{2}\) loss function in 3D

We rewrite (1.1) in 3D, using auxiliary variables \(p,q,r\), as

\[\left\{\begin{aligned} p-a^{\varepsilon}\cdot\frac{ \partial u}{\partial x}&=0&\mathbf{x}\in\Omega,\\ q-a^{\varepsilon}\cdot\frac{\partial u}{\partial y}& =0&\mathbf{x}\in\Omega,\\ r-a^{\varepsilon}\cdot\frac{\partial u}{\partial z}& =0&\mathbf{x}\in\Omega,\\ \frac{\partial p}{\partial x}+\frac{\partial q}{\partial y}+ \frac{\partial r}{\partial z}+f&=0&\mathbf{x}\in \Omega,\\ u&=g&\mathbf{x}\in\partial\Omega.\end{aligned}\right. \tag{2.10}\]

Then we use four RBFNNs (\(\mathbf{NN}_{1}\), \(\mathbf{NN}_{2}\), \(\mathbf{NN}_{3}\), \(\mathbf{NN}_{4}\)) to approximate \(u,p,q,r\) as

\[\widetilde{u}=\mathbf{NN}_{1}(\mathbf{x};\boldsymbol{\theta}^{u}),\quad \widetilde{p}=\mathbf{NN}_{2}(\mathbf{x};\boldsymbol{\theta}^{p}),\quad \widetilde{q}=\mathbf{NN}_{3}(\mathbf{x};\boldsymbol{\theta}^{q}),\quad \widetilde{r}=\mathbf{NN}_{4}(\mathbf{x};\boldsymbol{\theta}^{r}). \tag{2.11}\]

The \(L_{2}\) loss function for the system (2.10) is

\[L_{s}(\mathbf{x};\boldsymbol{\theta}^{u},\boldsymbol{\theta}^{p}, \boldsymbol{\theta}^{q},\boldsymbol{\theta}^{r})= \|a^{\varepsilon}\cdot\frac{\partial\widetilde{u}}{\partial x}- \widetilde{p}\|_{2,\Omega}^{2}+\|a^{\varepsilon}\cdot\frac{\partial \widetilde{u}}{\partial y}-\widetilde{q}\|_{2,\Omega}^{2}+\|a^{\varepsilon} \cdot\frac{\partial\widetilde{u}}{\partial z}-\widetilde{r}\|_{2,\Omega}^{2} \alpha^{+} \tag{2.12}\] \[\lambda_{1}\|\frac{\partial\widetilde{p}}{\partial x}+\frac{ \partial\widetilde{q}}{\partial y}+\frac{\partial\widetilde{r}}{\partial z}+ f\|_{2,\Omega}^{2}+\lambda_{2}\|\widetilde{u}-g\|_{2,\partial\Omega}^{2}.\]
### Network initialization

All RBFNNs use the same initialization strategy for the parameters \(\mathbf{w}\), \(\mathbf{c}\), \(\mathbf{D}\). The center \(\mathbf{c}\) is initialized with a uniform distribution over \(\Omega\). For the weight \(\mathbf{w}\), we initialize them with a uniform distribution U(0,1). The initialization of the shape parameter \(\mathbf{D}\) is extremely important for the network convergence in multiscale problems. According to the homogenization theory [3], the solution \(u^{\varepsilon}\) to (1.1) admits the following asymptotic expansion:

\[u^{\varepsilon}(\mathbf{x})=u_{0}(\mathbf{x})+\varepsilon u_{1}(\mathbf{x}, \frac{\mathbf{x}}{\varepsilon})+\varepsilon^{2}u_{2}(\mathbf{x},\frac{ \mathbf{x}}{\varepsilon})+\cdots, \tag{2.13}\]

where \(u_{0}\) is the homogenized solution, and \(u_{1}\), \(u_{2}\),... are high-order terms and depend on the derivatives of \(u_{0}\). Therefore, to resolve the smallest scale of the solution, we initialize \(\mathbf{D}\) with uniform distribution U(0, \(\frac{1}{\varepsilon}\)).

## 3 Training process of SRBFNN

In the training process of SRBFNN, we aim to find a sparse weight vector, \(\mathbf{w}^{u}\), while \(\mathbf{NN}_{i}\),\(i\)=1,...,\(n\)+1 are also good fits of the solution and the auxiliary variables. Therefore, an algorithm to optimize the loss function with \(\ell_{1}\) regularization during the training phase is proposed in Algorithm 1.
```
Input : The coefficient \(a^{e}\), the source term \(f\), and the boundary condition \(g\) in (1.1). Output : The list of parameters for \(\mathbf{NN}_{i},i=1,...,n+1\). Step 1. initialize \(\mathbf{NN}_{i}\), \(i=1,...,n+1\), \(n\) is the dimension; Step 2. generate two sets of points uniformly distributed over \(\Omega\) and \(\partial\Omega\): \(\left\{\mathbf{x}_{m}\right\}_{m=1}^{M_{r}}\) in \(\Omega\) and \(\left\{\widehat{\mathbf{x}}_{m}\right\}_{m=1}^{M_{b}}\) on \(\partial\Omega\), \(M_{r}\) and \(M_{b}\) are the number of sample points; Step 3. initialize some hyperparameters; Step 3.1. set the penalty coefficients of the loss functions, i.e., \(\lambda_{1}\), \(\lambda_{2}\) and \(\lambda_{3}\); Step 3.2. set the maximum iteration \(MaxNiter\) and the number of sparse optimization iterations \(SparseNiter\). Set \(Check_{iter}\) for deleting the neglectable weight every \(Check_{iter}\) iterations; Step 3.3. initialize the threshold, \(thres=0.0\), for adding \(\ell_{1}\) regularization and the recorded loss, \(L^{rec}=0.0\), for recording latest \(L_{2}\) loss per \(Check_{iter}\) iterations. Initialize the learning rate; Step 3.4. set the tolerance, \(tol_{1}\), to determine whether \(L_{2}\) loss is convergent and the tolerance, \(tol_{2}\), for neglectable weight in \(\mathbf{NN}_{1}\). Set \(batchsize\), \(Niter=0\) and \(j=0\); Step 4. start following training process; while\(Niter\leq MaxNiter\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do Step 4.1. calculate the \(L_{2}\) loss function, \(L_{s}\) (in 2.6, 2.9, 2.12), with \(batchsize\) sample points; Step 4.2. if\(L_{s}>thres\)or\(Niter>SparseNiter\)then \(\lambda_{3}=0\); else \(\lambda_{3}>0\);  end if Step 4.3. update the network parameters by ADAM optimizer;  end for if\(Niter\%Check_{iter}==0\)then Step 4.4. if\(|L^{rec}-L_{s}|<tol_{1}\)and\(j==0\)then  set \(thres=L_{s}+tol_{1}\); \(j+=1\);  end if Step 4.5. if\(j>0\)and\(Niter<SparseNiter\)then  delete neglectable basis function, whose coefficient is close to zero, i.e., \(|w_{i}|<tol_{2}\);  end if Step 4.6. record the \(L_{2}\) loss as the \(L^{rec}\) every \(Check_{iter}\) iterations, i.e., \(L^{rec}=L_{s}\) ;  end for  end if
```

**Algorithm 1**Training process of the SRBFNN

``` Input : The coefficient \(a^{e}\), the source term \(f\), and the boundary condition \(g\) in (1.1). Output : The list of parameters for \(\mathbf{NN}_{i},i=1,...,n+1\). Step 1. initialize \(\mathbf{NN}_{i}\), \(i=1,...,n+1\), \(n\) is the dimension; Step 2. generate two sets of points uniformly distributed over \(\Omega\) and \(\partial\Omega\): \(\left\{\mathbf{x}_{m}\right\}_{m=1}^{M_{r}}\) in \(\Omega\) and \(\left\{\widehat{\mathbf{x}}_{m}\right\}_{m=1}^{M_{b}}\) on \(\partial\Omega\), \(M_{r}\) and \(M_{b}\) are the number of sample points; Step 3. initialize some hyperparameters; Step 3.1. set the penalty coefficients of the loss functions, i.e., \(\lambda_{1}\), \(\lambda_{2}\) and \(\lambda_{3}\); Step 3.2. set the maximum iteration \(MaxNiter\) and the number of sparse optimization iterations \(SparseNiter\). Set \(Check_{iter}\) for deleting the neglectable weight every \(Check_{iter}\) iterations; Step 3.3. initialize the threshold, \(thres=0.0\), for adding \(\ell_{1}\) regularization and the recorded loss, \(L^{rec}=0.0\), for recording latest \(L_{2}\) loss per \(Check_{iter}\) iterations. Initialize the learning rate; Step 3.4. set the tolerance, \(tol_{1}\), to determine whether \(L_{2}\) loss is convergent and the tolerance, \(tol_{2}\), for neglectable weight in \(\mathbf{NN}_{1}\). Set \(batchsize\), \(Niter=0\) and \(j=0\); Step 4. start following training process; while\(Niter\leq MaxNiter\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+=batchsize\)do \(Niter+=1\); for\(i=0,i\leq M_{r},i+batchsize\)do \(Niter+=1\
Step 1 initializes the parameters of SRBFNN. Step 2 generates two sets of training points \(\{\mathbf{x}_{m}\}_{m=1}^{M_{r}}\) and \(\{\widetilde{\mathbf{x}}_{m}\}_{m=1}^{M_{b}}\) uniformly distributed over \(\Omega\) and \(\partial\Omega\), respectively. Step 3.1 sets the penalty coefficients \(\lambda_{1}\), \(\lambda_{2}\), \(\lambda_{3}\). Step 3.2 sets the number of maximum iterations and the number of sparse optimization iterations, i.e., the max number of iterations with the loss function having the \(\ell_{1}\) regularization term. Step 3.3 sets a threshold to the \(L_{2}\) loss, which denotes the sign for adding the \(\ell_{1}\) regularization, i.e., setting \(\lambda_{3}\!>\!0\). Step 3.4 sets a tolerance for the stability condition of \(L_{2}\) loss. Step 4.1 calculates the \(L_{2}\) loss with \(batchsize\) sample points. When the \(L_{s}\) is lower than \(thres\) and \(Niter\!<\!SparseNiter\), we add the \(\ell_{1}\) regularization in step 4.2. Step 4.3 selects the ADAM optimizer [12] to update the parameters and its pipeline is as follows

\[m_{k} =\!\beta_{1}\!\cdot\!m_{k-1}\!+\!(1\!-\!\beta_{1})\!\cdot\! \nabla Loss_{k},\] \[v_{k} =\!\beta_{2}\!\cdot\!v_{k-1}\!+\!(1\!-\!\beta_{2})\!\cdot\! \nabla(Loss_{k})^{2},\] \[\boldsymbol{\theta}_{k} =\!\boldsymbol{\theta}_{k-1}\!-\!lr\!\times\!\frac{\frac{m_{k}}{ 1-\beta_{1}^{k}}}{\sqrt{\frac{v_{k}}{1-\beta_{2}^{k}}}\!+\!\epsilon},\]

where \(lr\) is the learning rate. \(\beta_{1}\), \(\beta_{2}\) and \(\epsilon\) are set as default values (\(\beta_{1}\) = 0.9, \(\beta_{2}\) = 0.999 and \(\epsilon\) = \(10^{-8}\)), \(m_{k}\) is the \(k\)th biased first moment estimate (\(m_{0}\) = 0). \(v_{k}\) is the \(k\)th biased second raw estimate (\(v_{0}\) = 0). Step 4.4 determines whether the difference between the \(L_{2}\) loss at current iteration and the recorded loss is less than the tolerance, \(tol_{1}\). If it does, \(L_{2}\) loss is convergent and we let \(thres\) equal to the \(L_{2}\) loss of the current iteration plus \(tol_{1}\). Step 4.5 determines whether \(L_{2}\) loss is convergent and \(Niter\!<\!SparseNiter\). If it does, we delete the neglectable basis function in \(\boldsymbol{\theta}^{\boldsymbol{u}}\), whose coefficient is close to 0. Step 4.6 records the value of \(L_{2}\) loss per \(Check_{iter}\) iterations.

## 4 Numerical experiments

In this section, we present several numerical examples to illustrate the accuracy and effectiveness of our method. The coefficients are periodic functions in Example 1 and Example 5, two-scale problems with scale separation in Example 2, Example 3, and Example 6, a discontinuous periodic function in Example 4, a multiscale case in Example 7, respectively. Example 8 is a three dimensional example with the periodic coefficient.

We shall use the following three relative errors to measure the approximation accuracy.

\[err_{2}\!=\!\frac{\|u^{S}\!-\!u^{F}\|_{L_{2}(\Omega)}}{\|u^{F}\|_{L_{2}( \Omega)}},\quad err_{\infty}\!=\!\frac{\|u^{S}\!-\!u^{F}\|_{L_{\infty}(\Omega )}}{\|u^{F}\|_{L_{\infty}(\Omega)}},\quad err_{H_{1}}\!=\!\frac{\|u^{S}\!-\!u^ {F}\|_{H_{1}(\Omega)}}{\|u^{F}\|_{H_{1}(\Omega)}}, \tag{4.1}\]

where \(u^{S}\) is the approximate solution of SRBFNN and \(u^{F}\) is the reference solution of finite difference method (FDM) with a very fine mesh size. The reference solution is not
available in 3D since FDM is too expensive to solve a multiscale problem. \(\|\cdot\|_{L_{2}(\Omega)}\) is the \(L_{2}\) norm, \(\|\cdot\|_{L_{\infty}(\Omega)}\) is the \(L_{\infty}\) norm and \(\|\cdot\|_{H_{1}(\Omega)}\) is defined as

\[\|u^{S}\|_{H_{1}(\Omega)}=\sqrt{\|u^{S}\|^{2}_{L_{2}(\Omega)}+\| \nabla u^{S}\|^{2}_{L_{2}(\Omega)}}. \tag{4.2}\]

where \(\nabla u^{S}\) is calculated by the auxiliary variables, i.e., \(p/a^{\varepsilon}\) and \(q/a^{\varepsilon}\).

### Numerical experiments in 1D

**Setting:** In all one-dimensional cases, unless otherwise specified, the initial learning rate is set to be 0.1 and reduced by 1/10 every 300 iterations until it is less than \(10^{-5}\). We sample 10000 training points in the domain with a uniform distribution. The mesh size of FDM is set to 0.0001. \(batchsize\) is 2048, \(MaxNiter\) is 3000, \(SparseNiter\) is 2000, \(tol_{1}\) is \(10^{-3}\), \(tol_{2}\) is \(10^{-5}\), \(Check_{iter}\) is 100. The hyperparameters \(\lambda_{1}\), \(\lambda_{2}\) and \(\lambda_{3}\) are initialized to 1.0, 100.0 and 0.001, respectively.

We use six different values, \(\varepsilon=0.5\),0.1,0.05,0.01,0.005,0.002, to check the accuracy of our method. The initial number of RBFs is set to 100, 200, 300, 500, 1000, 1500, respectively. The computational domain is the unit interval. \(f(x)\) and \(g(x)\) in (1.1) are set to 1.0 in all 1D cases. All experiments are conducted on the GPU device RTX 2080Ti.

Example 1Consider the case that \(a^{\varepsilon}(x)\) is a periodic function

\[a^{\varepsilon}(x)=2+\sin(\frac{2\pi x}{\varepsilon}). \tag{4.3}\]

Fig. 2 plots the solution and its derivative obtained by SRBFNN and FDM when \(\varepsilon=0.01\), 0.002, respectively. One can see that SRBFNN has a good approximation for both the solution and the derivative.
Table 1 shows three relative errors and the number of basis functions in the final solution from SRBFNN after the training process.

From Table 1, one general observation is that the number of RBFs needed in the final solution increases as the \(\varepsilon\) decreases. In addition, the relative \(L_{2}\) and \(L_{\infty}\) errors from

\begin{table}
\begin{tabular}{c c c c c} \hline \(\varepsilon\) & \(N\) & \(err_{2}\) & \(err_{\infty}\) & \(err_{H_{1}}\) \\ \hline
0.5 & **17** & 3.100e-5 & 8.552e-5 & 1.987e-4 \\ \hline
0.1 & 38 & 6.718e-5 & 2.348e-4 & 8.560e-4 \\
0.05 & **66** & **8.311e-5** & **3.291e-4** & **7.847e-4** \\ \hline
0.01 & 186 & 1.086e-4 & 4.052e-4 & 8.433e-4 \\ \hline
0.005 & **367** & **7.481e-5** & **4.107e-4** & **1.017e-3** \\
0.002 & 750 & 2.541e-4 & 3.940e-4 & 1.899e-3 \\ \hline \end{tabular}
\end{table}
Table 1: Results of SRBFNN for Example 1. The second column records the number of basis functions in the final solution and the last three columns show the relative \(L_{2}\), \(L_{\infty}\) and \(H_{1}\) errors.

Figure 2: Numerical solution and its derivative obtained by SRBFNN and FDM when \(\varepsilon\)=0.01, 0.002 for Example 1.
SRBFNN are of \(10^{-4}\) for all scales and the relative \(H_{1}\) error is about \(10^{-3}\).

Example 2In this case, we consider the case that \(a^{\varepsilon}(x)\) is a two-scale function with scale separation

\[a^{\varepsilon}(x)=2+\sin(\frac{2\pi x}{\varepsilon})\cos(2\pi x). \tag{4.4}\]

We conduct a detailed comparison between SRBFNN and some deep neural network methods, including PINN [22], DGM [24], DRM [27] and MscaleDNN [13]. The training and test points are the same as those used in SRBFNN. PINN has five linear layers. DGM and DRM contain two residual blocks. The scale vector for MscaleDNN is [1, 2,..., \(\frac{1}{\varepsilon}\)]. The hidden layers are 16, 16, 32, 32, 64, 64 when \(\varepsilon=\)0.5, 0.1, 0.05, 0.01, 0.005, 0.002 for PINN, DGM and DRM. The hidden layers of MscaleDNN are (1000, 200, 150, 150, 100, 50, 50) for all \(\varepsilon\). \(tanh\) is used as the activation function for DRM, DGM and PINN. MscaleDNN employs a smooth and localized activation function \(s2ReLU\),

\[s2ReLU(x)=\sin(2\pi x)*ReLU(x)*ReLU(1-x).\]

Table 2 shows the results of relative \(L_{2}\) and \(H_{1}\) errors for 5 methods at different scales.

It is observed from Table 2 that more accurate results can be obtained by SRBFNN and MscaleDNN. This may be explained by the fact that SRBFNN and MscaleDNN use the activation functions with good local approximation capability, which is an essential factor for multiscale problems. Next, we compare the number of parameters in final solution for all methods and the average time per iteration, i.e., \(\dfrac{\text{Total running time}}{\text{Iterations}}\) in Table 3. It is found that SRBFNN contains fewer parameters and costs less time per iteration than other methods. This can be explained by the fact that SRBFNN has a simple network structure.

\begin{table}
\begin{tabular}{c|c c c c c|c c c c c} \hline  & & \multicolumn{4}{c|}{\(err_{2}\)} & \multicolumn{4}{c}{\(err_{H_{1}}\)} \\ \(\varepsilon\) & SRBFNN & MscaleDNN & DRM & DGM & PINN & SRBFNN & MscaleDNN & DRM & DGM & PINN \\ \hline
0.5 & 6.602e-6 & 2.531e-4 & 4.040e-2 & 1.697e-3 & 1.957e-3 & 1.972e-4 & 3.983e-3 & 1.437e-1 & 2.192e-2 & 2.122e-2 \\ \hline
0.1 & 3.59e-5 & 3.235e-4 & 5.237e-2 & 6.200e-3 & 6.168e-3 & 2.058e-4 & 1.111e-2 & 1.746e-1 & 4.876e-2 & 4.869e-2 \\ \hline
0.05 & 4.774e-5 & 2.743e-3 & 2.874e-2 & 4.466e-3 & 6.262e-3 & 3.368e-4 & 4.177e-2 & 1.020e-1 & 4.895e-2 & 4.932e-2 \\ \hline
0.01 & 2.01e-4 & 3.525e-3 & 2.874e-2 & 5.414e-3 & 5.549e-3 & 2.801e-3 & 4.694e-2 & 1.021e-1 & 4.924e-2 & 4.937e-2 \\ \hline
0.005 & 7.424e-5 & 3.695e-3 & 2.656e-2 & 7.424e-3 & 6.347e-3 & 1.282e-3 & 4.726e-2 & 9.597e-2 & 4.968e-2 & 4.935e-2 \\ \hline
0.002 & 1.654e-4 & 3.424e-3 & 2.661e-2 & 6.077e-3 & 5.511e-3 & 2.125e-3 & 4.711e-2 & 9.598e-2 & 4.914e-2 & 4.912e-2 \\ \hline \end{tabular}
\end{table}
Table 2: Comparison of SRBFNN, MscaleDNN, DRM, DGM, and PINN in terms of relative \(L_{2}\) and \(H_{1}\) errors for Example 2. The penalty coefficients of the boundary condition are 20.0, 20.0, 200.0 and 200.0 for PINN, DGM, DRM and MscaleDNN, respectively. The maximum number of iterations is 3000 for PINN, DGM, DRM and MscaleDNN.
Example 3Consider another two-scale function

\[a^{\varepsilon}(x)=2+\sin(2\pi x+\frac{2\pi x}{\varepsilon}). \tag{4.5}\]

Fig. 3 records the training processes of three relative errors and the number of RBFs when \(\varepsilon=0.002\). According to Fig. 3, we can observe that three relative errors stabilize around the 500th iteration, but they are still large at this point. When the regularization term was added in the loss function at around the 1000th iteration, the three relative errors all decreased significantly, which means that the \(ell_{1}\) regularization term can improve the generalization performance of the network.

Table 4 records the results of three relative errors and the number of RBFs in the final solution for Example 3. It is recognized that our method works well for two-scale coefficients. The accuracy and the number of RBFs in the final solution are also similar to the first two examples.

\begin{table}
\begin{tabular}{c|c c c c c c c c c} \hline \hline  & \multicolumn{4}{c}{Parameters} & \multicolumn{4}{c}{Average time per iteration (Seconds)} \\ \(\varepsilon\) & SRBFNN & MscaleDNN & DRM & DGM & PINN & SRBFNN & MscaleDNN & DRM & DGM & PINN \\ \hline
0.5 & 351 & 277251 & 1137 & 1137 & 865 & 0.122 & 0.714 & 0.143 & 0.877 & 0.792 \\ \hline
0.1 & 699 & 277751 & 1137 & 1137 & 865 & 0.140 & 0.746 & 0.142 & 0.872 & 0.754 \\ \hline
0.05 & 1062 & 277751 & 4321 & 4321 & 3265 & 0.147 & 0.764 & 0.147 & 0.843 & 0.796 \\ \hline
0.01 & 1977 & 277751 & 4321 & 4321 & 3265 & 0.144 & 0.793 & 0.139 & 0.861 & 0.777 \\ \hline
0.005 & 4107 & 277751 & 16833 & 16833 & 12673 & 0.149 & 0.825 & 0.171 & 0.799 & 0.837 \\ \hline
0.002 & 6474 & 277751 & 16833 & 16833 & 12673 & 0.150 & 0.734 & 0.176 & 0.841 & 0.835 \\ \hline \hline \end{tabular}
\end{table}
Table 3: Comparison of SRBFNN, MscaleDNN, DRM, DGM, and PINN in terms of the number of parameters and the average time per iteration measured in seconds for Example 2.

Figure 3: Three relative errors and the number of RBFs in terms of the iteration number when \(\varepsilon=0.002\) in Example 3.
\begin{table}
\begin{tabular}{c c c c c} \hline \(\varepsilon\) & \(N\) & \(err_{2}\) & \(err_{\infty}\) & \(err_{H_{1}}\) \\ \hline
0.5 & 18 & 4.850e-5 & 1.336e-4 & 2.319e-4 \\ \hline
0.1 & 40 & 6.539e-5 & 3.710e-4 & 5.730e-4 \\ \hline
0.05 & 69 & 6.778e-5 & 4.674e-4 & 4.249e-4 \\ \hline
0.01 & 178 & 1.086e-4 & 3.675e-4 & 1.414e-3 \\ \hline
0.005 & 379 & 1.431e-4 & 2.798e-4 & 2.134e-3 \\ \hline
0.002 & 760 & 9.549e-5 & 3.698e-4 & 1.500e-3 \\ \hline \end{tabular}
\end{table}
Table 4: Relative \(L_{2}\), \(L_{\infty}\) and \(H_{1}\) errors and the number of basis functions in the final solution for Example 3.
The relative \(L_{2}\) and \(H_{1}\) errors for five methods are shown in Table 5. The network architectures for MscaleDNN, DRM, DGM and PINN have the same settings as those used in Example 2.

As depicted in Table 5, SRBFNN obtains more accurate results than MscaleDNN, DRM, DGM, and PINN. Again, in the same settings, all the errors in Example 4 are

Figure 4: Numerical solution and its derivative obtained by SRBFNN and FDM when \(\varepsilon\)=0.05, 0.002 for Example 4.

\begin{table}
\begin{tabular}{c|c c c c c c|c c c c} \hline \hline  & \multicolumn{6}{c|}{\(err_{2}\)} & \multicolumn{6}{c}{\(err_{h}\)} \\ \(\varepsilon\) & SRBFNN & MscaleDNN & DRM & DGM & PINN & SRBFNN & MscaleDNN & DRM & DGM & PINN \\ \hline
0.5 & 3.673e-4 & 1.241e-3 & 1.668e-2 & 1.544e-2 & 1.698e-2 & 6.391e-3 & 2.353e-2 & 9.213e-2 & 1.121e-1 & 1.152e-1 \\ \hline
0.1 & 1.342e-4 & 1.375e-2 & 3.060e-2 & 3.712e-2 & 3.871e-2 & 5.174e-3 & 8.739e-2 & 1.583e-1 & 1.739e-1 & 1.740e-1 \\ \hline
0.05 & 1.185e-4 & 2.396e-2 & 3.092e-2 & 3.771e-2 & 3.898e-2 & 7.461e-3 & 1.273e-1 & 1.599e-1 & **1.744e-1** & **1.746e-1** \\ \hline
0.01 & 3.359e-4 & 3.127e-2 & 3.108e-2 & 3.622e-2 & 3.811e-2 & 1.765e-2 & 1.575e-1 & 1.595e-1 & 1.735e-1 & 1.375e-1 \\ \hline
0.005 & 6.714e-4 & 3.189e-2 & 3.104e-2 & 3.839e-2 & 3.846e-2 & 2.458e-2 & 1.587e-1 & 1.858e-1 & 1.730e-1 & 1.731e-1 \\ \hline
0.002 & 2.218e-4 & 3.201e-2 & 3.107e-2 & 3.661e-2 & 3.794e-2 & 3.832e-2 & 1.560e-1 & 1.557e-1 & 1.698e-1 & 1.705e-1 \\ \hline \hline \end{tabular}
\end{table}
Table 5: Relative \(L_{2}\) and \(H_{1}\) errors of SRBFNN, MscaleDNN, DRM, DGM, and PINN for Example 4. The penalty coefficients of the boundary condition and the maximum number of iterations are the same as Example 2 for MscaleDNN, DRM, DGM, and PINN.
larger than those in Example 2 and Example 3, while the discontinuity has less impact on SRBFNN. These results indicate that SRBFNN can handle multiscale problems with discontinuous coefficients.

Next, we plot the number of RBFs needed in the final solution and \(\varepsilon\) for Examples 1-4 in Fig. 5. It is observed that \(\ln(N)\) and \(\ln(\varepsilon^{-1})\) is approximatively linear and more RBFs are needed in Example 4 than in other examples. Table 6 estimates the slope by using the least squares method. From Table 6, we see that \(\ln(N)\sim 0.67\ln(\varepsilon^{-1})\), which implies \(N=O(\varepsilon^{0.67})\).

### Numerical experiments in 2D

**Setting:** In all two-dimensional cases, the initial learning rate is set to 0.1 and reduces by 1/10 every 40 iterations until it is less than \(10^{-5}\). The training data is equidistantly sampled with mesh size \(h\)=0.002. The mesh size of FDM is set to 0.0005. The \(batchsize\) is 1024 in the domain and 512\(\times\)4 on the boundary. \(MaxNiter=300\), \(SparseNiter=250\), \(tol_{1}=0.1\), \(tol_{2}=10^{-5}\), \(Check_{iter}=10\). The hyperparameters \(\lambda_{2}\) and \(\lambda_{3}\) are initialized to 20.0, 0.001 respectively. The solution is more oscillatory as \(\varepsilon\) decreases, and this makes the magnitude of derivatives larger. So to balance the effect of different loss terms, the

\begin{table}
\begin{tabular}{c c c c} \hline Example 1 & Example 2 & Example 3 & Example 4 \\ \hline
0.693 & 0.687 & 0.683 & 0.621 \\ \hline \end{tabular}
\end{table}
Table 6: The least squares estimation of slopes of \(\ln(N)\) with respect to \(\ln(\varepsilon^{-1})\) for Examples 1-4.

Figure 5: The number of basis functions needed in the final solution as a function of \(\varepsilon\) for Examples 1-4.
penalty coefficient \(\lambda_{1}\) is gradually reduced when \(\varepsilon\) goes to smaller as follows.

\[\lambda_{1}=\left\{\begin{array}{l}0.1,\,\mbox{if}\ \varepsilon=0.5,0.2,0.1\\ 0.02,\,\mbox{if}\ \varepsilon=0.05,0.02\\ 0.004,\,\mbox{if}\ \varepsilon=0.01\end{array}\right.. \tag{4.7}\]

Consider the problems in the domain \(\Omega=[0,1]^{2}\) with six scales, \(\varepsilon=\)0.5, 0.2, 0.1, 0.05, 0.02, 0.01 for all 2D examples. The initial number of RBFs is 1000, 1000, 2000, 5000, 15000, and 30000, respectively.

**Example 5**: \(a^{\varepsilon}(x,y)\) is defined by

\[a^{\varepsilon}(x,y)=2+\sin(\frac{2\pi(x+y)}{\varepsilon}), \tag{4.8}\]

where \(f(x,y)=-1.0\), \(g(x,y)=1.0\). Fig. 6 plots the training processes of three relative errors and the number of RBFs when \(\varepsilon=0.01\) for Example 5. One can see that SRBFNN achieves a good accuracy around the \(100th\) iteration and the number of basis functions no longer decreases at the \(250th\) iteration.

A detailed assessment of the predicted solution is presented in Fig. 7. In particular, it shows the absolute point-wise errors \(|u^{S}-u^{F}|\) at different scales. We can observe from Fig. 7 that the absolute point-wise errors in most areas are the order of O(\(10^{-3}\)), which means our method provides a good approximation for the solution. Then three relative errors and the number of basis functions in the final solution are recorded in Table 7. One can see that the number of basis functions grows much faster than that in 1D when \(\varepsilon\) becomes smaller. In addition, the accuracy of relative \(L_{2}\) and \(L_{\infty}\) errors is over one order of magnitude than that in 1D.

Figure 6: Training process of three relative errors and the number of RBFs when \(\varepsilon=0.01\) for Example 5.
Example 6[17]Consider \(a^{\varepsilon}(x,y)\) is a function with two scales

\[a^{\varepsilon}(x,y)=\frac{1.5+\sin(\frac{2\pi x}{\varepsilon})}{1.5+\sin( \frac{2\pi y}{\varepsilon})}+\frac{1.5+\sin(\frac{2\pi y}{\varepsilon})}{1.5+ \cos(\frac{2\pi x}{\varepsilon})}+\sin(4x^{2}y^{2})+1, \tag{4.9}\]

where \(f(x,y)=-10.0\), \(g(x,y)=0.0\). We conduct a comparison for SRBFNN with MscaleDNN, DRM, DGM, and PINN. The DRM and DGM contain four residual blocks. The PINN has nine linear layers. The hidden layers are 64, 64, 128, 128, 256, 256 when \(\varepsilon=0.5\), 0.2, 0.1, 0.05, 0.02, 0.01 for DRM, DGM and PINN. The hidden layers of MscaleDNN are (1000,400,300,300,200,100,100) for all \(\varepsilon\). \(batchsize\) is 1024 in the domain and 512\(\times\)4 on the

\begin{table}
\begin{tabular}{c c c c c} \hline \(\varepsilon\) & \(N\) & \(err_{2}\) & \(err_{\infty}\) & \(err_{H_{1}}\) \\ \hline
0.5 & 42 & 4.798e-4 & 1.596e-3 & 5.442e-3 \\ \hline
0.2 & 146 & 6.080e-4 & 1.328e-3 & 1.721e-2 \\
0.1 & 531 & **3.164e-4** & **1.453e-3** & **1.772e-2** \\
0.05 & 2121 & 3.561e-4 & 1.582e-3 & 2.684e-2 \\
0.02 & 6272 & **9.640e-4** & **1.862e-3** & **3.189e-2** \\
0.01 & 14164 & 7.719e-4 & 1.793e-3 & 2.070e-2 \\ \hline \end{tabular}
\end{table}
Table 7: Relative \(L_{2}\), \(L_{\infty}\) and \(H_{1}\) errors and the number of basis functions in the final solution for Example 5.

Figure 7: The absolute point-wise errors for different \(\varepsilon\) for Example 5.
boundary for these four methods. The initial learning rate is 0.001 and reduces by \(1/10\) every 100 iterations until it is less than \(10^{-5}\). We set the maximum number of iterations as 1000 for MsacleDNN, DRM, DGM, and PINN. Table 8 shows the relative \(L_{2}\) and \(H_{1}\) errors in terms of different scales for five methods.

As depicted in Table 8, SRBFNN reaches relatively lower errors than MscaleDNN, DRM, DGM, and PINN. PINN and DGM perform well when \(\varepsilon\) is big but the errors are relatively larger when \(\varepsilon\) is small. Additionally, one can observe that MscaleDNN and DRM outperform DGM and PINN, while MscaleDNN, DGM, PINN outperform DRM in one-dimensional cases. This observation is similar to the conclusion in [1]. Then we compare the number of parameters and the average time per iteration in all methods.

From Table 9, it is recognized that SRBFNN achieves better approximation than MscaleDNN, DRM, DGM, and PINN with fewer parameters. However, the average time per iteration for SRBFNN is larger than MscaleDNN, DRM, DGM, and PINN when \(\varepsilon\) is small, but SRBFNN requires fewer iterations to converge. Fig. 8 plots the final number of RBFs in terms of \(\varepsilon\). Analogously, we use the least squares method to estimate the slope of \(\ln(N)\) as a function of \(\ln(\varepsilon)\) for Example 5 and Example 6 in Table 10. One can observe from Fig. 8 and Table 10 that \(\ln(N)\sim 1.55\ln(\varepsilon^{-1})\), which implies \(N=O(\varepsilon^{-0.78*2})\). This observation, together with the results in aforementioned 1D examples can be written in a more general form: \(N=O(\varepsilon^{-\tau n})\), \(N\) is the number of RBFs in final solution, \(n=1\), \(2\) is the dimensionality and \(\tau\) is smaller than \(1\).

\begin{table}
\begin{tabular}{c|c c c c c|c c c c c} \hline \hline  & & & \(err_{2}\) & & & & & \(err_{H_{1}}\) & & \\ \(\varepsilon\) & SRBFNN & MscaleDNN & DRM & DGM & PINN & SRBFNN & MscaleDNN & DRM & DGM & PINN \\ \hline
0.5 & 1.319e-3 & 7.885e-3 & 2.062e-2 & 6.615e-3 & 8.812e-3 & 1.276e-2 & 5.979e-2 & 1.074e-1 & 1.005e-1 & 1.048e-1 \\
0.2 & 2.116e-3 & 8.467e-3 & 2.280e-2 & 7.393e-2 & 1.658e-1 & 1.650e-2 & 9.922e-2 & 1.654e-1 & 1.182e-1 & 0.303e-1 \\ \hline
0.1 & 3.099e-3 & 8.229e-3 & 2.601e-2 & 7.531e-1 & 4.333e-1 & 3.499e-2 & 1.116e-1 & 2.071e-1 & 7.626e-1 & 7.164e-1 \\ \hline
0.05 & 3.150e-3 & 3.096e-2 & 2.681e-2 & 9.676e-1 & 9.133e-1 & 7.774e-2 & 2.164e-1 & 2.182e-1 & 9.466e-1 & 9.231e-1 \\ \hline
0.02 & 3.959e-3 & 2.066e-2 & 3.546e-2 & 9.966e-1 & 1.000e-0 & 9.972e-2 & 2.019e-1 & 2.021e-1 & 9.958e-1 & 9.966e-1 \\ \hline
0.01 & 2.481e-3 & 2.589e-2 & 3.667e-2 & 9.982e-1 & 1.009e-0 & 1.125e-1 & 2.078e-1 & 2.392e-1 & 9.993e-1 & 9.997e-1 \\ \hline \hline \end{tabular}
\end{table}
Table 8: Relative \(L_{2}\) and \(H_{1}\) errors of SRBFNN, MscaleDNN, DRM, DGM, and PINN for Example 6. The penalty coefficient of the boundary condition for MscaleDNN and DRM is 500.0, and is 10.0 for PINN and DGM.

\begin{table}
\begin{tabular}{c|c c c c c|c c c c} \hline \hline  & & \multicolumn{4}{c}{Parameters} & \multicolumn{4}{c}{Average time per iteration (Seconds)} \\ \(\varepsilon\) & SRBFNN & MscaleDNN & DRM & DGM & PINN & SRBFNN & MscaleDNN & DRM & DGM & PINN \\ \hline
0.5 & 5240 & 704501 & 33537 & 33537 & 29377 & 11.102 & 9.731 & 6.948 & 11.725 & 10.087 \\ \hline
0.2 & 5760 & 704501 & 33537 & 33537 & 29377 & 10.747 & 10.315 & 6.722 & 11.796 & 10.111 \\ \hline
0.1 & 12690 & 704501 & 132609 & 132609 & 116097 & 10.889 & 9.885 & 6.722 & 11.821 & 10.233 \\ \hline
0.05 & 36375 & 704501 & 132609 & 132609 & 116097 & 15.310 & 9.713 & 6.751 & 11.783 & 10.185 \\ \hline
0.02 & 106810 & 704501 & 527361 & 527361 & 461569 & 37.869 & 9.821 & 6.827 & 11.704 & 10.032 \\ \hline
0.01 & 236425 & 704501 & 527361 & 527361 & 461569 & 74.965 & 8.231 & 6.756 & 11.942 & 10.080 \\ \hline \hline \end{tabular}
\end{table}
Table 9: The network parameters and the average running time per iteration in different methods for Example 6.
Example 7 ([17])Consider \(a^{\varepsilon}(x,y)\) having six scales

\[a^{\varepsilon}(x,y) = \frac{1}{6}(\frac{1.1+\sin(2\pi x/\varepsilon_{1})}{1.1+\sin(2\pi y /\varepsilon_{1})}+\frac{1.1+\sin(2\pi y/\varepsilon_{2})}{1.1+\cos(2\pi x/ \varepsilon_{2})}+\frac{1.1+\cos(2\pi x/\varepsilon_{3})}{1.1+\sin(2\pi y/ \varepsilon_{3})}+\] \[\frac{1.1+\sin(2\pi y/\varepsilon_{4})}{1.1+\cos(2\pi x/ \varepsilon_{4})}+\frac{1.1+\cos(2\pi x/\varepsilon_{5})}{1.1+\sin(2\pi y/ \varepsilon_{5})}+\sin(4x^{2}y^{2})+1),\]

where \(f(x,y)=-10.0\), \(g(x,y)=0.0\), \(\varepsilon_{1}=1/5\),\(\varepsilon_{2}=1/13\),\(\varepsilon_{3}=1/17\),\(\varepsilon_{4}=1/31\),\(\varepsilon_{5}=1/65\). Particularly, the initial learning rate is set to \(0.1\) and reduces to \(1/10\) every \(40\) iterations until it is less than \(10^{-6}\). The penalty parameter \(\lambda_{1}\) is \(0.02\) and the initial number of RBFs is \(30000\).

The absolute point-wise error \(|u^{S}-u^{F}|\) is shown in Fig. 9. SRBFNN can recover the solution up to \(10^{-3}\) point-wise accuracy. To further analyze the performance of SRBFNN, we draw the cross-section view of the solution and its derivatives in Fig. 10. Table 11 records three relative errors and the number of RBFs in the final solution. According to Fig. 10 and Table 11, we can observe that SRBFNN achieves a good approximation in terms of the solution and its derivatives.

\begin{table}
\begin{tabular}{l l} \hline Example 5 & Example 6 \\ \hline
**1.560** & **1.544** \\ \hline \end{tabular}
\end{table}
Table 10: The least squares estimation of slope about \(\ln(N)\) and \(\ln(\frac{1}{\varepsilon})\) for Example 5 and Example 6.

Figure 8: The number of basis functions in terms of \(\varepsilon\) in the log-log scale for Example 5 and 6.
### Numerical experiment in 3D

In this section, we consider a three-dimensional multiscale elliptic equation to verify the applicability of SRBFNN. In this scenario, solving multiscale elliptic equations with FDM is challenging. Thus we only consider the SRBFNN solution with a series of \(\varepsilon=0.5,0.2,0.1,0.05\) and take \(u^{0.05}\) as the reference.

**Setting:** The training data is sampled equidistantly with mesh size \(h\)=0.01. \(batchsize\) is 1024 in the domain and 200\(\times\)8 on the boundary. \(MaxNiter\) is 150 and \(SparseNiter\) is 120. The initial learning rate is 0.1 and reduces by 1/10 every 30 iterations. The initial

\begin{table}
\begin{tabular}{c c c c} \hline \(N\) & \(err_{2}\) & \(err_{\infty}\) & \(err_{H_{1}}\) \\ \hline
13213 & 4.427e-3 & 1.249e-2 & 1.581e-1 \\ \hline \end{tabular}
\end{table}
Table 11: Relative errors and the number of basis functions in the final solution for Example 7.

Figure 10: The cross-section view of \(u\), \(\frac{\partial u}{\partial x}\) and \(\frac{\partial u}{\partial y}\) for Example 7. The dotted line represents the FDM solution. The solid line with an asterisk is the SRBFNN solution.

Figure 9: The absolute point-wise error for Example 7.
number of RBFs is 1000, 2000, 5000, 10000, respectively. \(tol_{1}\!=\!0.05\), \(tol_{2}\!=\!10^{-5}\), \(Check_{iter}\!=\!10\). \(\lambda_{1}\),\(\lambda_{2}\),\(\lambda_{3}\) are initialized to 0.1, 50.0, 0.001, respectively.

Example 8Consider the coefficient \(a^{\varepsilon}(x,y,z)\)

\[a^{\varepsilon}(x,y,z)\!=\!2\!+\!\sin(\frac{2\pi x}{\varepsilon})\sin(\frac{2 \pi y}{\varepsilon})\sin(\frac{2\pi z}{\varepsilon}), \tag{4.10}\]

where \(f(x,y,z)\!=\!10.0\), \(g(x,y)\!=\!0.0\). Fig. 11 shows the slices of the SRBFNN solution when \(z\!=\!0.5\). We can see that the range of the solution is gradually the same as \(\varepsilon\) tends to 0.05.

Table 12 and Fig. 12 show the \(L_{2}\) norm between \(u^{0.05}\) and the solutions of other scales. As expected, \(u^{\varepsilon}\) approaches \(u^{0.05}\) and the number of basis functions increase with the scale \(\varepsilon\) tending to 0.05. These results illustrate our method can achieve numerical convergence for three-dimensional multiscale elliptic equations.
## 5 Conslusion

This paper proposes a novel method to solve multiscale elliptic problems with oscillatory coefficients, such as the coefficients with scale separation, discontinuity and multiple scales. Unlike the general DNNs, SRBFNN has a simpler and more efficient architecture, which makes the network converge faster. The activation functions with good local approximation are important for multiscale problems. Furthermore, the \(\ell_{1}\) regularization can deal with the overfitting problem owing to the simplicity of SRBFNN. The experiments show that SRBFNN presents the solution with fewer basis functions and has better approximation accuracy than most other deep learning methods. Finally, it is found that there is a relation between the number of basis functions \(N\) in the final solution and the scale \(\varepsilon\) in multiscale elliptic problems: \(N\!=\!O(\varepsilon^{-\tau n})\), \(n\) is the dimensionality, and \(\tau\) is typically smaller than 1, which is better than classical numerical methods.

There are several interesting issues which deserve further consideration. First, the choice of the initial number of RBFs is essential, affecting the network's accuracy and training. Second, a wise selection of penalty parameters in different loss terms facilitates the training process and provides a better approximation. Another issue is to find a faster sparse optimization algorithm.

## Acknowledgments

This work is supported in part by National Key R&D Program of China via grant No. 2022YFA1005203 and Jiangsu Provincial Key Research and Development Program under Grant BE2022058-4.

Figure 12: The error between the solution when \(\varepsilon\!=\!0.05\) and the solutions with larger scales and the relationship between the number of RBFs and \(\varepsilon\).