A Neural Network with Plane Wave Activation for Helmholtz Equation

We propose and study a plane wave activation based neural network (PWNN) for solving the Helmholtz boundary value problem

where $k>0$ is the wavenumber, $g\in L^{2}(\partial\Omega)$ is the boundary condition and $\Omega\subset{\mathbb{R}}^{2}$ is a lipschitz domain. We remark that the proposed approach is a basic technique for Helmholtz equation. To clarify the discussion, we only consider a simple case of Helmholtz equation, i.e. the Helmholtz boundary value problem (1). We will develop PWNN for some general Helmholtz equations, e.g. unbounded Helmholtz equation, Helmholtz scattering equation, at the future.

As a famous and significant partial differential equation (PDE), Helmholtz equation (1) appears in diverse scientific and engineering applications, since its solutions can represent the phenomenon of wave propagation, e.g. acoustic wave, electromagnetic wave, and seismic wave. The numerical methods for solving Helmholtz equation have attracted great interests of computational mathematicians, many numerical method like finite element method [24, 18], finite difference method [21], boundary element method [23, 27], plane wave method [2, 6, 13, 15], are developed for solving Helmholtz equation in different ways.

Since deep learning has made great success in engineering, researchers have become increasingly interested in using neural network to solve PDEs. Many relevant papers have been published around different topics in recent years. A deep Ritz method is presented by [11], it can use the variational formulation of elliptic PDEs and use the resulting energy as an objective function for the optimization of the parameters of the neural network. [32] convert the problem of finding the weak solution of PDEs into an operator norm minimization problem induced from the weak formulation, then parameterized the weak solution and the test function in the weak formulation as the primal and adversarial networks respectively, this method is called weak adversarial network (WAN). The physics-informed neural networks (PINN) is presented by [25], which take initial conditions and boundary conditions as the penalties of the optimization objective loss function. To deal with the noisy data, [31] proposes the Bayesian Physics-Informed Neural Networks (B-PINN), it use the Hamiltonian Monte Carlo or the variational inference as a posteriori estimator. Compared with PINNs, B-PINNs obtain more accurate predictions in scenarios with large noise. Unlike using variational formulation, this network needs high-order derivative information, it calculated the high-order differential in the framework of TensorFlow by using the backpropagation algorithm and automatic differentiation. The network can give a high accuracy solution for both forward and inverse problems. [28] proposes a deep Galerkin method. Due to the excessive computational cost of high-order differentiation of the neural network, the network bypasses this issue through the Monte Carlo method. [29] proposes sinusoidal representation networks (SIREN), in which sin(x) is used as activation function, and it performs better than other activation function like Sigmoid, Tanh or ReLU in Poisson equation and the Helmholtz and wave equations. In addition, deep learning tools are used to solve various practical problems related to PDE recently [26, 14, 5]. In conclusion, these existing works all focused on the design of the objective of using neural network to solve the general PDEs, not to develop new neural network for a specific PDE.

Motivated by the traditional activation based neural network (TANN), with Sigmoid, Tanh or ReLU as acitvation function [25], we introduce the PWNN method which use the plane wave function $e^{\mathbf{i}{x}}$ as activation function, since the plane wave is the basic component of the solution of Helmholtz equation. In addition, we find that the PWNN with one hidden layer is a more generalized form of the plane wave partition of unity method  [2, 6, 13, 15] (PWPUM). In the plane wave method, the direction of the basis function is closely related to the accuracy of the solution. In PWPUM, the direction of basis function can only be selected averagely. [4, 3] proposes NMLA to calculate the local direction of the true solution. A ray-FEM method is presented by [12], it first uses the standard finite element method to solve the low-frequency problem, then uses NMLA to calculate the local direction of the low-frequency numerical solution, and finally uses it as the direction of plane wave basis function to solve the high-frequency problem. [20] investigate the use of energy function to approximate the Helmholtz equation (LSM), the bases used in the discrete method are either plane waves of Bessel functions. [1] proposes a wave tracking strategy (LSM-WT), in each element, a same rotation angle is added to all plane wave basis functions with uniformly selected direction, by solving the double minimization problem of energy function in [20], the error can be effectively reduced. Based on this idea, we use neural network to construct this double minimization problem, using Monte Carlo method instead of integration. We call this method PWPUM-WT. Unlike using the plane wave basis with fixed directions in PWPUM, PWNN may combine several plane waves with learned amplitudes and directions by using plane wave acivation. As the direction of each basis function can be rotated by different angles, PWNN is more flexible than PWPUM-WT. Our experiment show the proposed PWNN method is more efficient than TANN, SIREN and higher accuracy than PWPUM, PWPUM-WT method.

This section introduces plane wave activation based neural network (PWNN) from the view of traditional activation based neural network (TANN). We first introduce a well known neural network based framework for solving PDEs, i.e. the physics-informed neural networks [25]. The entire neural network consists of $L+1$ layers, where layer $0$ is the input layer and layer $L$ is the output layer. Layers $0<l<L$ are the hidden layers. All of the layers have an activation function, excluding the output layer. As we known, the activation function can take the form of Sigmoid, Tanh (hyperbolic tangents) or ReLU (rectified linear units).

We denote $d_{0},d_{1},...,d_{L}$ as a list of integers, with $d_{0},d_{L}$ representing the lengths of input signal and output signal of the neural network. Define function ${\bf T}_{l}:\mathbb{R}^{d_{l}}\rightarrow\mathbb{R}^{d_{l+1}}$, $0\leq l<L$,

where ${\bf W}_{l}\in\mathbb{R}^{d_{l+1}\times d_{l}}$ and ${\bf b}_{l}\in\mathbb{R}^{d_{l+1}}$. Thus, we can simply represent a deep fully connected feedforward neural network using the composite function $h(\cdot;\theta):\mathbb{R}^{d_{0}}\rightarrow\mathbb{R}^{d_{L}}$,

where $\sigma$ is the activation function and $\theta:=\{{\bf W}_{l},{\bf b}_{l}\}$ represents the collection of all parameters. We only consider $d_{0}=2,d_{L}=1$, a two-dimensional scalar problem for Helmholtz equation (1).

Solving a general PDE defined as

by DNN is a physics-informed minimization problem with the objective $\mathcal{M}(\theta)$ consisting of two terms as follows:

where $\{{\bf x}_{f}^{i}\}_{i=1}^{N_{f}}$ and $\{{\bf x}_{g}^{i}\}_{i=1}^{N_{g}}$ are the collocation points in the inside and on the boundary of domain $\Omega$, respectively. $\lambda$ is a hyperparameter, by adjusting $\lambda$ to make the internal and boundary initial errors of the same order, usually there will be a better result. The domain term $\mathcal{M}_{\Omega}$ and boundary term $\mathcal{M}_{\partial\Omega}$ enforce the condition that the desired optimized neural network $h(\cdot;\theta^{*})$ satisfies $\mathcal{L}(u)=f$ and $\mathcal{B}(u)=g$, respectively. The effective and efficient stochastic gradient descents [19] with minibatches are recommended to solve the optimization problem (3).

For Helmholtz equation, since the solution is with complex value, we make a slight modification on the network architecture by redefining ${\bf W}_{L-1}\in\mathbb{C}^{d_{L}\times d_{L-1}},{\bf b}_{L-1}\in\mathbb{C}^{d_{L}}$. We may consider directly use the complex neural network [30] at the future. By this definition, we have $h(\cdot;\theta):\mathbb{R}^{d_{0}}\rightarrow\mathbb{C}^{d_{L}}$. Now, we let $\mathcal{L}=\Delta+k^{2},\mathcal{B}=\partial/\partial n+\mathbf{i}k$ and $f=0$, and let the activation function $\sigma(x)$ be ${\rm tanh}(x)$ or $e^{\mathbf{i}x}$ in minimization problem (3) to produce a TANN or PWNN based numerical method for solving Helmholtz equation (1) respectively.

As the key issue of numerically solving PDEs, a basic and interesting question is as follows: is it possible to use $h({\bf x};\theta)$ to approximate the solution of PDE? A well-known answer is the universal approximation property (UAP): if the solution $u({\bf x})$ is bounded and continuous, then $h({\bf x};\theta)$ can approximate $u({\bf x})$ to any desired accuracy, given the increasing hidden neurons [17, 9, 16]. We should notice that this conclusion required a simple assumption on the activation function, i.e. $\sigma(x)$ should be bounded and monotonically-increasing function. However, it’s easy to check that the plane wave activation doesn’t satisfy this assumption. We open a theoretical problem that whether PWNN has the UAP at the future research.

In this section , we will try to bridges PWNN to the plane wave partition of unity method (PWPUM), which is the famous and classical numerical method with complexity and approximation analysis. If we consider a simpler PWNN with one hidden layer, let $\mathbf{T}_{0}:=(w_{1},w_{2},\cdots,w_{d_{1}})^{T}$, $\mathbf{b}_{0}:=(b_{1},b_{2},\cdots,b_{d_{1}})^{T}$, $\mathbf{T}_{1}:=(s_{1},s_{2},\cdots,s_{d_{1}})$, and $w_{i}\in{\mathbb{R}}^{2}$, $b_{i}\in{\mathbb{R}}$, $s_{i}\in\mathbb{C}$, we will have the output function with the form

Notice that $e^{\mathbf{i}b_{i}}$ can be equivalently learned by learning $s_{i}$, we can ignore $b_{i}$ in (4) to obtain an equivalent output formulation

with $\theta:=\{s_{i},w_{i}\}_{i=1}^{d_{1}}$. So $\theta$ of PWNN with one hidden layer is determined by solving the following problem,

Here we write the loss function in integral form. In fact, we use Monte Carlo method instead of integration in neural network. In the standard finite element method or other methods, there will also be errors in numerical integration. This part of the error is only a small part when not pursuing too high accuracy, e.g. relative $L^{2}$ error less than 1e-8.

Now we come to introduce the plane wave partition of unity method  [2, 6, 13, 15] (PWPUM). For the Helmholtz boundary value equation with wavenumber $k$, $1\leq i\leq d_{1}$, let

the PWPUM is to find a numerical solution with the form

where we define $\beta:=\{s_{i}\}_{i=1}^{d_{1}}$ the parameters which need to be optimized. We define the discrete function space with form (8) depending on $\beta$ as $V_{h}$. Instead of directly solving an physics-informed optimization problem (3) in PWNN, PWPUM instead find $u\in V_{h}$, such that, $\forall\,v\in V_{h}$,

By the formulation of $u$ in (8) and the definition of $k_{i}$ in (7), we have $\Delta u+k^{2}u=0$, yielding that the first volume integral term of (9) vanishes. Since the Helmholtz operator and the discrete function space are both linear, we can take the test function $v=e^{\mathbf{i}k_{i}^{T}\mathbf{x}}$, $1\leq i\leq d_{1}$. Then, to find $u\in V_{h}$ satisfying (9) is to solve a linear system $\mathbf{M}\beta=\mathbf{G}$ with $\mathbf{M}\in\mathbb{C}^{d_{1}\times d_{1}},\mathbf{G}\in\mathbb{C}^{d_{1}}$,and

Actually, one may check that the above linear system is equivalent to solving the following problem:

In fact, the basis function of PWPUM is plane wave function and the wave number is strictly equal to $k$. The internal integral error is strictly equal to 0, so this optimization problem is independent of $\lambda$.

For the PWPUM method, the direction of the plane wave basis function is uniformly selected because the prior information of the solution is not known. Because any direction may be the direction of the exact solution, PWPUM may need a lot of basis functions to achieve the required accuracy. In order for PWPUM to use the direction that can be changed, [1] proposes a wave tracking strategy (LSM-WT), in each element, a same rotation angle is added to all plane wave basis functions with uniformly selected direction. Here we can also use this wave tracking strategy in one element, the form of solution (8) becomes:

where

and the optimization problem (10) becomes:

In addition, if we limit $w_{i}$ in 5 to $w_{i}^{T}=k_{i}^{T}(\alpha)$, then this is the form of PWPUM-WT, the output formulation $h(\mathbf{x};\theta)$ is equal to $u(\mathbf{x};\beta,\alpha)$ in (11), with $\theta:=\{s_{i},\alpha\}_{i=1}^{d_{1}}$.

Apparently, comparing PWNN (5), PWPUM (8) and PWPUM-WT (11) will show that PWPUM-WT is a generalization of PWPUM. Furthermore, PWNN is a generalization of PWPUM-WT. Unlike using the plane wave basis with fixed directions (7) in PWPUM, PWNN will learn a group of optimized plane waves $e^{\mathbf{i}w_{i}^{T}\mathbf{x}}$ in (5) with learned both amplitude $||w_{i}||$ and direction ${\rm arccos}(w_{i,1}/||w_{i}||)$. Since the approximation analysis of PWPUM was well addressed, see Theorem 3.9 of [15], and the optimization problem of PWNN (6) is obviously a generalization of that in PWPUM (10), it’s trivial to prove that the approximation error of PWNN with one hidden layer is no more than that of PWPUM theoretically. Actually, we can only consider learning the directions by adding a regularization term, like $\sum_{i=1}^{d_{1}}(||w_{i}||-k)^{2}$ in the minimization objective (3). However, motivated by the dispersion correction technique [7, 8], i.e. using different wavenumber in the discretization scheme may reduce the pollution effect [10] of the Helmholtz problems with high wavenumber, we ignore adding this regularization term.

We experimental investigate the performance of PWNN from two aspects: the exact solution of equation (1) is a combination of plane waves with all known directions (KD) or unknown directions (UD). Apparently, the UD problem is more common in practical problems. For instance, if we consider a Helmholtz scattering equation with complicated scatters, i.e. the propagation of acoustic wave in a complicated situation, we would not expect to accurately predict the directions of the multi reflection waves in advance. Considering that the sin activation based neural network (SIREN) performs better than the traditional activation based neural network (TANN) in some cases, we add the comparison between PWNN and SIREN in the experiments. The experimental results on problems KD and UD demonstrate that: PWNN works much better than TANN and SIREN, i.e. the introduction of plane wave as activation unit is indeed helpful for Helmholtz equation; PWPUM-WT has limited advantages over PWPUM, i.e. adding only one degree of freedom of rotation angle is not enough. PWNN has competitive performance with PWPUM under problem KD but outperforms PWPUM under problem UD, i.e. PWNN is indeed a generalization of PWPUM but more practical. We introduce the experimental setting first.

Experimental setting. We simply set the network architectures as having equal units for each layer. The $Layers$ used later indicates the hidden layers in the network, i.e. $Layers:=L-1$. The $Units$ indicates the units per layer, i.e. $d_{1}=d_{2}=\cdots=d_{L-1}=Units$. We chose Tanh as activation function in TANN. In both TANN and PWNN, we use limited-memory BFGS algorithm [22] to update parameters. L-BFGS is an improved algorithm for quasi-Newton method: BFGS. In BFGS, the approximate Hesse matrix $B_{k}^{-1}$ is stored at every step, which wastes a lot of storage space in high dimensional cases. In L-BFGS, only the recent $m$ steps’ iterative information is saved for calculation $B_{k}^{-1}$ to reduce the storage space of data. Here we set $m=50$. The stop criterion of the inner iteration is $\|\nabla\mathcal{M}\|_{\infty}<2\times 10^{-16}$, or that we exceeded the maximum number of allowed iterations, set as 50000 here. $N_{f},N_{g}$ denotes the number of training points in the inside and on the boundary. The training points in the inside are randomly selected, while the training points on the boundary are uniformly selected. The test points used to estimate the relative $L^{2}$ errors are uniformly sampled by row and column, with 10000 in all. We remark that this paper ignores the detailed discussion of computational cost, since all the numerical examples only take up to several minutes to get satisfying results. We may consider comparing the computational cost of TANN, PWNN, and PWPUM for a 3D Helmholtz problem at the future.

This paper proposed a plane wave activation based neural network (PWNN) for Helmholtz equation. We compared PWNN to traditional activation based neural network (TANN), sin based neural network (SIREN) and plane wave partition of unity method (PWPUM). We draw lessons from LSM-WT and propose PWPUM-WT, which can be regarded as a simplified PWNN. The experimental results can guarantee the improvement of introducing plane wave activation for Helmholtz equation and the improvement of PWNN on PWPUM for more practical problems. This paper only provide a basic technique for Helmholtz equation. We will develop PWNN for some general Helmholtz equations at the future.

This research was partially supported by National Natural Science Foundation of China with grant 11831016 and was partially supported by Science Challenge Project, China TZZT2019-B1.1.