The Finite Neuron Method and Convergence Analysis

Let $\Omega\subset\mathbb{R}^{d}$ be a bounded domain with a sufficiently smooth boundary $\partial\Omega$. For any integer $m\geq 1$, we consider the following model $2m$-th order partial differential equation with certain boundary conditions:

$$\left\{\begin{array}[]{rccl}\displaystyle Lu&=&f&\mbox{in }\Omega,\\ B^{k}(u)&=&0&\mbox{on }\partial\Omega\quad(0\leq k\leq m-1),\end{array}\right.$$

(1.1)

where $L$ is a partial differential operator as follows

$$Lu=\sum_{|\alpha|=m}(-1)^{m}\partial^{\alpha}(a_{\alpha}(x)\,\partial^{\alpha}\,u)+a_{0}(x)u,$$

(1.2)

and ${\bm{\alpha}}$ denotes $n$-dimensional multi-index ${\bm{\alpha}}=(\alpha_{1},\cdots,\alpha_{n})$ with

$$|{\bm{\alpha}}|=\sum_{i=1}^{n}\alpha_{i},\quad\partial^{\bm{\alpha}}=\frac{\partial^{|{\bm{\alpha}}|}}{\partial x_{1}^{\alpha_{1}}\cdots\partial x_{n}^{\alpha_{n}}}.$$

For simplicity, we assume that $a_{\alpha}$ are strictly positive and smooth functions on $\Omega$ for $|\alpha|=m$ and $\alpha=0$, namely, $\exists\alpha_{0}>0$, such that

$$a_{\alpha}(x),a_{0}(x)\geq\alpha_{0},\,\,\forall x\in\Omega,\,\,|\alpha|=m.$$

(1.3)

Given a nonnegative integer $k$ and a bounded domain $\Omega\subset\mathbb{R}^{d}$, let

$$H^{k}(\Omega):=\left\{v\in L^{2}(\Omega),\partial^{\alpha}v\in L^{2}(\Omega),|\alpha|\leq k\right\}$$

be standard Sobolev spaces with norm and seminorm given respectively by

$$\|v\|_{k}:=\left(\sum_{|\alpha|\leq k}\|\partial^{\alpha}v\|_{0}^{2}\right)^{1/2},\quad|v|_{k}:=\left(\sum_{|\alpha|=k}\|\partial^{\alpha}v\|_{0}^{2}\right)^{1/2}.$$

For $k=0$, $H^{0}(\Omega)$ is the standard $L^{2}(\Omega)$ space with the inner product denoted by $(\cdot,\cdot)$. Similarly, for any subset $K\subset\Omega$, $L^{2}(K)$ inner product is denoted by $(\cdot,\cdot)_{0,K}$. We note that, under the assumption (1.3),

$$a(v,v)\gtrsim\|v\|^{2}_{m,\Omega},\forall v\in H^{m}(\Omega).$$

(1.4)

The boundary value problem (1.1) can be cast into an equivalent optimization or a variational problem as described below for some approximate subspace $V\subset H^{m}(\Omega)$.

Minimization Problem M:

Find $u\in V$ such that

$$J(u)=\min_{v\in V}J(v),$$

(1.5)

or

Variational Problem V:

Find $u\in V$ such that

$$a(u,v)=\langle f,v\rangle\quad\forall v\in V.$$

(1.6)

The bilinear form $a(\cdot,\cdot)$ in (1.6), the objective function $J(\cdot)$ in (1.5) and the functional space $V$ depend on the type of boundary condition in (1.1).

One popular type of boundary conditions are Dirichlet boundary condition when $B^{k}=B_{D}^{k}$ are given by the following Dirichlet type trace operators

$$B_{D}^{k}(u):=\left.\frac{\partial^{k}u}{\partial\nu^{k}}\right|_{\partial\Omega}\quad(0\leq k\leq m-1),$$

(1.7)

with $\nu$ being the outward unit normal vector of $\partial\Omega$.

For the aforementioned Dirichlet boundary condition, the elliptic boundary value problem (1.1) is equivalent to (1.5) or (1.6) with $V=H^{m}_{0}(\Omega)$ and

$$a(u,v):=\sum_{|\alpha|=m}(a_{\alpha}\partial^{\alpha}u,\partial^{\alpha}v)_{0,\Omega}+(a_{0}u,v)\quad\forall u,v\in V,$$

(1.8)

and

$$J(v)=\frac{1}{2}a(v,v)-\int_{\Omega}fvdx.$$

(1.9)

Other boundary conditions such as Neumann boundary and mixed boundary conditions are a little bit complicated to describe for general case when $m\geq 2$ and will be discussed later.

Here we briefly review some classic finite element and other relevant methods for numerical solution of elliptic boundary value problems (1.1) for all $d,m\geq 1$.

Classic finite element methods use piecewise polynomial functions based on a given a subdivision, namely a finite element grid, of the domain, to discretize the variational problem. We will mainly review three different types of finite element methods: (1) conforming element method; (2) nonconforming and discontinuous Galerkin method; and (3) virtual element method.

Deep neural network (DNN), a tool developed for machine learning [Goodfellow et al., 2016]. DNN provides a very special function class that have been used for numerical solution of partial different equations, c.f.  [Lagaris et al., 1998]. By using different activation function such as sigmoidal, deep neural network can give rise to a very wide range of functional classes that can be drastically different from the piecewise polynomial function classes used in classic finite element method. One advantage of DNN approach is that it is quite easy to obtain smooth, namely $H^{m}$-conforming for any $m\geq 1$, DNN functions by simply choosing smooth activation functions. These function classes, however, do not usually form a linear vector space and hence the usual variational principle in classic finite element method can not be applied easily and instead collocation type of methods are often used. DNN is known to have much less “curse of dimensionality” than the traditional functional classes (such as polynomials or piecewise polynomials), DNN based method is potentially efficient to high dimensional problems and has been studied, for example, in [E et al., 2017] and [Sirignano and Spiliopoulos, 2018].

One main motivation of this paper is to explore DNN type methods that are most closely related to the traditional finite element methods. Namely we are interested in DNN function classes that consist of piecewise polynomials. By exploring relationship between DNN and FEM, we hope, on one hand, to expand or extend the traditional FEM approach by using new tools from DNN, and, on the other hand, to gain and develop theoretical insights and algorithmic tools into DNN by combining the rich mathematical theories and techniques in FEM.

In an earlier work [He et al., 2020b], we studied the relationship between deep neural networks (DNNs) using ReLU as activation function and continuous piecewise linear functions. One conclusion that can be drawn from [He et al., 2020b] is that any ReLU-DNN function is an $H^{1}$-conforming linear finite element function, and verse versa. The current work can be considered an extension of [He et al., 2020b] by considering using ReLU^k-DNN for high order partial differential equations. One focus in the current work is to provide error estimates when ReLU^k-DNN is applied to solve high order partial differential equations. More specifically, we will study a special class of $H^{m}$-conforming generalized finite element methods (consisting of piecewise polynomials) for (1.1) for any $m\geq 1$ and $d\geq 1$ based on artificial neural network for numerical solution of arbitrarily high order elliptic boundary value problem (1.1) and then provide convergence analysis for this method. For this type of method, the underlying finite element grids are not given a priori and the discrete solution can be obtained for solving a non-linear and non-convex optimization problem. In the case that the boundary of $\Omega\subset\mathcal{R}^{d}$, namely $\partial\Omega$, is curved, it is often an issue how to put a good finite element grid to accurately approximate $\partial\Omega$. As it turns out, this is not an issue for the finite neuron method, which is probably one of the advantages of the finite neuron method analyzed in this paper.

We note that the numerical method studied in this paper for elliptic boundary value problems is closely related to the classic finite element method, namely it amounts to piecewise polynomials with respect to an implicitly defined grid. We can also argue that it can be viewed as a mesh-less or even vertex-less method. But comparing with the popular meshless method, this method does correspond to some underlying grid, but this grid is not given a priori. This underlying grid is determined by the artificial neurons, which mathematically speaking refers to hyperplanes $\omega_{i}\cdot x+b_{i}=0$, together with a given activation function. By combining the names for finite element method and artificial neural network, for convenience of exposition, we will name the method studied in this paper as the finite neuron method.

The rest of the paper is organized as follows. In Section 2, we describe Monte-Carlo sampling technique, stratified sampling technique and Barron space. In Section 3, we construct the finite neuron functions and prove their approximation properties. In Section 5, we propose the finite neuron method and provide the convergence analysis. Finally, in Section 6, we give some summaries and discussions on the results in this paper.

Following [Xu, 1992], we will use the notation “$x\lesssim y$” to denote “$x\leq Cy$” for some constant $C$ independent of crucial parameter such as mesh size.

Let $\lambda\geq 0$ be a probability density function on a domain $G\subset\mathbb{R}^{D}(D\geq 1)$ such that

$$\int_{G}\lambda(\theta)d\theta=1.$$

(2.1)

We define the expectation and variance as follows

$$\mathbb{E}g:=\int_{G}g(\theta)\lambda(\theta)d\theta,\qquad\mathbb{V}g:=\mathbb{E}((g-\mathbb{E}g)^{2})=\mathbb{E}(g^{2})-(\mathbb{E}g)^{2}.$$

(2.2)

We note that

$$\displaystyle\mathbb{V}g\leq\max_{\theta,\theta^{\prime}\in G}(g(\theta)-g(\theta^{\prime}))^{2}.$$

For any subset $G_{i}\subset G$, let

$$\lambda(G_{i})=\int_{G_{i}}\lambda(\theta)d\theta,\qquad\lambda_{i}(\theta)={\lambda(\theta)\over\lambda(G_{i})}.$$

It holds that

$$\mathbb{E}_{G}g=\sum_{i=1}^{M}\lambda(G_{i})\mathbb{E}_{G_{i}}g.$$

For any function $h(\theta_{1},\cdots,\theta_{N}):G_{1}\times G_{2}\cdots G_{N}\mapsto\mathbb{R}$, define

$$\mathbb{E}_{G_{i}}g=\int_{G_{i}}g(\theta)\lambda_{i}(\theta)d\theta$$

and

$$\mathbb{E}_{N}h:=\int_{G_{1}\times G_{2}\times\ldots\times G_{N}}h(\theta_{1},\cdots,\theta_{N})\lambda_{1}(\theta_{1})\lambda_{2}(\theta_{2})\ldots\lambda_{N}(\theta_{N})d\theta_{1}d\theta_{2}\ldots d\theta_{N}.$$

(2.3)

For the Monte Carlo method, let $G_{i}=G$ for all $1\leq i\leq n$, namely,

$$\mathbb{E}_{N}h:=\int_{G\times G\times\ldots\times G}h(\theta_{1},\cdots,\theta_{N})\lambda(\theta_{1})\lambda(\theta_{2})\ldots\lambda(\theta_{N})d\theta_{1}d\theta_{2}\ldots d\theta_{N}.$$

(2.4)

The following result is standard [Rubinstein and Kroese, 2016] and their proofs can be obtained by direct calculations.

Stratified sampling [Bickel and Freedman, 1984] gives a more refined version of the Monte Carlo method.

Lemma 2.1 and Lemma 2.2 represent two simple identities and subsequent inequalities that can be verified by a direct calculation. Actually Lemma 2.1 is a special case of Lemma 2.2 with $M=1$. Lemma 2.1 and Lemma 2.2 are the basis of Monte-Carlo sampling and stratified sampling in statistics. In the presentation of the this paper, we choose not to use any concepts related to random samplings.

Given another domain $\Omega\subset\mathbb{R}^{d}$, we consider the case that $g(x,\theta)$ is a function of both $x\in\Omega$ and $\theta\in G$. Given any function $\rho\in L^{1}(G)$, we consider

$$u(x)=\int_{G}g(x,\theta)\rho(\theta)d\theta$$

(2.21)

with $\|\rho\|_{L^{1}(G)}<\infty$. Let $\lambda(\theta)={\rho(\theta)\over\|\rho\|_{L^{1}(G)}}$. Thus,

$$u(x)=\|\rho\|_{L^{1}(G)}\int_{G}g(x,\theta)\lambda(\theta)d\theta$$

(2.22)

with $\|\lambda(\theta)\|_{L^{1}(G)}=1$.

We can apply the above two lemmas to the given function $u(x)$.

Let us use a simple example to motivate the Barron space. Consider the Fourier transform of a real function $u:\mathbb{R}^{d}\mapsto\mathbb{R}$

$$\hat{u}(\omega)=(2\pi)^{-d/2}\int_{\mathbb{R}^{d}}e^{-i\omega\cdot x}u(x)dx.$$

(2.31)

This gives the following integral representation of $u$ in terms of the cosine function

$$u(x)=Re\int_{\mathbb{R}^{d}}e^{i\omega\cdot x}\hat{u}(\omega)d\omega=\int_{\mathbb{R}^{d}}\cos(\omega\cdot x+b(\omega))|\hat{u}(\omega)|d\omega,$$

(2.32)

where $\hat{u}(\omega)=e^{ib(\omega)}|\hat{u}(\omega)|$. Let

$$g(x,\omega)=\cos(\omega\cdot x+b(\omega))\quad\mbox{ and }\quad\rho(\omega)=|\hat{u}(\omega)|.$$

(2.33)

Thus,

$$u(x)=\int_{\mathbb{R}^{d}}g(x,\omega)\rho(\omega)d\omega,$$

(2.34)

If

$$\int_{\mathbb{R}^{d}}|\hat{u}(\omega)|d\omega<\infty,$$

then $\|\rho\|_{L^{1}}<\infty$. By applying the Lemma 2.3, there exist $\omega_{i}\in\mathbb{R}^{d}$ such that

$$\|u-u_{N}\|_{0,\Omega}\leq N^{-1/2}\|\hat{u}\|_{L^{1}(\mathbb{R}^{d})}.$$

(2.35)

where

$$u_{N}(x)={\|\hat{u}\|_{L^{1}(\mathbb{R}^{d})}\over N}\sum_{i=1}^{N}\cos(\omega_{i}\cdot x+b(\omega_{i}))$$

(2.36)

More generally, we consider the approximation property in $H^{m}$-norm. By (2.32),

$$\partial^{\alpha}u(x)=\int_{\mathbb{R}^{d}}\cos^{|\alpha|}(\omega\cdot x+b(\omega))\omega^{\alpha}|\hat{u}(\omega)|d\omega,\quad\forall\ |\alpha|\leq m.$$

(2.37)

For any positive integer $m$, let

$$g_{m}(x,\omega)={\cos(\omega\cdot x+b(\omega))\over(1+\|\omega\|)^{m}}\quad\mbox{and}\quad\rho_{m}(\omega)=(1+\|\omega\|)^{m}|\hat{u}(\omega)|,$$

(2.38)

where

$$\|\rho_{m}\|_{L^{1}(\mathbb{R}^{d})}=\int_{\mathbb{R}^{d}}(1+\|\omega\|)^{m}|\hat{u}(\omega)|d\omega<\infty.$$

Then, $\displaystyle u(x)=\int_{\mathbb{R}^{d}}g_{m}(x,\omega)\rho_{m}d\omega=\|\rho_{m}\|_{L^{1}(\mathbb{R}^{d})}\mathbb{E}g_{m}(x,\omega)$. Define

$$u_{N}(x)={\|\rho_{m}\|_{L^{1}(\mathbb{R}^{d})}\over N}\sum_{i=1}^{N}g_{m}(x,\omega_{i})={\|\rho_{m}\|_{L^{1}(\mathbb{R}^{d})}\over N}\sum_{i=1}^{N}{\cos(\omega_{i}\cdot x+b(\omega_{i}))\over(1+\|\omega_{i}\|)^{m}}.$$

(2.39)

It holds that

$$\partial^{\alpha}(u(x)-u_{N}(x))={\|\rho_{m}\|_{L^{1}(\mathbb{R}^{d})}\over N}\sum_{i=1}^{N}\mathbb{E}\partial^{\alpha}(g_{m}(x,\omega)-g_{m}(x,\omega_{i})).$$

By Lemma 2.1,

	$\displaystyle\mathbb{E}_{N}\sum_{|\alpha|\leq m}\|\partial^{\alpha}(u(x)-u_{N}(x))\|_{0,\Omega}^{2}$	$\displaystyle\leq\|\rho_{m}\|_{L^{1}(\mathbb{R}^{d})}^{2}\mathbb{E}_{N}\sum_{|\alpha|\leq m}\frac{1}{N^{2}}\sum_{i=1}^{N}\left(\mathbb{E}\partial^{\alpha}(g_{m}(x,\omega)-g_{m}(x,\omega_{i}))\right)^{2}$		(2.40)
		$\displaystyle\leq{\|\rho_{m}\|_{L^{1}(\mathbb{R}^{d})}^{2}\over N}\sum_{|\alpha|\leq m}\mathbb{E}\left(\partial^{\alpha}g_{m}(x,\omega)\right)^{2}$		(2.41)

Note that the definitions of $g_{m}$ and $\rho_{m}$ in (2.38) guarantee that

$$|\partial^{\alpha}g_{m}(x,\omega)|\leq 1.$$

Thus,

$$\mathbb{E}_{N}\sum_{|\alpha|\leq m}\|\partial^{\alpha}(u(x)-u_{N}(x))\|_{0,\Omega}^{2}\lesssim{\|\rho_{m}\|_{L^{1}(\mathbb{R}^{d})}^{2}\over N}.$$

This implies that there exist $\omega_{i}\in\mathbb{R}^{d}$ such that

$$\|u-u_{N}\|_{H^{m}(\Omega)}\lesssim N^{-1/2}\int_{\mathbb{R}^{d}}(1+\|\omega\|)^{m}|\hat{u}(\omega)|d\omega.$$

(2.42)