# Subspace Decomposition based DNN algorithm for elliptic type multi-scale
PDEs

Xi-An Li  lixian9131<EMAIL_ADDRESS>Zhi-Qin John Xu
<EMAIL_ADDRESS>Lei Zhang<EMAIL_ADDRESS>School of Mathematical
Sciences, Shanghai Jiao Tong University, Shanghai 200240, China  Institute of
Natural Sciences, Shanghai Jiao Tong University, Shanghai 200240, China  MOE-
LSC, Shanghai Jiao Tong University, Shanghai 200240, China

###### Abstract

While deep learning algorithms demonstrate a great potential in scientific
computing, its application to multi-scale problems remains to be a big
challenge. This is manifested by the “frequency principle” that neural
networks tend to learn low frequency components first. Novel architectures
such as multi-scale deep neural network (MscaleDNN) were proposed to alleviate
this problem to some extent. In this paper, we construct a subspace
decomposition based DNN (dubbed SD2NN) architecture for a class of multi-scale
problems by combining traditional numerical analysis ideas and MscaleDNN
algorithms. The proposed architecture includes one low frequency normal DNN
submodule, and one (or a few) high frequency MscaleDNN submodule(s), which are
designed to capture the smooth part and the oscillatory part of the multi-
scale solutions, respectively. In addition, a novel trigonometric activation
function is incorporated in the SD2NN model. We demonstrate the performance of
the SD2NN architecture through several benchmark multi-scale problems in
regular or irregular geometric domains. Numerical results show that the SD2NN
model is superior to existing models such as MscaleDNN.

Mathematics Subject Classifications: 52B10; 65D18; 68U05; 68U07

###### keywords:

Multi-scale; DNN; Fourier; Subspace-decomposed; Activation function

##  1 Introduction

Machine learning algorithms, especially deep neural networks (DNNs), have not
only achieved great success in traditional artificial intelligence tasks such
as image recognition, natural language processing, and recommendation systems
[1, 2], but also attracted more and more attention in the field of scientific
computation including the numerical solution of ordinary/partial differential
equations, integral-differential equations and dynamical systems [3, 4, 5, 6,
7, 8], with the fast development of novel computing devices, as well as the
rapidly increasing volume and complexity of data. In this paper, we will
introduce a new DNN-based algorithm to solve the following multi-scale
equation,

$\left\\{\begin{array}[]{cc}\mathcal{L}u(\bm{x})=f(\bm{x}),&\bm{x}\in\Omega,\\\
\mathcal{B}u(\bm{x})=g(\bm{x}),&\bm{x}\in\partial\Omega,\end{array}\right.$
(1.1)

where $\Omega$ is a bounded subset of $\mathbb{R}^{d}$ with piecewise
Lipschitz boundary which satisfies the interior cone condition. $\mathcal{L}$
is a linear or non-linear elliptic type differential operator on $\Omega$
which contains possibly non-separable multiple scales, and $\mathcal{B}$ is a
boundary operator on $\partial\Omega$. We assume that $\mathcal{L}$ is
uniformly strongly elliptic such that $\int_{\Omega}v^{\prime}\mathcal{L}v\geq
c\|v\|_{\mathcal{V}}$ for some admissible space $\mathcal{V}$, with a constant
$c>0$.

The multi-scale equation (1.1) has many physical and engineering applications,
such as heat conduction in composite materials, reservoir modeling in porous
media, convection dominated flows and the Poisson-Boltzmann model for
dielectric systems [9, 10, 11, 12], There has been amount of work concerning
the design of multi-scale numerical methods to achieve the optimal balance of
accuracy and complexity, such as homogenization [13, 14], heterogeneous multi-
scale methods (HMM) [15, 16], multi-scale network approximations [17], multi-
scale finite element methods (MsFEM)[18], variational multi-scale methods
(VMS) [19, 20], flux norm homogenization [21, 22], rough polyharmonic splines
(RPS) [23, 24], generalized multi-scale finite element methods (GMsFEM) [25,
26], localized orthogonal decomposition (LOD) [27, 28], etc. A common theme of
the multi-scale methods is to construct coarse approximation spaces with
optimal error control through the identification of low dimensional structures
in the high dimensional multi-scale solution space. The solution space can be
decomposed into a direct sum of a coarse space (with “smooth” components) and
(possibly a few) fine spaces (with “oscillatory” components), which is a
natural extension of the Fourier decomposition. The idea of subspace
decomposition has been widely used in numerical analysis and multi-scale
modeling [19, 27, 29, 30], furthermore, it has surprisingly deep connections
with Bayesian inference, kernel learning and probabilistic numerics [31, 32,
33, 34].

Comparing with conventional numerical algorithms, DNN approximation can
overcome the curse of dimensionality, therefore it is ideal for high
dimensional PDEs. On the other hand, it can be used as a meshless method which
is suitable for PDEs in complex domains. Various DNN based algorithms have
been proposed in [4, 35, 6, 36] to solve PDEs. Despite the above mentioned
progress, the frequency principle (F-Principle) [37, 38, 39, 40] shows that
general DNN-based algorithms often encounter a “curse of high-frequency” as
they are inefficient to learn high-frequency information of multiscale
functions. A series of subsequent theoretical investigations further confirm
such empirical observation [38, 41, 42, 43, 44, 45, 46, 47]. The study of the
F-Principle has also been utilized to understand various phenomena emerging in
applications of deep learning [48, 49, 50, 51, 52]. Inspired by the
F-Principle, a series of algorithms are developed to solve multi-scale PDEs
and to overcome the high-frequency curse of general DNNs [53, 54, 55, 56, 57,
58]. For example, in MscaleDNN, high frequency components are shifted into low
frequency ones by radial scaling such that they can be learnt more
efficiently. In addition, some smooth and localized activation functions were
proposed for MscaleDNN algorithm [56, 57].

In this paper, motivated by the subspace decomposition technique in numerical
analysis, and the MscaleDNN framework [56, 57], we propose a subspace-
decomposition based DNN (called SD2NN) architecture to solve the multi-scale
equation (1.1). The SD2NN framework consists of two parts: one low-frequency
or normal DNN submodule and one (or a few) MscaleDNN submodule(s), to capture
the low-frequency and high-frequency components of the multi-scale solution,
respectively. We enforce the orthogonality between subspaces in the SD2NN
architecture to enforce the stability of the decomposition, by adding a
penalty term to the loss function. Such a framework enables simultaneous fast
learning of both the smooth part and the oscillatory part of multiscale
functions/solutions, thus avoids the curse of frequency. We also improve the
activation function by adding a relaxation factor and using trigonometric
functions [59, 53], which will be elaborated in the later sections.

The paper is organized as follows. In Section 2, we briefly introduce the
subspace decomposition based DNN framework to approximate multi-scale
functions. In Section 3, we use the SD2NN framework to solve multi-scale
problems in the variational formulation, and discuss various options to choose
activation functions. Benchmark numerical experiments are carried out in
Section 4 to evaluate the performance of SD2NN. We conclude the paper in
Section 5.

##  2 Formulation of the Subspace Decomposition based DNN

###  2.1 Multi-scale DNN (MscaleDNN)

We briefly review the formulation of MscaleDNN introduced in [56, 57], for
convenience of readers. A deep neural network defines a function mapping:
$\bm{x}\in\mathbb{R}^{d}\rightarrow y_{\bm{\theta}}(\bm{x})\in\mathbb{R}$,
where $d$ is the dimension of input. In particular, the DNN function is a
nested composition of linear functions and nonlinear activation functions,
which is of the form

$\begin{cases}\bm{y}_{\bm{\theta}}^{[0]}=\bm{x}\\\
\bm{y}_{\bm{\theta}}^{[\ell]}=\sigma\circ(\bm{W}^{[\ell]}\bm{y}_{\bm{\theta}}^{[\ell-1]}+\bm{b}^{[\ell]}),~{}~{}\text{for}~{}~{}\ell=1,2,3,\cdots\cdots,L\end{cases}$
(2.1)

where $\bm{W}^{[\ell]}\in\mathbb{R}^{m_{\ell+1}\times
m_{\ell}},\bm{b}^{[\ell]}\in\mathbb{R}^{m_{\ell+1}}$ are the weight matrix and
bias vector of the $\ell$-th hidden layer, respectively, $m_{0}=d$ and
$m_{L+1}$ is the dimension of the output, $``\circ"$ stands for the
elementary-wise operation. $\sigma(\cdot)$ is an element-wise activation
function. We also denote the output of DNNs by $\bm{y}(\bm{x};\bm{\theta})$
with $\bm{\theta}$ standing for the parameter set
$\left(\bm{W}^{[1]},\cdots\bm{W}^{[L]},\bm{b}^{[1]},\cdots\bm{b}^{[L]}\right)$.

A normal DNN is usually inadequate to solve multi-scale PDEs as it tends to
stagnate on low frequency components. The MscaleDNN architecture was proposed
[56, 57] based on the frequency principle to alleviate this problem. It is
illustrated in Figure 1 and described in the following: We first divide the
neurons in the first hidden-layer into $Q$ parts, and construct the following
vector

$\Lambda=(\overbrace{a_{1},\cdots,a_{1}}^{N_{1}},\overbrace{a_{2},\cdots,a_{2}}^{N_{2}},\cdots,\overbrace{a_{Q},\cdots,a_{Q}}^{N_{Q}})^{T}$
(2.2)

where $\\{a_{q}\\}_{q=1}^{Q}$ are positive scale-factors (usually in ascending
order), and $\\{N_{q}\\}_{q=1}^{Q}$ are the number of the duplicated $q$-th
scale-factor. The input data $\bm{x}$ is scaled by $\Lambda$ such that
$\tilde{\bm{x}}=\Lambda\odot(\bm{W}^{1}\bm{x}+\bm{b}^{1})$, then
$\tilde{\bm{x}}$ is fed into the subsequent module of MscaleDNN, here and
hereinafter $\odot$ stands for the element-wise product.

$\bm{x}$input$a_{Q}x$$\cdots$$\cdots$$a_{2}x$$a_{1}x$DNN
$y(\bm{x};\bm{\theta})$ Figure 1: A schematic diagram for MScaleDNN with
$N_{i}=1$, $i=1,\cdots,Q$.

From the viewpoint of function approximation, the first layer of the MscaleDNN
model can be regarded as a series of basis functions with various scales and
the output of MscaleDNN is the (nonlinear) combination of those basis
functions [56, 57, 53].

We use the residual neural network (ResNet) [60] to overcome the vanishing
gradient phenomenon in backpropagation by introducing skip connections between
nonadjacent layers. For example, the ResNet unit with one-step connection
produces a filtered version $\bm{y}^{[\ell+1]}(\bm{x};\bm{\theta})$ for the
input $\bm{y}^{[\ell]}(\bm{x};\bm{\theta})$ is as follows

$\bm{y}^{[\ell+1]}(\bm{x};\bm{\theta})=\bm{y}^{[\ell]}(\bm{x};\bm{\theta})+\sigma\circ\bigg{(}\bm{W}^{[\ell+1]}\bm{y}^{[\ell]}(\bm{x};\bm{\theta})+\bm{b}^{[\ell+1]}\bigg{)}.$

ResNet can accelerate the training process and improve the performance of
DNNs. For scientific computation tasks, it can help improve the capability of
DNNs to approximate high-order derivatives and solutions of PDEs [61, 62].

###  2.2 Subspace decomposition

It is natural to decompose multi-scale solutions into smooth components and
oscillatory components. For example, in periodic homogenization with
$\mathcal{L}=-\nabla\cdot a(x/\epsilon)\nabla$, the solution $u^{\epsilon}$ of
$\mathcal{L}u^{\epsilon}=f$ can be approximated by $u_{0}+\epsilon u_{1}$,
such that $u_{0}$ is the solution of the homogenized equation
$\mathcal{L}^{0}u_{0}=f$, where $\mathcal{L}^{0}=-\nabla\cdot A\nabla$ with
$A$ being the homogenized coefficient. $u^{0}$ is much smoother than
$u^{\epsilon}$, and
$u_{1}=\epsilon\sum_{i=1}^{d}\chi_{i}(x/\epsilon)\frac{\partial
u_{0}}{\partial x_{i}}$ characterizes the oscillations. $\chi_{i}$ is the so-
called corrector of the homogenization problem or solution of the cell problem
[63, 14], such that

$\left\\{\begin{array}[]{c}-\nabla_{y}\cdot
a(y)(\nabla_{y}\xi_{i}(y)+e_{i})=0,\\\ \xi_{i}\text{ periodic in
}Y.\end{array}\right.$ (2.3)

where $y=x/\epsilon$, and $Y$ is the unit periodic cell.

In many numerical homogenization type multi-scale methods such as VMS, RPS and
LOD methods [19, 20, 23, 24, 27, 28, 64], a coarse solution is computed in a
low dimensional approximation space $V_{c}$ with quasi-optimal approximation,
stability and localization properties. The bases in $V_{c}$ contain fine scale
information, and they can be pre-computed in localized subdomains in parallel.
For some multi-scale problems, numerical homogenization methods can be proved
to achieve guaranteed quantitative error estimate with respect to coarse
resolutions. Furthermore, multi-scale methods can be considered as two level
numerical methods, and can be generalized as multilevel subspace decomposition
to construct multigrid type preconditioners for the efficient resolution of
fine scale problem, see [65, 30].

##  3 SD2NN model to multi-scale problems and the options for activation
function

###  3.1 Unified SD2NN architecture to solve multi-scale problems

We introduce subspace-decomposed DNN (called SD2NN) architecture in the
following. For simplicity, we employ the deep Ritz method [61] for the
solution of multi-scale PDEs (1.1). Other architectures such as PINN [66],
deep Galerkin [35] etc. can be adapted.

##### Continuous variational formulation

The solution of (1.1) can be obtained by minimizing the following Dirichlet
energy

$\mathcal{J}(v)=\frac{1}{2}\int_{\Omega}v^{\prime}\mathcal{L}vd\bm{x}-\int_{\Omega}fvd\bm{x},$
(3.1)

where $v=v(\bm{x})\in\mathcal{V}$ is a trial function, where $\mathcal{V}$ is
the admissible function space for $v$. We are looking for the solution

$u=\arg\min_{v\in\mathcal{V}}\mathcal{J}(v).$

Suppose that the solution $u$ has the following coarse/fine decomposition,
$u=u_{c}+u_{f}$, in which $u_{c}$ contains the coarse shape and $u_{f}$
contains the fine details of the multi-scale solution $u$, respectively.
Formally, (3.1) can be rewritten as

$\mathcal{J}(v_{c},v_{f})=\frac{1}{2}\int_{\Omega}(v_{c}+v_{f})^{\prime}\mathcal{L}(v_{c}+v_{f})d\bm{x}-\int_{\Omega}f(v_{c}+v_{f})d\bm{x}.$
(3.2)

To ensure the well-posedness of the variational problem with respect to
$v_{c}$ and $v_{f}$, we need to define the respective subspaces
$\mathcal{V}_{c}$ and $\mathcal{V}_{f}$, in where $v_{c}$ and $v_{f}$ ”live”,
which will be specified in the particular examples. The coarse part $v_{c}$
can be represented by a low-frequency or a normal DNN
$y_{1}(\cdot,\bm{\theta}_{1})$, and the fine-part (or high-frequency part)
$u_{f}$ can be represented by a multi-scale DNN
$y_{2}(\cdot,\bm{\theta}_{2})$. Here, $\bm{\theta}_{1}\in\Theta_{1}$ and
$\bm{\theta}_{2}\in\Theta_{2}$ denote the parameters of the underlying DNNs.
In this case, we have

$\begin{array}[]{cc}\mathcal{V}_{c}&=\bigg{\\{}y_{1}(\bm{x};\bm{\theta}_{1})\bigg{|}y_{1}(\bm{x};\bm{\theta}_{1})=g,\bm{x}\in\partial\Omega;\bm{\theta}_{1}\in\Theta_{1}\bigg{\\}},\\\
\mathcal{V}_{f}&=\bigg{\\{}y_{2}(\bm{x};\bm{\theta}_{2})\bigg{|}y_{2}(\bm{x};\bm{\theta}_{2})=0,\bm{x}\in\partial\Omega;\bm{\theta}_{2}\in\Theta_{2}\bigg{\\}}.\end{array}$
(3.3)

######  Remark.

Here we present the SD2NN model as a two level model, and it can be easily
generalized to a multilevel model by adding more MscaleDNN submodules with
appropriate frequency factors. We will illustrate this point through Example 2
in the numerics section.

We further introduce a hyper-parameter $\alpha>0$ to control the contribution
of fine-part, i.e.,
$y(\bm{x};\bm{\theta}_{1},\bm{\theta}_{2})=y_{1}(\bm{x};\bm{\theta}_{1})+\alpha
y_{2}(\bm{x};\bm{\theta}_{2})$. Let

$\displaystyle\mathcal{J}_{\alpha}\bigg{(}y_{1}(\bm{x};\bm{\theta}_{1}),y_{2}(\bm{x};\bm{\theta}_{2})\bigg{)}$
$\displaystyle=\frac{1}{2}\int_{\Omega}\bigg{(}y_{1}(\bm{x};\bm{\theta}_{1})+\alpha
y_{2}(\bm{x};\bm{\theta}_{2})\bigg{)}^{\prime}\mathcal{L}\bigg{(}y_{1}(\bm{x};\bm{\theta}_{1})+\alpha
y_{2}(\bm{x};\bm{\theta}_{2})\bigg{)}d\bm{x}$ (3.4)
$\displaystyle-\int_{\Omega}f\bigg{(}y_{1}(\bm{x};\bm{\theta}_{1})+\alpha
y_{2}(\bm{x};\bm{\theta}_{2})\bigg{)}d\bm{x},$

we then obtain the following variational problem

$u_{c},u_{f}=\underset{y_{1}(\bm{x},\bm{\theta}_{1})\in\mathcal{V}_{c},y_{2}(\bm{x},\bm{\theta}_{2})\in\mathcal{V}_{f}}{\arg\min}~{}\mathcal{J}_{\alpha}\bigg{(}y_{1}(\bm{x};\bm{\theta}_{1}),y_{2}(\bm{x};\bm{\theta}_{2})\bigg{)}$

##### Discrete Variational Formulation

The integral in (3.4) can be discretized by Monte Carlo method [67], namely,
we define

$L_{\mathcal{J}_{\alpha}}(S_{I};\bm{\theta}_{1},\bm{\theta}_{2})=\frac{1}{n_{in}}\sum_{i=1}^{n_{in}}\Bigg{[}\frac{1}{2}\bigg{(}y_{1}(\bm{x}_{I}^{i};\bm{\theta}_{1})+\alpha
y_{2}(\bm{x}_{I}^{i};\bm{\theta}_{2})\bigg{)}^{\prime}\mathcal{L}\bigg{(}y_{1}(\bm{x}_{I}^{i};\bm{\theta}_{1})+\alpha
y_{2}(\bm{x}_{I}^{i};\bm{\theta}_{2})\bigg{)}-f\bigg{(}y_{1}(\bm{x}_{I}^{i};\bm{\theta}_{1})+\alpha
y_{2}(\bm{x}_{I}^{i};\bm{\theta}_{2})\bigg{)}\Bigg{]},$

here and hereinafter $S_{I}$ stands for the sampling points with uniform
distribution in $\Omega$, and

$\bm{\theta}_{1}^{*},\bm{\theta}_{2}^{*}=\underset{\bm{\theta}_{1}\in{\Theta}_{1},\bm{\theta}_{2}\in{\Theta}_{2}}{\arg\min}~{}L_{\mathcal{J}_{\alpha}}(S_{I};\bm{\theta}_{1},\bm{\theta}_{2}),$

such that $u_{c}=y_{1}(\bm{x};\bm{\theta}_{1}^{*})$, $u_{f}=\alpha
y_{2}(\bm{x};\bm{\theta}_{2}^{*})$.

##### Orthogonality constraints

In order to separate the coarse-part and fine-part of solution, we add the
$L^{2}$ orthogonality constraint
$\int_{\Omega}y_{1}(\bm{x};\bm{\theta})\cdot\alpha
y_{2}(\bm{x};\bm{\theta})d\bm{x}=0$ as a penalty term to the loss function,
namely, let

$L_{orth}(S_{I};\bm{\theta}_{1},\bm{\theta}_{2})=\bigg{|}\frac{1}{n_{in}}\sum_{i=1}^{n_{in}}y_{1}(\bm{x}_{I}^{i};\bm{\theta}_{1})\cdot\alpha
y_{2}(\bm{x}_{I}^{i};\bm{\theta}_{2})\bigg{|}^{2},\bm{x}_{I}^{i}\in S_{I}.$
(3.5)

##### Boundary condition

Boundary conditions are important constraints for numerical solution of PDEs.
According to the definition of coarse and fine spaces (3.3), we have the
following penalties for the coarse part
$y_{1}\big{(}\bm{x};\bm{\theta}_{1}\big{)}$ and the fine-part
$y_{2}\big{(}\bm{x};\bm{\theta}_{2}\big{)}$,

$L_{bdc}(S_{B};\bm{\theta}_{1})=\frac{1}{n_{bd}}\sum_{j=1}^{n_{bd}}\bigg{[}\mathcal{B}y_{1}\big{(}\bm{x}_{B}^{j};\bm{\theta}_{1}\big{)}-g(\bm{x}_{B}^{j})\bigg{]}^{2}~{}~{}\text{for}~{}~{}\bm{x}_{B}^{j}\in
S_{B},$ (3.6)

and

$L_{bdf}(S_{B};\bm{\theta}_{2})=\frac{1}{n_{bd}}\sum_{\mathcal{}=1}^{n_{bd}}\bigg{[}\mathcal{B}\alpha
y_{2}\big{(}\bm{x}_{B}^{j};\bm{\theta}_{2}\big{)}-0\bigg{]}^{2}~{}~{}\text{for}~{}~{}\bm{x}_{B}^{j}\in
S_{B}.$ (3.7)

For comparison, we also introduce the penalty for the ”unified” boundary
condition of
$y\big{(}\bm{x};\bm{\theta}_{1},\bm{\theta}_{2}\big{)}=y_{1}(\bm{x};\bm{\theta}_{1})+\alpha
y_{2}(\bm{x};\bm{\theta}_{2})$, i.e.,

$L_{bdu}(S_{B};\bm{\theta}_{1})=\frac{1}{n_{bd}}\sum_{j=1}^{n_{bd}}\bigg{[}\mathcal{B}y\big{(}\bm{x}_{B}^{j};\bm{\theta}_{1},\bm{\theta}_{2}\big{)}-g(\bm{x}_{B}^{j})\bigg{]}^{2}~{}~{}\text{for}~{}~{}\bm{x}_{B}^{j}\in
S_{B},$ (3.8)

with $S_{B}$ standing for the collection of sampling points with uniform
distribution on $\partial\Omega$.

We conclude with the following loss function:

$L_{\alpha}({S_{I},S_{B}};\bm{\theta}_{1},\bm{\theta}_{2})=\underbrace{L_{\mathcal{J}_{\alpha}}(S_{I};\bm{\theta}_{1},\bm{\theta}_{2})}_{loss\\_in}+\underbrace{\gamma\bigg{(}L_{bdc}(S_{B};\bm{\theta}_{1})+L_{bdf}(S_{B};\bm{\theta}_{2})\bigg{)}}_{loss\\_bd}+\underbrace{\beta
L_{orth}(S_{I};\bm{\theta}_{1},\bm{\theta}_{2})}_{loss\\_orth}$ (3.9)

with $S_{I}=\\{\bm{x}_{I}^{i}\\}_{i=1}^{n_{it}}$ and
$S_{B}=\\{\bm{x}_{B}^{j}\\}_{j=1}^{n_{bd}}$ being the sets of uniformly
distributed sample points on $\Omega$ and $\partial\Omega$, respectively. The
first term $loss\\_{in}$ minimizes the residual of the PDE, the second term
$loss\\_{bd}$ enforces the given boundary condition, and the third term
$loss\\_dot$ imposes the orthogonality constraint. In addition, we introduce
two penalty parameter $\gamma$ and $\beta$ to control the contribution of
$loss\\_bd$ and $loss\\_orth$ for loss function, in which $\gamma$ is
increasing gradually during the training process [68] and $\beta$ is a fixed
constant.

Our goal is to find two sets of parameter $\bm{\theta}_{1},\bm{\theta}_{2}$
which minimize the loss function
$L_{\alpha}({S_{I},S_{B}};\bm{\theta}_{1},\bm{\theta}_{2})$, i.e.,

$\bm{\theta}_{1}^{*},\bm{\theta}_{2}^{*}=\arg\min
L_{\alpha}(S_{I},S_{B};\bm{\theta}_{1},\bm{\theta}_{2})~{}~{}\Longrightarrow~{}~{}u_{c}(\bm{x}),u_{f}(\bm{x})=y_{1}(\bm{x};\bm{\theta}_{1}^{*}),\alpha
y_{2}(\bm{x};\bm{\theta}_{2}^{*}).$

If $L_{\alpha}({S_{I},S_{B}};\bm{\theta}_{1},\bm{\theta}_{2})$ is small
enough, then $y(\bm{x};\bm{\theta}_{1},\bm{\theta}_{2})$ will be very close to
the solution of (1.1). The parameter $\bm{\theta}_{1}$ and $\bm{\theta}_{2}$
can be computed by stochastic gradient descent method with a pre-defined
learning rate schedule. We use the Xavier initialization for the weights and
biases [69, 70]. In other word, the weights and biases are sampled from normal
distribution $\mathcal{D}=\mathcal{N}\left(0,\frac{2}{m_{in}+m_{out}}\right)$,
where $m_{in}$ and $m_{out}$ are the input and output dimensions of the
corresponding layer, respectively.

###  3.2 Options for activation function

The choice of activation functions is crucial for the performance of DNN based
algorithms. Two localized activation functions, namely,
$\text{sReLU}(\bm{x})=\text{ReLU}(\bm{x})*\text{ReLU}(1-\bm{x})$ and
$\text{s2ReLU}(\bm{x})=\sin(2\pi\bm{x})*\text{ReLU}(\bm{x})*\text{ReLU}(1-\bm{x})$
were proposed for MscaleDNN in [56, 57]. Heuristic analysis and numerical
results show that the latter one is smoother and more robust. However, these
two activation functions are supported in range $[0,1]$, therefore the output
will become zero if the input is outside the interval $[0,1]$. This will
affect the performance of MscaleDNN in the training process.

In this work, we propose the following _soften Fourier mapping_ (_SFM_)
activation function inspired by the trigonometric activation function proposed
in [59, 54, 53], namely,

$\sigma(\bm{z})=s\times\left[\begin{array}[]{c}\cos(\bm{z})\\\
\sin(\bm{z})\end{array}\right],$

where the _relaxation_ parameter $s\in(0,1]$ is used to control the range of
output. Empirically, we find that $s=0.5$ is a good choice, see Fig. 2. With
the SFM activation function, the MscaleDNN can be regard as a pipeline of
Fourier-like transform. The first hidden layer of the MscaleDNN architecture
mimic the Fourier expansion, and the remaining structure learns the
approximate Fourier coefficients, which are often relatively less oscillate
with respect to the input. In this sense, the learning can be effectively
accelerated.

The above mentioned activation functions are illustrated in Fig. 3.

Figure 2: sReLU and s2ReLU functions (left), SFM functions with $s=0.5$
(right).

We conclude the section with the schematic of the SD2NN architecture in Figure
3. The output of SD2NN is obtained by a linear layer with no activation
function.

$\bm{x}_{I}$$\bm{x}_{B}$$\bm{x}$$\int$$\int$$\vdots$$\int$$\int$$\int$$\int$$\vdots$$\int$$\int$Low-
freqency or normal NN$\int$$\int$$\vdots$$\int$$\int$
$y_{1}(\bm{x};\bm{\theta}_{1})$ $\nabla
y_{1}(\bm{x};\bm{\theta}_{1})$$f(\bm{x})$$L_{\mathcal{J}_{\alpha}}\bigg{(}y_{1}(\bm{x};\bm{\theta}_{1}),y_{2}(\bm{x};\bm{\theta}_{2})\bigg{)}$$|y_{1}(\bm{x};\bm{\theta}_{1})-g(\bm{x})|^{2}~{}$
$\begin{matrix}s\cdot\sin(W^{[1]}_{1}\tilde{\bm{x}})\\\
s\cdot\cos(W^{[1]}_{1}\tilde{\bm{x}})\end{matrix}$$\vdots$$\begin{matrix}s\cdot\sin(W^{[1]}_{Q}\tilde{\bm{x}})\\\
s\cdot\cos(W^{[1]}_{Q}\tilde{\bm{x}})\end{matrix}$$\tilde{\bm{x}}=a_{1}\bm{x}$$\tilde{\bm{x}}=a_{Q}\bm{x}$$\int$$\int$$\vdots$$\int$$\int$High-
freqency NN$\int$$\int$$\vdots$$\int$$\int$ $y_{2}(\bm{x};\bm{\theta}_{2})$
$\nabla y_{2}(\bm{x};\bm{\theta}_{2})$$|\alpha
y_{2}(\bm{x};\bm{\theta}_{2})-0|^{2}~{}$
loss$\bm{\theta}_{1}^{*},\bm{\theta}_{2}^{*}$minimize$\alpha^{2}\big{|}y_{1}(\bm{x};\bm{\theta}_{1})\cdot
y_{2}(\bm{x};\bm{\theta}_{2})\big{|}^{2}$ Figure 3: Schematic of SD2NN for
solving multi-scale problems

##  4 Numerical experiments

In this section, we test the performance of SD2NN method for multi-scale
problems with highly oscillatory coefficients and Poisson-Boltzmann equation
in perforated domains. We demonstrate that the SD2NN method can efficiently
and simultaneously solve the coarse part and the fine part of the multi-scale
solution. We also investigated the effects of different activation functions
(for example, SFM and s2ReLU). We show that SD2NN has better performance
compared with existing methods such as MscaleDNN in [57, 53].

###  4.1 Model and training setup

####  4.1.1 Model setup

The details of all the models in the numerical experiments are elaborated in
the following.

  * 1.

_SD 2NN1_: We use a normal fully-connected DNN as the coarse submodule and a
high-frequency MscaleDNN as the fine submodule. The first layer of the fine
submodule has scaling factors $\Lambda=(21,22,23,\cdots,120)$ with $Q=100$, as
in (2.2). The activation function for each hidden layer in the normal DNN is
chosen as $\rm{tanh}(x)$, the activation functions for the first hidden layer
and the subsequent layers of MscaleDNN are chosen as SFM with $s=0.5$ and
s2ReLU, respectively.

  * 2.

_SD 2NN2_: We use a low-frequency MscaleDNN as the coarse submodule and a
high-frequency MscaleDNN as the fine submodule. The first layers of coarse
submodule and fine submodule have scaling factors
$\Lambda_{c}=(0.5,1,1.5,2,\cdots,19.5,20)$ with $Q=40$, and
$\Lambda_{f}=(21,22,23,\cdots,120)$ with $Q=100$, respectively. The activation
function for their first hidden layers are all chosen as SFM with $s=1.0$ and
$s=0.5$(denoted as SFM(1.0) and SFM(0.5)), but the activation functions for
their subsequent layers are chosen as ${\rm tanh}(x)$ and s2ReLU,
respectively.

  * 3.

_SD 2NN3_: The activation functions and scaling factors are same as _SD 2NN2_.
But we take away the orthogonality penalty term (3.9) from the loss function,
and the boundary condition is the unified one in (3.8).

  * 4.

_Mscale_ : A MscaleDNN model [57] with s2ReLU activation function for all
layers. The neurons of the first layer are divided into $Q=120$ groups with
scaling factors $\Lambda=(1,2,3,\cdots,119,120)$ in (2.2).

  * 5.

_WWP_ : A MscaleDNN model [53] with hybrid activation function including SFM
with $s=1.0$ (denoted as SFM(1.0)) and ${\rm tanh}(x)$. The scaling factors of
the first layer are divided as four subgroups. Each subgroup has 30 samples
from Gaussian distribution $\mathcal{N}(0;\tau_{i}^{2})$ with $\tau_{1}=1$,
$\tau_{1}=20$, $\tau_{1}=50$ and $\tau_{4}=100$. This model is proposed by
Wang, Wang, and Perdikaris in [53], and we denote it as _WWP_.

Table 1:  Comparisons for the above models Model | Submodules | Activation | Scale factor | Boundary | Orthogonality
---|---|---|---|---|---
SD2NN1 | Normal | tanh | —– | individual boundary | Yes
MscaleDNN | SFM(0.5)+s2ReLU | $(21,22,23,\cdots,120)$ | (3.6) and (3.7)
SD2NN2 | MscaleDNN | SFM(1.0)+tanh | $(0.5,1,1.5,2,\cdots,19.5,20)$ | individual boundary | Yes
MscaleDNN | SFM(0.5)+s2ReLU | $(21,22,23,\cdots,120)$ | (3.6) and (3.7)
SD2NN3 | MscaleDNN | SFM(1.0)+tanh | $(0.5,1,1.5,2,\cdots,19.5,20)$ | unified boundary (3.8) | No
MscaleDNN | SFM(0.5)+s2ReLU | $(21,22,23,\cdots,120)$
Mscale | —– | s2ReLU | $(1,2,3,\cdots,119,120)$ | unified boundary (3.8) | —–
WWP | —– | SFM(1.0)+tanh | $\begin{matrix}\Lambda=[\Lambda_{1};\Lambda_{2};\Lambda_{3};\Lambda_{4}]~{}~{}\textup{with}\\\ \Lambda_{i}\sim\mathcal{N}(0,\tau_{i}),~{}~{}\tau_{1}=1,\\\ \tau_{2}=20,\tau_{3}=50,\tau_{4}=100\end{matrix}$ | unified boundary (3.8) | —–

####  4.1.2 Training setup

We use the relative square error to evaluate the accuracy of different models:

$REL=\sum_{i=1}^{N^{\prime}}\frac{|\tilde{u}(\bm{x}^{i})-u^{*}(\bm{x}^{i})|^{2}}{|u^{*}(\bm{x}^{i})|^{2}}$

where $\tilde{u}(\bm{x}^{i})$ and $u^{*}(\bm{x}^{i})$ are the approximate DNN
solution and the exact solution, respectively,
$\\{\bm{x}^{i}\\}_{i=1}^{N^{\prime}}$ are testing points, and $N^{\prime}$ is
the number of testing points.

In our numerical experiments, all training and testing data are sampled
uniformly in $\Omega$ (or $\partial\Omega$), and all networks are trained by
Adam optimizer. The initial learning rate is set as $2\times 10^{-4}$ with a
decay rate $5\times 10^{-5}$ for each training epoch. For visualization of the
training process, we test our model every 1000 epochs in the training process.
The penalty parameter $\beta$ for the orthogonality constraint (3.5) is set as
20, $\gamma$ for the boundary constraint (3.6) and (3.7) is set as

$\gamma=\left\\{\begin{aligned}
\gamma_{0},\quad&\textup{if}~{}~{}i_{\textup{epoch}}<0.1T_{\max}\\\
10\gamma_{0},\quad&\textup{if}~{}~{}0.1T_{\max}<=i_{\textup{epoch}}<0.2T_{\max}\\\
50\gamma_{0},\quad&\textup{if}~{}~{}0.2T_{\max}<=i_{\textup{epoch}}<0.25T_{\max}\\\
100\gamma_{0},\quad&\textup{if}~{}~{}0.25T_{\max}<=i_{\textup{epoch}}<0.5T_{\max}\\\
200\gamma_{0},\quad&\textup{if}~{}~{}0.5T_{\max}<=i_{\textup{epoch}}<0.75T_{\max}\\\
500\gamma_{0},\quad&\textup{otherwise}\end{aligned}\right.$ (4.1)

where $\gamma_{0}=100$ in all our tests and $T_{\max}$ represents the total
epoch number. We implement our code in TensorFlow (version 1.14.0) on a work
station (256-GB RAM, single NVIDIA GeForce GTX 2080Ti 12-GB).

###  4.2 Numerical examples and results

In this section, we test the SD2NN models for multi-scale diffusion and
Poisson-Boltzmann equations.

######  Example 1 (1D multi-scale elliptic problem).

We use this problem to benchmark our models. Let us consider the following
multi-scale elliptic equation in domain $\Omega=[a,b]^{d}$

$\begin{cases}-\textrm{div}\bigg{(}A(\bm{x})|\nabla u(\bm{x})|^{p-2}\nabla
u(\bm{x})\bigg{)}=f(\bm{x}),~{}~{}\bm{x}\in\Omega,\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}u(\bm{x})=g(\bm{x}),~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\bm{x}\in\partial\Omega.\end{cases}$
(4.2)

which has the following Dirichlet energy

$\mathcal{J}(v)=\frac{1}{p}\int_{\Omega}A|\nabla
v|^{p}d\bm{x}-\int_{\Omega}fvd\bm{x}.$ (4.3)

Linear Case:

We first consider the linear case $p=2$ with $\Omega=[0,1]$ and boundary
condition $u(0)=u(1)=0$,

$A(x)=\left(2+\cos\left(2\pi\frac{x}{\epsilon}\right)\right)^{-1}$ (4.4)

with a small parameter $\epsilon>0$ such that
$\epsilon^{-1}\in\mathbb{N}^{+}$, the unique solution is given by

$u(x)=x-x^{2}+\epsilon\left(\frac{1}{4\pi}\sin\left(2\pi\frac{x}{\epsilon}\right)-\frac{1}{2\pi}x\sin\left(2\pi\frac{x}{\epsilon}\right)-\frac{\epsilon}{4\pi^{2}}\cos\left(2\pi\frac{x}{\epsilon}\right)+\frac{\epsilon}{4\pi^{2}}\right).$
(4.5)

for $f(x)=1$.

We use MscaleDNNs and SD2NNs with aforementioned setups to solve (4.2) when
$\epsilon=0.1$ and $\epsilon=0.01$, respectively. For comparison, a normal DNN
model with _tanh_ activation function (denoted by DNN) is also employed to
solve this multi-scale problem. In SD2NN models, the balance parameter
$\alpha$ for fine-part is set as 0.01. The number of parameters are comparable
for different models (see Table 4 in Appendix A). At each training step, we
randomly sample 3000 interior points and 500 boundary points to evaluate the
loss function. All models are trained for 60000 epochs. In the testing step,
we uniformly sample 1000 points in [0, 1]. We show the testing results for
$\epsilon=0.1$ and $\epsilon=0.01$ in Figures 4 and 5, respectively.

Table 2:  The relative error of different models for Example 1when $\alpha=0.01$. | DNN | Mscale | WWP | SD2NN1 | SD2NN2 | SD2NN3
---|---|---|---|---|---|---
$\epsilon=0.1$ | 2.40e-2 | 3.42e-6 | 3.84e-6 | 1.36e-5 | 6.22e-7 | 7.81e-7
$\epsilon=0.01$ | 2.05e-2 | 6.95e-3 | 1.94e-2 | 2.80e-4 | 7.69e-7 | 3.30e-5

(a) difference of exact and DNN solutions

(b) relative error of DNN, Mscale and WWP

(c) relative error of WWP and SD2NNs

(d) difference of coarse solutions

(e) difference of fine solutions

Figure 4: Testing results with $\epsilon=0.1$ for Example 1.

(a) difference of exact and DNN solutions

(b) relative error of DNN, Mscale and WWP

(c) relative error of WWP and SD2NNs

(d) difference of coarse solutions

(e) difference of fine solutions

Figure 5: Testing results when $\epsilon=0.01$ for Example 1.

From Figures 4(a) and 5(a), we observe that all three SD2NN models capture the
exact solutions with $\epsilon=0.1$ and $\epsilon=0.01$ pretty well. By
comparison, the performances of Mscale and WWP are okay for the case
$\epsilon=0.1$, but deteriorate for the more oscillatory case $\epsilon=0.01$.
The normal DNN model fails in both cases. Figures 4(b), 4(c), 5(b) and 5(c)
and the relative errors in Table 2 further demonstrate that SD2NN2 is the best
method. For the $\epsilon=0.01$ case, the error of SD2NN2 is smaller than
other methods by at least two orders of magnitude.

We choose the coarse part of the exact solution $u(x)$ as $u^{c}(x):=x-x^{2}$,
and the fine part $u^{f}(x)$ as the remainder. Though it is not unique to
choose the fine (and coarse) parts, they are “equivalent” if the difference is
smooth. We draw the differences of coarse solutions of SD2NN1, SD2NN2 and
SD2NN3 with $u^{c}(x)$ in Figures 4(d) and 5(d), and the differences of fine
solutions with $u^{f}(x)$ in Figures 4(e) and 5(e). The differences are
smooth, which shows that SD2NN models can capture the correct coarse and fine
components of the solution.

_Influence of hyper-parameter $\alpha$:_ we study the influence of $\alpha$
for SD2NN model. In the test, we set $\alpha=0.05$, $\epsilon=0.01$, keep all
other parameters fixed, then train all models for 60000 epochs. Based on the
results in Figure 6, we observe that SD2NN2 is more stable and accurate
compared to all other models.

(a) difference of coarse solutions

(b) difference of fine solutions

(c) difference of exact and DNN solutions

Figure 6: Testing results when $\epsilon=0.01$ for Example 1.

_Influence of relaxation-parameter $s$:_ we study the influence of relaxation-
parameter $s$ for activation function _SFM_ in _SD 2NN_ model, especially for
_SD 2NN2_. Here, we list the setup of _SD 2NN2_ in Table 1, as well as two
alternative setups:

  * 1.

_SD 2NN2_: SFM(1.0)+tanh for coarse submodule and SFM(0.5)+S2ReLU for fine
submodule;

  * 2.

_SD 2NN2_(a): SFM(0.5)+tanh for coarse submodule and SFM(0.5)+S2ReLU for fine
submodule;

  * 3.

_SD 2NN2_(b): SFM(1.0)+tanh for coarse submodule and SFM(1.0)+S2ReLU for fine
submodule.

we set $\alpha=0.05$, $\epsilon=0.01$, and train both models for 60000 epochs.
From the results in Figure 7, it seems for the fine submodule of SD2NN2, a
relaxed parameter $s=0.5$ offers better performance.

Figure 7: Relative error for SD2NN2 when $\epsilon=0.01$ for Example 1.

Nonlinear Case:

We further consider the nonlinear case with $p=8$ for (4.2) with
$\Omega=[0,1]$ and
$A(x)=\left(2+\cos\left(2\pi\frac{x}{\epsilon}\right)\right)^{-1}$. By
choosing appropriate $f(x)$, the (unique) solution is given by

$u(x)=x-x^{2}+\epsilon\left(\frac{1}{4\pi}\sin\left(2\pi\frac{x}{\epsilon}\right)-\frac{1}{2\pi}x\sin\left(2\pi\frac{x}{\epsilon}\right)-\frac{\epsilon}{4\pi^{2}}\cos\left(2\pi\frac{x}{\epsilon}\right)+\frac{\epsilon}{4\pi^{2}}\right).$
(4.6)

with $u_{\epsilon}(0)=u_{\epsilon}(1)=0$.

(a) difference of exact and DNN solutions

(b) relative error for different models

(c) difference of coarse solution

(d) difference of fine solution

Figure 8: Testing results for $\epsilon=0.01$ and $p=8$.

We employ Mscale, WWP and SD2NN2 models to solve (4.2). Based on the results
in Figure 8, the SD2NN2 model still outperform Mscale, WWP, SD2NN1 and SD2NN3
(fails to converge) for nonlinear multi-scale problem. In addition, the
differences of coarse and fine solutions in Figures 8(c) and 8(d) show that
the submodules of SD2NN2 also can well capture the coarse and fine parts of
the multi-scale solution.

######  Example 2 (Three-scale problem).

We solve the problem (4.2) with the following three-scale coefficient

$A(x)=\left(2+\cos\left(2\pi\frac{x}{\epsilon_{1}}\right)\right)\left(2+\cos\left(2\pi\frac{x}{\epsilon_{2}}\right)\right)$
(4.7)

with two small parameter $1\gg\epsilon_{1}\gg\epsilon_{2}>0$, also for $p=2$
and $\Omega=[0,1]$. We impose the exact solution

$u(x)=x-x^{2}+\frac{\epsilon_{1}}{4\pi}\sin\left(2\pi\frac{x}{\epsilon_{1}}\right)+\frac{\epsilon_{2}}{4\pi}\sin\left(2\pi\frac{x}{\epsilon_{2}}\right).$
(4.8)

with $u(0)=u(1)=0$, such that $f(x)$ can be obtained by direct computation.

We employ the aforementioned Mscale, WWP, SD2NN1, SD2NN2 and SD2NN3 to solve
(4.2) with (4.7) with $\epsilon_{1}=0.1$ and $\epsilon_{2}=0.01$. We use 2
MscaleDNN submodules in the SD2NN models, and the balance parameters
$\alpha_{1}=0.1$ and $\alpha_{2}=0.01$ for the mesoscale submodule and the
fine submodule, respectively. Table 5 in Appendix A shows that the number of
parameters for different models are comparable. All models are trained for
60000 epochs. For each training step, we randomly sample 3000 interior points
and 500 boundary points as the training data, and uniformly sample 1000 points
in [0,1] as the testing data.

(a) difference for exact solution and approximate solutions

(b) coarse part in the frequency domain

(c) mesoscale part in the frequency domain

(d) finescale part in the frequnecy domain

Figure 9: Numerical results for Example 2 Table 3:  The relative error of different models for Example 2 Mscale | WWP | SD2NN1 | SD2NN2 | SD2NN3
---|---|---|---|---
2.70e-3 | 1.53e-3 | 1.25e-3 | 6.75e-5 | 1.55e-4

In Figure 9(a), we draw the difference of the approximate solution and the
numerical solution, together with the relative errors from Table 3, it shows
that SD2NN2 is the best method. In Figures 9(b) – 9(d) we compare the coarse,
mesoscale, fine scale parts of the SD2NN2 solution with the exact solution, in
the frequency domain. We demonstrate that three submodules of the SD2NN2 model
can correctly capture the corresponding components of the solution as they are
supposed to be.

From now on, we will only consider the Mscale, WWP and SD2NN2 models in the
following examples.

######  Example 3.

We consider the following two-dimensional problem (4.2) for $p=2$ and
$\Omega=[-1,1]\times[-1,1]$. In this example, we set $f=1$ and the multi-scale
coefficient

$A(x_{1},x_{2})=\Pi_{i=1}^{6}\bigg{(}1+0.5\cos\left(2^{i}\pi(x_{1}+x_{2})\right)\bigg{)}\bigg{(}1+0.5\sin\left(2^{i}\pi(x_{2}-3x_{1})\right)\bigg{)}.$

from [71, 72, 73]. The reference solution $u(x_{1},x_{2})$ can be computed by
the finite element method on a square grid of mesh-size $h=1/129$.

We compute the solution of (4.2) by Mscale, WWP and SD2NN2, respectively. In
the SD2NN2 model, the balance parameter $\alpha$ for fine-submodule is 0.05.
The network size and parameter number for different models are comparable and
listed in Table 6 in Appendix A. All models are trained for 100000 epochs. At
each training step, we sample 3000 interior points and 500 boundary points.
The testing points are taken from the finite element grid with $h=1/129$.

(a) point-wise error for Mscale

(b) point-wise error for WWP

(c) point-wise error for SD2NN2

(d) relative error for Mscale, WWP and SD2NN2

(e) SD2NN2 coarse solution

(f) SD2NN2 fine solution

Figure 10: Testing results for Example 3

The numerical results in Figures 10(a) – 10(d) show that the SD2NN2 is still
the method of choice in terms of both the point-wise error and the mean square
error, and it is much better than Mscale and WWP models. In addition, the
results in Figures 10(e) and 10(f) show that the submodules of SD2NN2
successfully separate the coarse and fine parts of the solution.

######  Example 4.

We now consider the following Poisson-Boltzmann equation with Dirichlet
boundary condition,

$\begin{cases}-\textbf{div}(A(\bm{x})\nabla
u(\bm{x}))+\kappa(\bm{x})u(\bm{x})=f(\bm{x}),~{}~{}\bm{x}\in\Omega\subset\mathbb{R}^{d},\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}u(\bm{x})=g(\bm{x}),~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\bm{x}\in\partial\Omega\end{cases}$
(4.9)

where $A_{\epsilon}(\bm{x})$ is the dielectric constant and $\kappa(\bm{x})$
the inverse Debye-Huckel length of an ionic solvent. It is natural to
introduce the energy

$\mathcal{J}(v)=\frac{1}{2}\int_{\Omega}\bigg{(}A|\nabla v|^{2}+\kappa
v^{2}\bigg{)}d\bm{x}-\int_{\Omega}fvd\bm{x}.$ (4.10)

and the corresponding variational formulation.

We solve the elliptic equation (4.9) in the cube $\Omega=[0,1]^{3}$ with 8 big
holes (blue) and 27 smaller holes (red), see Figure 11(a). We take
$\kappa(x_{1}.x_{2},x_{3})=\pi^{2}$,

$A(x_{1},x_{2},x_{3})=0.5\bigg{(}2+\cos(10\pi x_{1})\cos(20\pi
x_{2})\cos(30\pi x_{3})\bigg{)},$ (4.11)

and impose an exact solution

$u(x_{1},x_{2},x_{3})=\sin(\pi x_{1})\sin(\pi x_{2})\sin(\pi
x_{3})+0.05\sin(10\pi x_{1})\sin(20\pi x_{2})\sin(30\pi x_{3}).$ (4.12)

The exact solution prescribes the corresponding boundary condition and $f(x)$
in (4.9).

In this example, the network sizes and the number of parameter for Mscale, WWP
and SD2NN2 are listed in Table 7 in Appendix A, their parameters are
comparable. In addition, the balance parameter $\alpha=0.05$ for the fine-part
of SD2NN2. All models are trained for 100000 epochs. At each training step,
the training data set including 6000 interior points and 1000 boundary points
randomly sampled from $\Omega$ and $\partial\Omega$. The testing dataset is
1600 random samples in $\Omega$. Testing results are plotted in Fig.11.

(a) Domain with holes

(b) point-wise error for Mscale

(c) point-wise error for WWP

(d) point-wise error for SD2NN2

(e) relative error of Mscale, WWP and SD2NN2

Figure 11: Testing results for Example 4.

From the results in 11(b) – 11(d), the SD2NN2 model still keeps its good
performance for the Poisson-Boltzmann equation in 3D perforated domain.

##  5 Conclusion

Deep learning algorithms have demonstrated great potential in scientific
computing tasks. However, the efficient solution of multi-scale problems
remains to be a big challenge for DNN based numerical methods. In this paper,
we combine subspace decomposition ideas from traditional numerical analysis
and multi-scale deep neural network, and propose the _SD 2NN_ method for
multi-scale PDEs. This new architecture consists of one low-frequency or
normal DNN submodule and one (or several) high-frequency MscaleDNN
submodule(s). In addition, we incorporate the trigonometric _SFM_ activation
function for the SD2NN model. Computational results show that this new method
is feasible and efficient for multi-scale problems in 1d, 2d and 3d, and in
both regular or perforated domains. In the future, we plan to apply this
method to multi-scale PDEs with nonlinearity and/or randomness, as well as
operator learning for multi-scale PDEs.

## Acknowledgements

X.A.L and L.Z are partially supported by the National Natural Science
Foundation of China (NSFC 11871339, 11861131004). Z.X. is supported by the
National Key R&D Program of China Grant No. 2019YFA0709503, the Shanghai
Sailing Program, the Natural Science Foundation of Shanghai Grant No.
20ZR1429000, the National Natural Science Foundation of China Grant No.
62002221, Shanghai Municipal of Science and Technology Project Grant No.
20JC1419500, Shanghai Municipal of Science and Technology Major Project NO.
2021SHZDZX0102, and the HPC of School of Mathematical Sciences and the Student
Innovation Center at Shanghai Jiao Tong University. Thanks to Xin-Liang Liu
for his some codes and helpful suggestions.

## References

  * LeCun et al. [2015] Y. LeCun, Y. Bengio, G. Hinton,  Deep learning,  Nature 521 (2015) 436–444.
  * Goodfellow et al. [2016] I. Goodfellow, Y. Bengio, A. Courville, Deep Learning, MIT press, Cambridge, 2016.
  * E [2017] W. E,  A Proposal on Machine Learning via Dynamical Systems,  Communications in Mathematics and Statistics 5 (2017) 1–11.
  * E et al. [2017] W. E, J. Han, A. Jentzen,  Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations,  Communications in Mathematics and Statistics 5 (2017) 349–380.
  * Han et al. [2018] J. Han, L. Zhang, R. Car, et al.,  Deep Potential: A General Representation of a Many-Body Potential Energy Surface,  Communications in Computational Physics 23 (2018).
  * Raissi et al. [2019] M. Raissi, P. Perdikaris, G. E. Karniadakis,  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations,  Journal of Computational Physics 378 (2019) 686–707.
  * Weinan [2020] E. Weinan,  Machine learning and computational mathematics,  arXiv preprint arXiv:2009.14596 (2020).
  * Shi and Xu [2021] E. Shi, C. Xu,  A comparative investigation of neural networks in solving differential equations,  Journal of Algorithms & Computational Technology 15 (2021) 1–15.
  * Owhadi [2003] H. Owhadi,  Anomalous slow diffusion from perpetual homogenization,  Ann. Probab. 31 (2003) 1935–1969.
  * Homeyer and Gohlke [2012] N. Homeyer, H. Gohlke,  Free Energy Calculations by the Molecular Mechanics Poisson-Boltzmann Surface Area Method.,  Molecular Informatics 31 (2012) 114–122.
  * Hagelaar and Pitchford [2005] G. J. M. Hagelaar, L. C. Pitchford,  Solving the Boltzmann equation to obtain electron transport coefficients and rate coefficients for fluid models,  Plasma Sources Science and Technology 14 (2005) 722–733.
  * Cohen et al. [2012] A. Cohen, W. Dahmen, G. Welper,  Adaptivity and variational stabilization for convection-diffusion equations,  Mathematical Modelling and Numerical Analysis 46 (2012) 1247–1273.
  * Papanicolau et al. [1978] G. Papanicolau, A. Bensoussan, J.-L. Lions, Asymptotic analysis for periodic structures, Elsevier, 1978.
  * Jikov et al. [2012] V. Jikov, S. M. Kozlov, O. A. Oleinik, Homogenization of differential operators and integral functionals, Springer Science & Business Media, 2012\.
  * E and Engquist [2003] W. E, B. Engquist,  Heterogeneous multiscale methods,  Comm. Math. Sci. 1 (2003) 87–132.
  * Ming et al. [2005] P. Ming, P. Zhang, et al.,  Analysis of the heterogeneous multiscale method for elliptic homogenization problems,  Journal of the American Mathematical Society 18 (2005) 121–156.
  * Berlyand et al. [2013] L. Berlyand, A. G. Kolpakov, A. Novikov, Introduction to the network approximation method for materials modeling, volume 148, Cambridge University Press, 2013.
  * Hou and Wu [1997] T. Y. Hou, X.-H. Wu,  A multiscale finite element method for elliptic problems in composite materials and porous media,  Journal of Computational Physics 134 (1997) 169–189.
  * Hughes et al. [1998] T. J. R. Hughes, G. Feijoo, L. Mazzei, J. Quincy,  The variational multiscale method—a paradigm for computational mechanics,  Comput. Methods Appl. Mech. Engrg. 166 (1998) 3–24.
  * Bazilevs et al. [2007] Y. Bazilevs, V. Calo, J. Cottrell, T. Hughes, A. Reali, G. Scovazzi,  Variational multiscale residual-based turbulence modeling for large eddy simulation of incompressible flows,  Computer Methods in Applied Mechanics and Engineering 197 (2007) 173–201.
  * Berlyand and Owhadi [2010] L. Berlyand, H. Owhadi,  Flux norm approach to finite dimensional homogenization approximations with non-separated scales and high contrast,  Archive for rational mechanics and analysis 198 (2010) 677–721.
  * Owhadi and Zhang [2008] H. Owhadi, L. Zhang,  Homogenization of parabolic equations with a continuum of space and time scales,  SIAM Journal on Numerical Analysis 46 (2008) 1–36.
  * Owhadi et al. [2014] H. Owhadi, L. Zhang, L. Berlyand,  Polyharmonic homogenization, rough polyharmonic splines and sparse super-locatization,  ESIAM 48 (2014) 517–552.
  * Liu et al. [2021] X. Liu, L. Zhang, S. Zhu,  Generalized rough polyharmonic splines for multiscale pdes with rough coefficients,  Numerical Mathematics: Theory, Methods and Applications 14 (2021) 862–892.
  * Efendiev et al. [2013] Y. Efendiev, J. Galvis, T. Hou,  Generalized multiscale finite element methods,  Journal of Computational Physics 251 (2013) 116–135.
  * Chung et al. [2018] E. T. Chung, Y. Efendiev, W. T. Leung,  An adaptive generalized multiscale discontinuous galerkin method for high-contrast flow problems,  Multiscale Modeling & Simulation 16 (2018) 1227–1257.
  * Målqvist and Peterseim [2014] A. Målqvist, D. Peterseim,  Localization of elliptic multiscale problems,  Math. Comp. 83 (2014) 2583–2603.
  * Henning et al. [2014] P. Henning, A. Målqvist, D. Peterseim,  A localized orthogonal decomposition method for semi-linear elliptic problems,  ESAIM Math. Model. Numer. Anal. 48 (2014) 1331–1349.
  * Xu [1992] J. Xu,  Iterative methods by space decomposition and subspace correction,  SIAM review 34 (1992) 581–613.
  * Xie et al. [2019] H. Xie, L. Zhang, H. Owhadi,  Fast eigenpairs computation with operator adapted wavelets and hierarchical subspace correction,  SIAM Journal on Numerical Analysis 57 (2019) 2519–2550.
  * Owhadi [2015] H. Owhadi,  Bayesian numerical homogenization,  Multiscale Model. Simul. 13 (2015) 812–828.
  * Owhadi [2017] H. Owhadi,  Multigrid with rough coefficients and Multiresolution operator decomposition from Hierarchical Information Games,  SIAM Rev. 59 (2017) 99–149.
  * Owhadi et al. [2019] H. Owhadi, C. Scovel, G. R. Yoo,  Kernel mode decomposition and programmable/interpretable regression networks,  arXiv preprint arXiv:1907.08592 (2019).
  * Owhadi and Yoo [2019] H. Owhadi, G. R. Yoo,  Kernel flows: From learning kernels from data into the abyss,  Journal of Computational Physics 389 (2019) 22–47.
  * Sirignano and Spiliopoulos [2018] J. Sirignano, K. Spiliopoulos,  DGM: A deep learning algorithm for solving partial differential equations,  Journal of Computational Physics 375 (2018) 1339–1364.
  * Zang et al. [2020] Y. Zang, G. Bao, X. Ye, H. Zhou,  Weak adversarial networks for high-dimensional partial differential equations,  Journal of Computational Physics 411 (2020) 109409.
  * Xu et al. [2019] Z.-Q. J. Xu, Y. Zhang, Y. Xiao,  Training behavior of deep neural network in frequency domain,  International Conference on Neural Information Processing (2019) 264–274.
  * Xu et al. [2020] Z.-Q. J. Xu, Y. Zhang, T. Luo, Y. Xiao, Z. Ma,  Frequency Principle: Fourier Analysis Sheds Light on Deep Neural Networks,  Communications in Computational Physics 28 (2020) 1746–1767.
  * Rahaman et al. [2019] N. Rahaman, D. Arpit, A. Baratin, F. Draxler, M. Lin, F. A. Hamprecht, Y. Bengio, A. Courville,  On the spectral bias of deep neural networks,  International Conference on Machine Learning (2019).
  * Zhang et al. [2021] Y. Zhang, T. Luo, Z. Ma, Z.-Q. J. Xu,  A linear frequency principle model to understand the absence of overfitting in neural networks,  Chinese Physics Letters 38 (2021) 038701.
  * Luo et al. [2019] T. Luo, Z. Ma, Z.-Q. J. Xu, Y. Zhang,  Theory of the Frequency Principle for General Deep Neural Networks,  arXiv preprint arXiv:1906.09235 (2019).
  * Ronen et al. [2019] B. Ronen, D. Jacobs, Y. Kasten, S. Kritchman,  The convergence rate of neural networks for learned functions of different frequencies,  in: Advances in Neural Information Processing Systems, volume 32, 2019, pp. 4761–4771.
  * E et al. [2019] W. E, C. Ma, L. Wu,  Machine learning from a continuous viewpoint,  arXiv preprint arXiv:1912.12777 (2019).
  * Cao et al. [2019] Y. Cao, Z. Fang, Y. Wu, D.-X. Zhou, Q. Gu,  Towards understanding the spectral bias of deep learning,  arXiv preprint arXiv:1912.01198 (2019).
  * Yang and Salman [2019] G. Yang, H. Salman,  A fine-grained spectral perspective on neural networks,  arXiv preprint arXiv:1907.10599 (2019).
  * Bordelon et al. [2020] B. Bordelon, A. Canatar, C. Pehlevan,  Spectrum dependent learning curves in kernel regression and wide neural networks,  in: Proceedings of the 37th International Conference on Machine Learning, ICML 2020, 13-18 July 2020, Virtual Event, volume 119 of Proceedings of Machine Learning Research, PMLR, 2020, pp. 1024–1034.
  * Luo et al. [2020] T. Luo, Z. Ma, Z.-Q. J. Xu, Y. Zhang,  On the exact computation of linear frequency principle dynamics and its generalization,  arXiv preprint arXiv:2010.08153 (2020).
  * Ma et al. [2020] C. Ma, L. Wu, W. E,  The slow deterioration of the generalization error of the random feature model,  in: Mathematical and Scientific Machine Learning, PMLR, 2020, pp. 373–389.
  * Sharma and Ross [2020] R. Sharma, A. Ross,  D-netpad: An explainable and interpretable iris presentation attack detector,  in: 2020 IEEE International Joint Conference on Biometrics (IJCB), IEEE, 2020, pp. 1–10.
  * Zhu et al. [2019] H. Zhu, Y. Qiao, G. Xu, L. Deng, Y. Yu-Feng,  Dspnet: A lightweight dilated convolution neural networks for spectral deconvolution with self-paced learning,  IEEE Transactions on Industrial Informatics (2019).
  * Chakrabarty and Maji [2019] P. Chakrabarty, S. Maji,  The spectral bias of the deep image prior,  arXiv preprint arXiv:1912.08905 (2019).
  * Xu and Zhou [2021] Z.-Q. J. Xu, H. Zhou,  Deep frequency principle towards understanding why deeper learning is faster,  in: Proceedings of the AAAI Conference on Artificial Intelligence, volume 35, 2021.
  * Wang et al. [2021] S. Wang, H. Wang, P. Perdikaris,  On the eigenvector bias of fourier feature networks: From regression to solving multi-scale pdes with physics-informed neural networks,  Computer Methods in Applied Mechanics and Engineering 384 (2021) 113938\.
  * Cai et al. [2020] W. Cai, X. Li, L. Liu,  A phase shift deep neural network for high frequency approximation and wave problems,  SIAM Journal on Scientific Computing 42 (2020) A3285–A3312.
  * Jagtap et al. [2020] A. D. Jagtap, K. Kawaguchi, G. E. Karniadakis,  Adaptive activation functions accelerate convergence in deep and physics-informed neural networks,  Journal of Computational Physics 404 (2020) 109136.
  * Liu et al. [2020] Z. Liu, W. Cai, Z.-Q. J. Xu,  Multi-scale deep neural network (mscalednn) for solving poisson-boltzmann equation in complex domains,  Communications in Computational Physics 28 (2020) 1970–2001.
  * Li et al. [2020] X.-A. Li, Z.-Q. J. Xu, L. Zhang,  A multi-scale dnn algorithm for nonlinear elliptic equations with multiple scales,  Communications in Computational Physics 28 (2020) 1886–1906.
  * Wang et al. [2020] B. Wang, W. Zhang, W. Cai,  Multi-scale deep neural network (mscalednn) methods for oscillatory stokes flows in complex domains,  Communications in Computational Physics 28 (2020) 2139–2157.
  * Tancik et al. [2020] M. Tancik, P. P. Srinivasan, B. Mildenhall, S. Fridovich-Keil, N. Raghavan, U. Singhal, R. Ramamoorthi, J. T. Barron, R. Ng,  Fourier features let networks learn high frequency functions in low dimensional domains,  arXiv preprint arXiv:2006.10739 (2020).
  * He et al. [2016] K. He, X. Zhang, S. Ren, J. Sun,  Deep residual learning for image recognition,  in: 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2016, pp. 770–778.
  * E and Yu [2018] W. E, B. Yu,  The Deep Ritz Method: A Deep Learning-Based Numerical Algorithm for Solving Variational Problems,  Communications in Mathematics and Statistics 6 (2018) 1–12.
  * Zou et al. [2020] Z. Zou, H. Zhang, Y. Guan, J. Zhang,  Deep residual neural networks resolve quartet molecular phylogenies.,  Molecular Biology and Evolution 37 (2020) 1495–1507.
  * Lions et al. [1978] J. L. Lions, A. Bensoussan, G. Papanicolaou, Asymptotic analysis for periodic structures, volume 5, North Holland, Amsterdam, 1978.
  * Hughes et al. [1998] T. J. R. Hughes, G. R. Feijóo, L. Mazzei, J.-B. Quincy,  The variational multiscale method—a paradigm for computational mechanics,  Computer Methods in Applied Mechanics and engineering 166 (1998) 3–24.
  * Owhadi [2017] H. Owhadi,  Multigrid with rough coefficients and multiresolution operator decomposition from hierarchical information games,  Siam Review 59 (2017) 99–149.
  * Raissi et al. [2017] M. Raissi, P. Perdikaris, G. E. Karniadakis,  Inferring solutions of differential equations using noisy multi-fidelity data,  Journal of Computational Physics 335 (2017) 736–746.
  * Robert and Casella [1999] C. P. Robert, G. Casella, Monte Carlo Statistical Methods, 1999.
  * Quarteroni et al. [2007] A. Quarteroni, R. Sacco, F. Saleri, Numerical Mathematics, Numerical Mathematics, 2007.
  * Glorot and Bengio [2010] X. Glorot, Y. Bengio,  Understanding the difficulty of training deep feedforward neural networks,  in: Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, 2010, pp. 249–256.
  * Rotskoff and Vanden-Eijnden [2018] G. M. Rotskoff, E. Vanden-Eijnden,  Parameters as interacting particles: long time convergence and asymptotic error scaling of neural networks,  in: 32nd Conference on Neural Information Processing Systems, NeurIPS 2018, 2018, pp. 7146–7155.
  * Owhadi and Zhang [2007] H. Owhadi, L. Zhang,  Homogenization of Parabolic Equations with a Continuum of Space and Time Scales,  SIAM Journal on Numerical Analysis 46 (2007) 1–36.
  * Owhadi and Zhang [2008] H. Owhadi, L. Zhang,  Numerical homogenization of the acoustic wave equations with a continuum of scales,  Computer Methods in Applied Mechanics and Engineering 198 (2008) 397–406.
  * Owhadi et al. [2014] H. Owhadi, L. Zhang, L. Berlyand,  Polyharmonic homogenization, rough polyharmonic splines and sparse super-localization,  Mathematical Modelling and Numerical Analysis 48 (2014) 517–552.

## Appendix

## Appendix A

Table 4:  Network sizes and number of parameters for models in Example 1 | Network Size | Number of Parameters
---|---|---
DNN | (250, 100, 80, 80, 60) | 44510
Mscale | (250, 100, 80, 80, 60) | 44510
WWP | (125, 100, 80, 80, 60) | 44385
SD2NN1 | coarse:(100, 80, 60, 60, 40) | 44440
fine:(125, 60, 60, 60, 50)
SD2NN2 | coarse:(50, 80, 60, 60, 40) | 44315
fine:(100, 60, 60, 50, 40)
SD2NN3 | coarse:(50, 80, 60, 60, 40) | 44315
fine:(100, 60, 60, 50, 40)
Table 5:  Network sizes and number of parameters for models in Example 2 | Network Size | Para. | Scaling factors
---|---|---|---
Mscale | (350, 300, 200, 200, 100) | 225450 | ——
WWP | (175, 300, 200, 200, 100) | 225275 | ——
SD2NN1 | coarse:(100, 80, 60, 60, 40) | 207730 | ——
mesoscale:(125, 80, 60, 60, 40) | $(30,31,32,\cdots,69,70)$
fine:(225, 200, 150, 150, 100) | $(251,252,253,\cdots,360)$
SD2NN2 | coarse:(50, 80, 60, 60, 40) | 207680 | $(0.5,1,1.5,\cdots,24.5,25)$
mesoscale:(125, 80, 60, 60, 40) | $(30,31,32,\cdots,69,70)$
fine:(225, 200, 150, 150, 100) | $(251,252,253,\cdots,360)$
SD2NN3 | coarse:(50, 80, 60, 60, 40) | 207680 | $(0.5,1,1.5,\cdots,24.5,25)$
mesoscale:(125, 80, 60, 60, 40) | $(30,31,32,\cdots,69,70)$
fine:(225, 200, 150, 150, 100) | $(251,252,253,\cdots,360)$
Table 6:  Network sizes and parameter of models for Example 3 | Network Size | Para.
---|---|---
Mscale | (250, 200, 200, 100, 100, 80) | 128330
WWP | (125, 200, 200, 100, 100, 80) | 128205
SD2NN2 | coarse:(50, 100, 80, 80, 80, 60) | 127410
fine:(120, 150, 150, 100, 100, 80)
Table 7:  Network sizes and parameter of different models for Example 4 | Network Size | Para.
---|---|---
Mscale | (500, 400, 400, 200, 200, 150) | 510650
WWP | (250, 400, 400, 200, 200, 150) | 510400
SD2NN2 | coarse:(70, 200, 200, 150, 150, 150) | 508870
fine:(250, 300, 290, 200, 200, 150)


* [1] M. F. Adamer, M. Helmer, Complexity of model testing for dynamical systems with toric steady states, Adv. in Appl. Math. 110 (2019), 42–75.
  * [2] C. Améndola, L. Gustafsson, K. Kohn, O. Marigliano, A. Seigal, The maximum likelihood degree of linear spaces of symmetric matrices, Matematiche (Catania), 76 (2021), no. 2, 535–557.
  * [3] C. Améndola, P. Zwiernik, Likelihood geometry of correlation models, Matematiche (Catania), 76 (2021), no. 2, 559–583.
  * [4] S. Agarwal, N. Snavely, I. Simon, S. Seitz, R. Szeliski, Building Rome in a day, Communications of the ACM, 54 (2011), 105–112.
  * [5] E. L. Allgower, K. Georg, Numerical continuation methods. An introduction. Springer Series in Computational Mathematics, 13. Springer-Verlag, Berlin, 1990.
  * [6] P. Aluffi, Differential forms with logarithmic poles and Chern-Schwartz-MacPherson classes of singular varieties, C. R. Acad. Sci. Paris Sér. I Math. 329 (1999), 619–624.
  * [7] P. Aluffi, Limits of Chow groups, and a new construction of Chern-Schwartz-MacPherson classes, Pure Appl. Math. Q. 2 (2006), no. 4, Special Issue: In honor of Robert D. MacPherson. Part 2, 915–941.
  * [8] P. Aluffi, Euler characteristics of general linear sections and polynomial Chern classes, Rend. Circ. Mat. Palermo 62 (2013), 3–26.
  * [9] P. Aluffi, Projective duality and a Chern-Mather involution, Transactions of the American Mathematical Society, 370 (2018), 1803–1822.
  * [10] P. Aluffi, C. Harris, The Euclidean distance degree of smooth complex projective varieties, Algebra Number Theory 12 (2018), no. 8, 2005–2032.
  * [11] P. Aluffi, L. Mihalcea, J. Schürmann, C. Su, Shadows of characteristic cycles, Verma modules, and positivity of Chern-Schwartz-MacPherson classes of Schubert cells, arXiv:1709.08697.
  * [12] B. D. O. Anderson, U. Helmke, Counting critical formations on a line, SIAM J. Control Optim. 52 (2014), 219–242.
  * [13] T. Boege, J. I. Coons, C. Eur, A. Maraj, F. Röttger, Reciprocal maximum likelihood degrees of Brownian motion tree models, Matematiche (Catania), 76 (2021), no. 2, 383–398.
  * [14] J.-P. Brasselet, Characteristic classes, in Handbook of geometry and topology of singularities III, 303–417, Springer, Cham, 2022.
  * [15] J.-P. Brasselet, D. Massey, A. J. Parameswaran, J. Seade, Euler obstruction and defects of functions on singular varieties, J. London Math. Soc. (2), 70 (2004), 59–76.
  * [16] P. Breiding, F. Rydell, E. Shehu, A. Torres, Line Multiview Varieties, arXiv:2203.01694.
  * [17] E. Brieskorn, Die Monodromie der isolierten Singularitäten von Hyperflächen, Manuscripta Math. 2 (1970) 103–161.
  * [18] F. Catanese, S. Hoşten, A. Khetan, B. Sturmfels, The maximum likelihood degree, Amer. J. Math. 128 (2006), no. 3, 671–697.
  * [19] T. O. Çelik, A. Jamneshan, G. Montúfar, B. Sturmfels, L. Venturello, Wasserstein distance to independence models, J. Symbolic Comput., 104 (2021), 855–873.
  * [20] M. Compagnoni, R. Notari, F. Antonacci, A. Sarti, A comprehensive analysis of the geometry of TDOA maps in localization problems, Inverse Problems 30 (2014), 035004 (49pp).
  * [21] J. Draisma, E. Horobeţ, G. Ottaviani, B. Sturmfels, R. Thomas, The Euclidean distance degree of an algebraic variety, Found. Comput. Math. 16 (2016), no. 1, 99–149.
  * [22] S. Di Rocco, L. Gustafsson, L. Schaffler, Gaussian likelihood geometry of projective varieties, arXiv:2208.12560.
  * [23] S. Di Rocco, D. Eklund, M. Weinstein, The bottleneck degree of algebraic varieties, SIAM J. Appl. Algebra Geom, 4 (2020), no. 1, 227–253.
  * [24] A. Dimca, Sheaves in topology, Universitext. Springer-Verlag, Berlin, 2004.
  * [25] M. Drton, B. Sturmfels, S. Sullivant, Lectures on Algebraic Statistics, Oberwolfach Seminars, Vol 39, Birkhäuser, Basel, 2009.
  * [26] S. Dye, K. Kohn, F. Rydell, R. Sinn, Maximum likelihood estimation for nets of conics, Matematiche (Catania), 76 (2021), no. 2, 399–414.
  * [27] C. Eur, T. Fife, J. A. Samper, T. Seynnaeve, Reciprocal maximum likelihood degrees of diagonal linear concentration models, Matematiche (Catania), 76 (2021), no. 2, 447–459.
  * [28] G. Ewald, Combinatorial convexity and algebraic geometry, Graduate Texts in Mathematics, 168. Springer-Verlag, New York, 1996.
  * [29] T. Gaffney, N. Grulha, M. Ruas, The local Euler obstruction and topology of the stabilization of associated determinantal varieties, Math. Z. 291 (2019), no. 3-4, 905–930.
  * [30] V. Ginsburg, Characteristic varieties and vanishing cycles, Invent. Math. 84 (1986), 327–402.
  * [31] E. Gross, H. A. Harrington, Z. Rosen, B. Sturmfels, Algebraic systems biology: a case study for the Wnt pathway, Bull. Math. Biol., 78 (2016), no.1, 21–51.
  * [32] D. Grayson, M. Stillman, Macaulay 2, available at http://www.math.uiuc.edu/Macaulay2.
  * [33] C. Harris, D. Lowengrub, The Chern-Mather class of the multiview variety, Comm. Algebra 46 (2018), no. 6, 2488–2499.
  * [34] D. L. Hartl, A. G. Clark, Principles of population genetics, Sinauer Associates, Inc. Publishers, Sunderland, Massachusetts (Fourth Edition).
  * [35] R. Hartley, A. Zisserman, Multiple view geometry in computer vision, Second edition. With a foreword by Olivier Faugeras. Cambridge University Press, Cambridge, 2003.
  * [36] J. D. Hauenstein, Numerically computing real points on algebraic sets, Acta Appl. Math., 125 (2013), 105–119.
  * [37] D. Hauenstein, J. I. Rodriguez, B. Sturmfels, Maximum likelihood for matrices with rank constraints, J. Algebr. Stat. 5 (2014), 18–38.
  * [38] M. Helmer, V. Nanda, Conormal Spaces and Whitney Stratifications, Foundations of Computational Mathematics (2022).
  * [39] M. Helmer, B. Sturmfels, Nearest points on toric varieties, Math. Scand. 122 (2018), no. 2, 213–238.
  * [40] E. Horobeţ, The data singular and the data isotropic loci for affine cones, Comm. Algebra, 45 (2017), no. 3, 1177–1186.
  * [41] E. Horobeţ, J. I. Rodriguez, Data loci in algebraic optimization, J. Pure Appl. Algebra, 226 (2022), no. 12, Paper No. 107144, 15 pp.
  * [42] E. Horobeţ, M. Weinstein, Offset hypersurfaces and persistent homology of algebraic varieties, Comput. Aided Geom. Design 74 (2019), 101767, 14 pp.
  * [43] J. Huh, The maximum likelihood degree of a very affine variety, Compos. Math. 149 (2013), no. 8, 1245–1266.
  * [44] J. Huh, B. Sturmfels, Likelihood geometry, in Combinatorial algebraic geometry, 63–117, Lecture Notes in Math., 2108, Fond. CIME/CIME Found. Subser., Springer, Cham, 2014.
  * [45] Y. Jiang, K. Kohn, R. Winter, Linear spaces of symmetric matrices with non-maximal maximum likelihood degree, Matematiche (Catania), 76 (2021), no. 2, 461–481.
  * [46] K. Kubjas, O. Kuznetsova, L. Sodomaco, Algebraic degree of optimization over a variety with an application to p-norm distance degree, arXiv:2105.07785.
  * [47] N. Kowsuwan, P. Kanongchaiyos, 3d cloud animation using ca based method, in 2009 International Symposium on Intelligent Signal Processing and Communication Systems (ISPACS), pages 387–392. IEEE, 2009.
  * [48] J. Lindberg, N. Nicholson, J.I. Rodriguez, Z. Wang, The Maximum Likelihood Degree of Sparse Polynomial Systems, SIAM J. Appl. Algebra Geometry 7 (2023), no. 1, 159–171.
  * [49] R. MacPherson, Chern classes for singular algebraic varieties, Ann. of Math. 100 (1974), 423–432.
  * [50] J. Milnor, Singular points of complex hypersurfaces, Annals of Mathematics Studies, No. 61, Princeton University Press, Princeton, N.J.; University of Tokyo Press, Tokyo, 1968.
  * [51] A. Martín del Campo, J. I. Rodriguez, Critical points via monodromy and local methods, J. Symbolic Comput., 79 (2017), 559–574.
  * [52] J. Mather, Notes on topological stability, Bull. Amer. Math. Soc. (N.S.), 49 (2012), 475–506.
  * [53] A. Milella, G. Reina, 3d reconstruction and classification of natural environments by an autonomous vehicle using multi-baseline stereo, Intelligent Service Robotics, 7 (2014), 79–92.
  * [54] L. Maxim, Intersection Homology & Perverse Sheaves, with Applications to Singularities, Graduate Texts in Mathematics, Vol. 281, Springer, 2019.
  * [55] L. Maxim, J. I. Rodriguez, B. Wang, Euclidean distance degree of the multiview variety, SIAM J. Appl. Algebra Geometry 4 (2020), no. 1, 28–48.
  * [56] L. Maxim, J. I. Rodriguez, B. Wang, Euclidean distance degree of projective varieties, Int. Math. Res. Not. IMRN 2021, no. 20, 15788–15802.
  * [57] L. Maxim, J. I. Rodriguez, B. Wang, Defect of Euclidean distance degree, Adv. in Appl. Math. 121 (2020), 102101, 22 pp.
  * [58] L. Maxim, J. I. Rodriguez, B. Wang, A Morse theoretic approach to non-isolated singularities and applications to optimization, J. Pure Appl. Algebra 226 (2022), no. 3, Paper No. 106865, 23 pp.
  * [59] L. Maxim, J. I. Rodriguez, B. Wang, L. Wu, Logarithmic cotangent bundles, Chern-Mather classes, and the Huh-Sturmfels Involution conjecture, Comm. Pure Appl. Math. (accepted); arXiv:2202.00554.
  * [60] L. Maxim, J. I. Rodriguez, B. Wang, L. Wu, Linear optimization on varieties and Chern-Mather classes, arXiv:2208.09073.
  * [61] L. Maxim, J. Schürmann, Constructible sheaf complexes in complex geometry and Applications, in Handbook of Geometry and Topology of Singularities, vol. III, pages 679–791, Springer, Cham, 2022.
  * [62] L. Maxim, M. Tibăr, Euclidean distance degree and limit points in a Morsification, arXiv:2210.06022.
  * [63] T. Ohmoto, An elementary remark on the integral with respect to Euler characteristics of projective hyperplane sections, Rep. Fac. Sci. Kagoshima Univ. No. 36 (2003), 37–41.
  * [64] G. Ottaviani, P.-J. Spaenlehauer, B. Sturmfels, Exact solutions in structured low-rank approximation, SIAM J. Matrix Anal. Appl., 35 (4):1521–1542, 2014.
  * [65] R. S. Palais, S. Smale, A generalized Morse theory, Bull. Amer. Math. Soc. 70 (1964), 165–172.
  * [66] A. Parusinski, P. Pragacz, A formula for the Euler characteristic of singular hypersurfaces, J. Algebraic Geom., 4 (1995), 337–351.
  * [67] J. I. Rodriguez, X. Tang, Data-discriminants of likelihood equations, in ISSAC’15—Proceedings of the 2015 ACM International Symposium on Symbolic and Algebraic Computation, pages 307–314. ACM, New York, 2015.
  * [68] J. I. Rodriguez, X. Tang, A probabilistic algorithm for computing data-discriminants of likelihood equations, J. Symbolic Comput., 83 (2017), :342–364.
  * [69] J. I. Rodriguez, B. Wang, The maximum likelihood degree of mixtures of independence models, SIAM J. Appl. Algebra Geom. 1 (2017), no. 1, 484–506.
  * [70] J. I. Rodriguez, B. Wang, Computing Euler obstruction functions using maximum likelihood degrees, Int. Math. Res. Not. IMRN 2020, no. 20, 6699–6712.
  * [71] P. Rostalski, B. Sturmfels, Dualities, in Semidefinite optimization and convex algebraic geometry, MOS-SIAM Ser. Optim. vol. 13, pag. 203–249.
  * [72] F. Rydell, E. Shehu, A. Torres, Theoretical and Numerical Analysis of 3D Reconstruction Using Point and Line Incidences, arXiv:2303.13593.
  * [73] C. Sabbah, Quelques remarques sur la géométrie des espaces conormaux, in Differential systems and singularities (Luminy, 1983), Astérisque No. 130 (1985), 161–192.
  * [74] J. Schürmann, Topology of singular spaces and constructible sheaves, Instytut Matematyczny Polskiej Akademii Nauk. Monografie Matematyczne (New Series), vol. 63, Birkhäuser Verlag, Basel, 2003.
  * [75] J. Seade, On Milnor’s fibration theorem and its offspring after 50 years, Bull. Amer. Math. Soc. (N.S.) 56 (2019), no. 2, 281–348.
  * [76] J. Seade, M. Tibăr, A. Verjovsky, Global Euler obstruction and polar invariants, Math. Ann. 333 (2005), 393–403.
  * [77] H. Stewenius, F. Schaffalitzky, D. Nister, How hard is 3-view triangulation really?, in Tenth IEEE International Conference on Computer Vision (ICCV’05), volume 1, pages 686–693. IEEE, 2005.
  * [78] B. Sturmfels, Beyond linear algebra, arXiv:2108.09494.
  * [79] B. Sturmfels, C. Uhler, Multivariate Gaussian, semidefinite matrix completion, and convex algebraic geometry, Ann. Inst. Statist. Math. 62 (2010), no. 4, 603–638.
  * [80] Z. Sun, U. Helmke, B. D. O. Anderson, Rigid formation shape control in general dimensions: an invariance principle and open problems, in 2015 54th IEEE Conference on Decision and Control (CDC), 6095–6100, Dec 2015.
  * [81] C. Uhler, Geometry of maximum likelihood estimation in Gaussian graphical models, Ann. Statist. 40 (2012), no. 1, 238–261.
  * [82] S. Sullivant, Algebraic statistics, Graduate Studies in Mathematics, volume 194, American Mathematical Society, Providence, RI, 2018.
  * [83] L. D. Tráng, B. Teissier, Variétés polaires locales et classes de Chern des variétés singulières, Ann. of Math. (2), 114 (1981), no 3, 457–491.
  * [84] H. Whitney, Local properties of analytic varieties, 1965 Differential and Combinatorial Topology (A Symposium in Honor of Marston Morse) pp. 205–244 Princeton Univ. Press, Princeton, N.J.
  * [85] H. Whitney, Tangents to an analytic variety, Ann. of Math. (2) 81 (1965), 496–549.
  * [86] L. Wu, P. Zhou, Log D-modules and index theorems, Forum Math. Sigma 9 (2021), Paper No. e3, 32 pp.