A multiscale neural network based on hierarchical matrices

Introduction

In the past decades, there has been a great combined effort in developing efficient algorithms to solve linear problems issued from discretization of integral equations (IEs), and partial differential equations (PDEs). In particular, multiscale methods such as multi-grid methods , the fast multipole method , wavelets , and hierarchical matrices , have been strikingly successful in reducing the complexity for solving such systems. In several cases, for operators of pseudo-differential type, these algorithms can achieve linear or quasi-linear complexity. In a nutshell, these methods aim to use the inherent multiscale structure of the underlying physical problem to build efficient representations at each scale, thus compressing the information contained in the system. The gains in complexity stem mainly from processing information at each scale, and merging it in a hierarchical fashion.

Even though these techniques have been extensively applied to linear problems with outstanding success, their application to nonlinear problems has been, to the best of our knowledge, very limited. This is due to the high complexity of the solution maps. In particular, building a global approximation of such maps would normally require an extremely large amount of parameters, which in return, is often translated to algorithms with a prohibitive computational cost. The development of algorithms and heuristics to reduce the cost is an area of active research . However, in general, each method is application-dependent, and requires a deep understanding of the underlying physics.

On the other hand, the surge of interest in machine learning methods, in particular, deep neural networks, has dramatically improved speech recognition , visual object recognition , object detection, etc. This has fueled breakthroughs in many domains such as drug discovery , genomics , and automatic translation , among many others . Deep neural networks have empirically shown that it is possible to obtain efficient representations of very high-dimensional functions. Even though the universality theorem holds for neural networks , i.e., they can approximate arbitrarily well any function with mild regularity conditions, how to efficiently build such approximations remains an open question. In particular, the degree of approximation depends dramatically on the architecture of the neural network, i.e. how the different layers are connected. In addition, there is a fine balance between the number of parameters, the architecture, and the degree of approximation .

This paper aims to combine the tools developed in deep neural networks with ideas from multiscale methods. In particular, we extend hierarchical matrices ($\mathcal{H}$-matrices) to nonlinear problems within the framework of neural networks. Let

be a nonlinear generalization of pseudo-differential operators, issued from an underlying physical problem, described by either an integral equation or a partial differential equation, where $v$ can be considered as a parameter in the equation, $u$ is either the solution of the equation or a function of it, and $n$ is the number of variables.

We build a neural network with a novel multiscale structure inspired by hierarchical matrices. We interpret the application of an $\mathcal{H}$-matrix to a vector using a neural network structure as follows. We first reduce the dimensionality of the vector, or restrict it, by multiplying it by a short and wide structured matrix. Then we process the encoded vector by multiplying it by a structured square matrix. Then we return the vector to its original size, or interpolate it, by multiplying it by a structured tall and skinny matrix. Such operations are performed separately at different spatial scales. The final approximation to the matrix-vector multiplication is obtained by adding the contributions from all spatial scales, including the near-field contribution, which is represented by a near-diagonal matrix. Since every step is linear, the overall operation is also a linear mapping. This interpretation allows us to directly generalize the $\mathcal{H}$-matrix to nonlinear problems by replacing the structured square matrix in the processing stage by a structured nonlinear network with several layers. The resulting artificial neural network, which we call multiscale neural network, only requires a relatively modest amount of parameters even for large problems.

We demonstrate the performance of the multiscale neural network by approximating the solution to the nonlinear Schrödinger equation , as well as the Kohn-Sham map . These mappings are highly nonlinear, and are still well approximated by the proposed neural network, with a relative accuracy on the order of $10^{-4}\sim 10^{-3}$. Furthermore, we do not observe overfitting even in the presence of a relatively small set of training samples.

Although machine learning and deep learning literature is vast, the application of deep learning to numerical analysis problems is relatively new, though that is rapidly changing. Research using deep neural networks with multiscale architectures has primarily focused on image and video processing.

Deep neural networks have been recently used to solve PDEs and classical inverse problems . For general applications of machine learning to nonlinear numerical analysis problems, the work of Raissi and Karnidiakis used machine learning, in particular, Gaussian processes, to find parameters in nonlinear equations ; Chan, and Elsheikh predicted the basis function on the coarse grid in multiscale finite volume method by neural network; Khoo, Lu and Ying used neural network in the context of uncertainty quantification ; Zhang et al used neural network in the context of generating high-quality interatomic potentials for molecular dynamics . Wang et al. applied non-local multi-continuum neural network on time-dependent nonlinear problems . Khrulkov et al. and Cohen et al. developed deep neural network architectures based on tensor-train decomposition . In addition, we note that deep neural networks with related multi-scale structures have been proposed mainly for applications such as image processing, however, we are not aware of any applications of such architectures to solving nonlinear differential or integral equations.

2 Organization

The reminder of the paper is organized as follows. Section 2 reviews the $\mathcal{H}$-matrices and interprets the $\mathcal{H}$-matrices using the framework of neural networks. Section 3 extends the neural network representation of $\mathcal{H}$-matrices to the nonlinear case. Section 4 discusses the implementation details and demonstrates the accuracy of the architecture in representing nonlinear maps, followed by the conclusion and future directions in Section 5.

Neural network architecture for ℋ\mathcal{H}-matrices

In this section, we aim to represent the matrix-vector multiplication of $\mathcal{H}$-matrices within the framework of neural networks. For the sake of clarity, we succinctly review the structure of $\mathcal{H}$-matrices for the one dimensional case in Section 2.1. We interpret $\mathcal{H}$-matrices using the framework of neural networks in Section 2.2, and then extend it to the multi-dimensional case in Section 2.3.

Hierarchical matrices ($\mathcal{H}$-matrices) were first introduced by Hackbusch et al. in a series of papers as an algebraic framework for representing matrices with a hierarchically off-diagonal low-rank structure. This framework provides efficient numerical methods for solving linear systems arising from integral equations (IE) and partial differential equations (PDE) . In the sequel, we follow the notation in to provide a brief introduction to the framework of $\mathcal{H}$-matrices in a simple uniform and Cartesian setting. The interested readers are referred to for further details.

where $u$ and $v$ are periodic in $\Omega$ and $g(x,y)$ is smooth and numerically low-rank away from the diagonal. A discretization with an uniform grid with $N=2^{L}m$ discretization points yields the linear system given by

We introduce a hierarchical dyadic decomposition of the grid in $L+1$ levels as follows. We start by the $0$-th level of the decomposition, which corresponds to the set of all grid points defined as

Fig. 2 illustrates the block partition of $A$ induced by the dyadic partition, and the decomposition induced by the different interaction lists at each level that follows (2.4).

Thus the matrix-vector multiplication (2.2) can be expressed as

as illustrated in Fig. 3, which constitutes the basis for writing $\mathcal{H}$-matrices as a neural network.

2 Matrix-vector multiplication as a neural network

An artificial neural network, in particular, a feed-forward network, can be thought of the composition of several simple functions, usually called propagation functions, in which the intermediate one-dimensional variables are called neurons, which in return, are organized in vector, or tensor, variables called layers. For example, an artificial feed-forward neural network

with $K+1$ layer can be written using the following recursive formula

Within this context, training a network refers to finding the weights and biases, whose entries are collectively called parameters, in order to approximate a given map. This is usually done by minimizing a loss function using a stochastic optimization algorithm.

We interpret the structure of $\mathcal{H}$-matrices (2.6) using the framework of neural networks. The different factors in (2.7) possess a distinctive structure, which we aim to exploit by using locally connected (LC) network. LC networks are propagating functions whose weights have a block-banded constraint. For the one-dimensional example, we also treat $\xi$ as a 2-tensor of dimensions $\alpha\times N_{x}$, where $\alpha$ is the channel dimension and $N_{x}$ is the spatial dimension, and $\zeta$ be a 2-tensor of dimensions $\alpha^{\prime}\times N_{x}^{\prime}$. We say that $\xi$ is connected to $\zeta$ by a LC networks if

Each LC network requires $6$ parameters, $N_{x}$, $\alpha$, $N_{x}^{\prime}$, $\alpha^{\prime}$, $w$ and $s$ to be characterized. Next, we define three types of LC network by specifying some of their parameters,

Restriction network: we set $s=w=\frac{N_{x}}{N_{x}^{\prime}}$ and $\alpha=1$ in LC. This network represents the multiplication of a block diagonal matrix with block sizes $\alpha^{\prime}\times s$ and a vector with size $N_{x}\alpha$, as illustrated by Fig. 4 (a). We denote this network using $\mathsf{LCR}[\phi;N_{x},N_{x}^{\prime},\alpha^{\prime}]$. The application of $\mathsf{LCR}[\mathsf{linear};32,8,3]$ is depicted in Fig. 4 (a).

Kernel network: we set $s=1$ and $N_{x}^{\prime}=N_{x}$. This network represents the multiplication of a cyclic block band matrix of block size $\alpha^{\prime}\times\alpha$ and band size $\frac{w-1}{2}$ times a vector of size $N_{x}\alpha$, as illustrated by the upper portion of Fig. 4 (b). To account for the periodicity we pad the input layer $\xi_{c,j}$ on the spatial dimension to the size $(N_{x}+w-1)\times\alpha$. We denote this network by $\mathsf{LCK}[\phi;N_{x},\alpha,\alpha^{\prime},w]$. This network has two steps: the periodic padding of $\xi_{c,j}$ on the spatial dimension, and the application of (2.10). The application of $\mathsf{LCK}[\mathsf{linear};8,3,3,3]$ is depicted in Fig. 4 (b).

Interpolation network: we set $s=w=1$ and $N_{x}^{\prime}=N_{x}$ in LC. This network represents the multiplication of a block diagonal matrix with block size $\alpha^{\prime}\times\alpha$, times a vector of size $N_{x}\alpha$, as illustrated by the upper figure in Fig. 4 (c). We denote the network $\mathsf{LCI}[\phi;N_{x},\alpha,\alpha^{\prime}]$, which has two steps: the application of (2.10), and the reshaping of the output to a vector by column major indexing. The application of $\mathsf{LCI}[\mathsf{linear};8,3,4]$ is depicted in Fig. 4 (c).

2.2 Neural network representation

Following (2.7), in order to construct the neural network (NN) architecture for (2.7), we need to represent the following four operations:

From 1.1 and the definition of $\mathsf{LCR}$ and $\mathsf{LCI}$, the operations (2.11a) and (2.11c) are equivalent to

respectively. Analogously, 1.3 indicates that (2.11b) is equivalent to

Combining (2.12), (2.13) and (2.14), we obtain Algorithm 1, whose architecture is illustrated in Fig. 5. In particular, Fig. 5 is the translation to the neural network framework of (2.7) (see Fig. 3) using the building blocks depicted in Fig. 4.

Moreover, the memory footprints of the neural network architecture and $\mathcal{H}$-matrices are asymptotically the same with respect the spatial dimension of $u$ and $v$. This can be readily shown by computing the total number of parameters. For the sake of simplicity, we only count the parameters in the weights, ignoring those in the biases. A direct calculation yields the number of parameters in $\mathsf{LCR}$, $\mathsf{LCK}$ and $\mathsf{LCI}$:

respectively. Hence, the number of parameters in Algorithm 1 is

3 Multi-dimensional case

Following the previous section, the extension of Algorithm 1 to the $d$-dimensional case can be easily deduced using the tensor product of one-dimensional cases. Consider $d=2$ below for instance, and the generalization to higher dimensional case will be straight-forward. Suppose that we have an IE in 2D given by

we discretize the domain $\Omega$ with a uniform grid with $n=N^{2}$ ($N=2^{L}m$) discretization points, and let $A$ be the resulting matrix obtained from discretizing (2.17). We denote the set of all grid points as

where $n_{b}=(n_{b,1},n_{b,2})$ is called the band size for tensor. Thus 1 can be generalized to tensors yielding the following properties.

Next, we characterize LC networks for the 2D case. An NN layer for 2D can be represented by a 3-tensor of size $\alpha\times N_{x,1}\times N_{x,2}$, in which $\alpha$ is the channel dimension and $N_{x,1}$, $N_{x,2}$ are the spatial dimensions. If a layer $\xi$ with size $\alpha\times N_{x,1}\times N_{x,2}$ is connected to a locally connected layer $\zeta$ with size $\alpha^{\prime}\times N_{x,1}^{\prime}\times N_{x,2}^{\prime}$, then

where $(i-1)s=((i_{1}-1)s_{1},(i_{2}-1)s_{2})$. As in the 1D case, the channel dimension corresponds to the rank $r$, and the spatial dimensions correspond to the grid points of the discretized domain. Analogously to the 1D case, we define the LC networks $\mathsf{LCR}$, $\mathsf{LCK}$ and $\mathsf{LCI}$ and use them to express the four operations in (2.11) which constitute the building blocks of the neural network. The extension is trivial, the parameters $N_{x}$, $s$ and $w$ in the one-dimensional LC networks are replaced by their 2-dimensional counterpart $N_{x}=(N_{x,1},N_{x,2})$, $s=(s_{1},s_{2})$ and $w=(w_{1},w_{2})$, respectively. We point out that $s=w=\frac{N_{x}}{N_{x}^{\prime}}$ for the 1D case is replaced by $s_{j}=w_{j}=\frac{N_{x,j}}{N_{x,j}^{\prime}}$, $j=1,2$ for the 2D case in the definition of LC.

Using the notations above we extend Algorithm 1 to the 2D case in Algorithm 2. We crucially remark that the $\mathsf{Reshape}[r^{2},n_{1},n_{2}]$ function in Algorithm 2 is not the usual major column based reshaping. It reshapes a 2-tensor $T$ with size $rn_{1}\times rn_{2}$ to a 3-tensor $S$ with size $r^{2}\times n_{1}\times n_{2}$, by treating the former as a block tensor with block size $r\times r$, and reshaping each block as a vector following the formula $S(k,i,j)=T((i-1)r+k_{1},(j-1)r+k_{2})$ with $k=(k_{1}-1)r+k_{2}$, for $k_{1},k_{2}=1,\dots,r$, $i=1,\dots,n_{1}$ and $j=1,\dots,n_{2}$. Fig. 6 provides an example for the case $\mathsf{Reshape}[2^{2},3,3]$. The $\mathsf{Flatten}$ is its inverse.

Multiscale neural network

In this section, we extend the aforementioned NN architecture to represent a nonlinear generalization of pseudo-differential operators of the form

Due to its multiscale structure, we refer to the resulting NN architecture as the multiscale neural network (MNN). We consider the one-dimensional case below for simplicity, and the generalization to higher dimensions follows directly as in Section 2.3.

NN can represent nonlinearities by choosing the activation function, $\phi$, to be nonlinear, such as ReLU or sigmoid. The range of the activation function also imposes constraints on the output of the NN. For example, the range of “ReLU” in $[0,\infty)$ and the range of the sigmoid function is $$. Thus, the last layer is often chosen to be a linear layer to relax such constraint. Algorithm 1 is then revised to Algorithm 3, and the architecture is illustrated in Fig. 7. We remark that the nonlinear activation function is only used in the$\mathsf{LCK}$network. The$\mathsf{LCR}$and$\mathsf{LCI}$ networks in Algorithm 1 are still treated as restriction and interpolation operations between coarse grid and fine grid, respectively, so we use the linear activation functions in these layers. Particularly, we also use the linear activation function for the last layer of the adjacent part, which is marked in line 13 in Algorithm 3.

As in the linear case, we calculate the number of parameters of MNN and obtain (neglecting the number of parameters in $b$ in (2.10))

2 Translation-invariant case

Compared to the LC network, the only difference is that the parameters $W$ and $b$ are independent of $i$. Hence, inheriting the definition of $\mathsf{LCR}$, $\mathsf{LCK}$ and $\mathsf{LCI}$, we define the layers $\mathsf{CR}$, $\mathsf{CK}$ and $\mathsf{CI}$, respectively. By replacing the LC layers in Algorithm 1 by the corresponding CNN layers, we obtain the neural network architecture for the translation invariant kernel.

For the nonlinear case, the translation invariant kernel for the linear case can be extended to kernels that are equivariant to translation, i.e. for any translation $\mathcal{T}$,

For this case, all the LC layers in Algorithm 3 can be replaced by its corresponding CNN layers. The number of parameters of $\mathsf{CR}$, $\mathsf{CK}$ and $\mathsf{CI}$ are

Thus, the number of parameters in Algorithm 3 using CNN is

Numerical results

In this section we discuss the implementation details of MNN. We demonstrate the accuracy of the MNN architecture using two nonlinear problems: the nonlinear Schrödinger equation (NLSE), and the Kohn-Sham map (KS map) in the Kohn-Sham density functional theory (KSDFT).

where $u$ is the target solution generated by a numerical discretization of the PDEs and $u_{NN}$ is the predicted solution by MNN. The optimization is performed using the NAdam optimizer . The weights and biases in MNN are initialized randomly from the normal distribution and the batch size is always set between $1/100$th and $1/50$th of the number of train samples.

2 NLSE with inhomogeneous background potential

The nonlinear Schrödinger equation (NLSE) is widely used in quantum physics to describe the single particle properties of the Bose-Einstein condensation phenomenon . Here we study the NLSE with inhomogeneous background potential $V(x)$:

with period boundary condition. We aim to find its ground state denoted by $u_{G}(x)$. We take a strongly nonlinear case $\beta=10$ in this work and thus consider a defocusing cubic Schrödinger equation. Due to the cubic term, an iterative method is required to solve (4.2) numerically. We employ the method in for the numerical solution, which solves a time-dependent NLSE by a normalized gradient flow. The MNN is used to learn the map from the background potential to the ground state

This map is equivariant to translation, and thus MNN is implemented using the CNN layers. The constraints in (4.2) can be guaranteed by adding a post-correction in the network as

where $\hat{u}$ is the prediction of the neural network. In the following, we study the performance of MNN on 1D and 2D cases.

For the one-dimensional case, the number of discretization points is $N=320$, and we set $L=6$ and $m=\frac{N}{2^{L}}=5$. The potential $V$ is chosen as

We also compare the MNN with an instance of the classical convolutional neural networks (CNN). The CNN architecture used in is adopted as the reference architecture. Table 4 presents the setup of the networks and the training and validation errors for CNN. Clearly, by comparing the results of MNN in Tables 2 and 3 with the results of CNN in Table 4, one can observe that the MNN not only reduces the number of parameters, but also improve the accuracy.

throughout the results shown in Tables 2 and 3, the validation errors are very close to the corresponding training errors, thus no overfitting is observed. Fig. 8 presents a sample for the potential $V$ and its corresponding solution and prediction solution by MNN. We can observe that the prediction solution agrees with the target solution very well.

2.2 Two-dimensional case

For the two-dimensional case, the number of discretization points is $n=80\times 80$, and we set $L=4$ and $m=5$. The potential $V$ is chosen as

Tables 5 and 6 present the numerical results for different number of channels, $r$, and different number of $\mathsf{CK}$ layers ,$K$, respectively. Similarly to the 1D case, the choice of parameters $r=6$ and $K=7$ also yield accurate results in the 2D case. Fig. 9 presents a sample of the potential $V$ in the test set and its corresponding solution, prediction solution and its error with respect to the reference solution.

3 Kohn-Sham map

Kohn-Sham density functional theory is the most widely used electronic structure theory. It requires the solution of the following set of nonlinear eigenvalue equations (real arithmetic assumed for all quantities):

Here $n_{e}$ is the number of electrons (spin degeneracy omitted), $d$ is the spatial dimension, and $\delta_{ij}$ stands for the Kronecker delta. In addition, all eigenvalues $\{\varepsilon_{i}\}$ are real and ordered non-decreasingly, and $\rho(x)$ is the electron density, which satisfies the constraint

The Kohn-Sham equations (4.7) need to be solved self-consistently, which can also viewed as solving the following fixed point map

Here the mapping $\mathcal{F}_{\text{KS}}[\cdot]$ from $V$ to $\rho$ is called the Kohn-Sham map, which for a fixed potential is reduced to a linear eigenvalue problem, and it constitues the most computationally intensive step for solving (4.7). We seek to approximate the Kohn-Sham map using a multiscale neural network, whose output was regularized so it satisfies (4.8).

In the following numerical experiments the potential, $V$, is given by

where $d$ is the dimension and $r_{i}\in[-1,1)^{d}$. We set $\sigma=0.05$ or $\sigma=0.15$ for 1D and $\sigma=0.2$ for 2D. The coefficients $c_{i}$ are randomly chosen following the uniform distribution $\mathcal{U}([0.8,1.2])$, and the centers of the Gaussian wells $r_{i}$, are chosen randomly under the constraint that $|r_{i}-r_{i^{\prime}}|>4\sigma$, unless explicitely specified. The Kohn-Sham map is discretized using a pseudo-spectral method , and solved by a standard eigensolver.

We set $\sigma=0.05$ and we generated $4$ data sets using different number of wells, $n_{g}$, which in this case is also equal to the number of electrons $n_{e}$, ranging from $2$ to $8$.

The number of discretization points is $N=320$. We trained the architecture defined in Section 3 for each $n_{g}$, setting the number of levels $L=6$, using different values for $r$ and $K$.

Table 7 shows that there is no overfitting, even at this level of accuracy and number of parameters. This behavior is found in all the numerical examples, thus we only report the test error in what follows.

From Table 8 we can observe that as we increase $r$ the error decreases sharply. Fig. 10 depict this behavior. In Fig. 10 we have that if $r=2$, then the network output $\rho_{NN}$, fails to approximate $\rho$ accurately; however, by modestly increasing $r$, the network is able to accurately approximate $\rho$.

However, the accuracy of the network stagnates rapidly. In fact, increasing $r$ beyond $10$ does not provide any considerable gains. In addition, Table 8 shows that the accuracy of the network is agnostic to the number of Gaussian wells present in the system.

In addition, we studied the relation between the quality of the approximation and $K$. We fixed $r=6$, and we trained several networks using different values of $K$, ranging from $2$, i.e., a very shallow network, to $10$. The results are summarized in Table 9. We can observe that the error decreases sharply as the depth of the network increases, but it rapidly stagnates as $K$ becomes large.

The Kohn-Sham map is a very non-linear mapping, and we demonstrate below that the non-linear activation functions in the network play a crucial role. Consider a linear network, in which we have completely eliminated the non-linear activation functions. This linear network can be easily implemented by setting $K=1$, and $\phi(x)=x$. We then train the resulting linear network using the same data as before and following the same optimization procedure. Table 10 shows that the approximation error can be very large, even for relatively small $n_{g}$.

In addition, we note that the favorable behavior of MNN with respect to $n_{g}$ shown in the prequel is rooted in the fact that the band gap (here the band gap is equal to $\epsilon_{n_{g}+1}-\epsilon_{n_{g}}$) remains approximately the same as we increase $n_{g}$. According to the density functional perturbation theory (DFPT), the Kohn-Sham map is relatively insensitive to the change of the external potential. This setup mimics an insulating system. On the other hand, we may choose $\sigma$ in (4.10) to be $0.15$, and relax the constraint between Gaussian centers to $|r_{i}-r_{i^{\prime}}|>2\sigma$. The rest of the coefficients are randomly chosen using the same distributions as before. In this case, we generate a new data set, in which the average gap for the generated data sets are $0.2$, $0.08$, $0.05$ and $0.03$ for $n_{g}$ equal to $2$, $4$, $6$, and $8$, respectively. The decrease of the band gap with respect to the increase of $n_{g}$ resembles the behavior of a metallic system, in which the the Kohn-Sham map becomes more sensitive to small perturbations of the potential. After generating the samples, we trained two different networks: the MNN with $K=6$, $r=10$ and a regular CNN with $15$ layers, $10$ channels, window size $13$. The results are shown in Tables 11 and 12. From Tables 11 and 12 we can observe that MNN outperforms CNN and that as the band gap decreases, the performance gap between MNN and CNN widens. We point out that it is possible to partially alleviate this adverse dependence on the band gap in the MNN by introducing a nested hierarchical structure to the interpolation and restriction operators as shown in .

3.2 Two-dimensional case

The discretization is the standard extension to 2D using tensor products, using a $64\times 64$ grid. In this case we only used $n_{g}=2$ and we followed the same number of training and test samples as that in the $1$D case. We fixed $K=6$, $L=4$, and we trained the network for different number of channels, $r$. The results are displayed in Table 13, which shows the same behavior as for the 1D case, in which the error decays sharply and then stagnates, and there is no over fitting. In particular, the network is able to effectively approximate the Kohn-Sham map as shown in Fig. 11. Fig. 11a shows the output of neural network for a test sample and Fig. 11b shows the approximation error with respect to the reference.

Conclusion

We have developed a multiscale neural network (MNN) architecture for approximating nonlinear mappings, such as those arising from the solution of integral equations (IEs) or partial differential equations (PDEs). In order to control the number of parameters, we first rewrite the widely used hierarchical matrix into the form of a neural network, which mainly consists of three sub-networks: restriction network, kernel network, and interpolation network. The three sub-networks are all linear, and correspond to the components of a singular value decomposition. We demonstrate that such structure can be directly generalized to nonlinear problems, simply by replacing the linear kernel network by a multilayer kernel network with nonlinear activation functions. Such “nonlinear singular value decomposition operation” is performed at different spatial scales, which can be efficiently implemented by a number of locally connected (LC) networks, or convolutional neural networks (CNN) when the mapping is equivariant to translation. Using the parameterized nonlinear Schrödinger equation and the Kohn-Sham map as examples, we find that MNN can yield accurate approximation to such nonlinear mappings. When the mapping has $N$ degrees of freedom, the complexity of MNN is only $O(N\log N)$. Thus the resulting MNN can be further used to accelerate the evaluation of the mapping, especially when a large number of evaluations are needed within a certain range of parameters.

Acknowledgements

The authors thank Yuehaw Khoo for constructive discussions. The work of Y.F. and L.Y. is partially supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Scientific Discovery through Advanced Computing (SciDAC) program and the National Science Foundation under award DMS-1818449, and the GCP Research Credits Program from Google. The work of Y.F. is also partially supported by AWS Cloud Credits for Research program from Amazon. The work of L.L and L. Z. is partially supported by the Department of Energy under Grant No. DE-SC0017867, by the Department of Energy under the CAMERA project, and by the Air Force Office of Scientific Research under award number FA9550-18-1-0095.

References