Machine learning has been gaining recognition rapidly as a powerful computational technique to address some of the most challenging problems arising from scientific and engineering computations (SEC) with promising results in simulations of biological and quantum systems, fluid dynamics, wave scattering, high dimensional PDEs, and inverse problems, etc. This special issue contains 1 survey paper and 17 original research articles on recent developments in machine learning, especially deep neural networks, concerning both its theoretical and algorithmic aspects pertinent to SEC.

The dying ReLU refers to the problem when ReLU neurons become inactive and only output 0 for any input. There are many empirical and heuristic explanations of why ReLU neurons die. However, little is known about its theoretical analysis. In this paper, we rigorously prove that a deep ReLU network will eventually die in probability as the depth goes to infinite. Several methods have been proposed to alleviate the dying ReLU. Perhaps, one of the simplest treatments is to modify the initialization procedure. One common way of initializing weights and biases uses symmetric probability distributions, which suffers from the dying ReLU. We thus propose a new initialization procedure, namely, a randomized asymmetric initialization. We show that the new initialization can effectively prevent the dying ReLU. All parameters required for the new initialization are theoretically designed. Numerical examples are provided to demonstrate the effectiveness of the new initialization procedure.

We study a family of $H^m$-conforming piecewise polynomials based on the artificial neural network, referred to as the finite neuron method (FNM), for numerical solution of $2m$-th-order partial differential equations in $\mathbb{R}^d$ for any $m,d≥1$ and then provide convergence analysis for this method. Given a general domain Ω$⊂\mathbb{R}^d$ and a partition $\mathcal{T}_h$ of Ω, it is still an open problem in general how to construct a conforming finite element subspace of $H^m$(Ω) that has adequate approximation properties. By using techniques from artificial neural networks, we construct a family of $H^m$-conforming functions consisting of piecewise polynomials of degree $k$ for any $k≥m$ and we further obtain the error estimate when they are applied to solve the elliptic boundary value problem of any order in any dimension. For example, the error estimates that $‖u−u_N‖_{H^m(\rm{Ω})}=\mathcal{O}(N^{−\frac{1}{2}−\frac{1}{d}})$ is obtained for the error between the exact solution $u$ and the finite neuron approximation $u_N$. A discussion is also provided on the difference and relationship between the finite neuron method and finite element methods (FEM). For example, for the finite neuron method, the underlying finite element grids are not given a priori and the discrete solution can be obtained by only solving a non-linear and non-convex optimization problem. Despite the many desirable theoretical properties of the finite neuron method analyzed in the paper, its practical value requires further investigation as the aforementioned underlying non-linear and non-convex optimization problem can be expensive and challenging to solve. For completeness and the convenience of the reader, some basic known results and their proofs are introduced.

We study the training process of Deep Neural Networks (DNNs) from the Fourier analysis perspective. We demonstrate a very universal Frequency Principle (F-Principle) — DNNs often fit target functions from low to high frequencies — on high-dimensional benchmark datasets such as MNIST/CIFAR10 and deep neural networks such as VGG16. This F-Principle of DNNs is opposite to the behavior of Jacobi method, a conventional iterative numerical scheme, which exhibits faster convergence for higher frequencies for various scientific computing problems. With theories under an idealized setting, we illustrate that this F-Principle results from the smoothness/regularity of the commonly used activation functions. The F-Principle implies an implicit bias that DNNs tend to fit training data by a low-frequency function. This understanding provides an explanation of good generalization of DNNs on most real datasets and bad generalization of DNNs on parity function or a randomized dataset.

This paper quantitatively characterizes the approximation power of deep feed-forward neural networks (FNNs) in terms of the number of neurons. It is shown by construction that ReLU FNNs with width $\mathcal{O}$(max{$d⌊N^{1/d}⌋$,$N$+1}) and depth $\mathcal{O}(L)$ can approximate an arbitrary Hölder continuous function of order $α∈(0,1]$ on $[0,1]^d$ with a nearly tight approximation rate $\mathcal{O}(\sqrt{d}N^{−2α/d}L^{−2α/d})$ measured in $L^p$ -norm for any $N,L∈\mathbb{N}^+$ and $p∈[1,∞]$. More generally for an arbitrary continuous function $f$ on $[0,1]^d$ with a modulus of continuity $ω_f (·)$, the constructive approximation rate is $\mathcal{O}(\sqrt{d}ω_f(N^{−2α/d}L^{−2α/d}))$. We also extend our analysis to $f$ on irregular domains or those localized in an ε-neighborhood of a $d_\mathcal{M}$-dimensional smooth manifold $\mathcal{M}⊆[0,1]^d$ with $d_\mathcal{M}≪d$. Especially, in the case of an essentially low-dimensional domain, we show an approximation rate $\mathcal{O}(ω_f(\frac{ε}{1−δ}\sqrt{\frac{d}{d_δ}}+ε)+\sqrt{d}ω_f(\frac{\sqrt{d}}{1−δ\sqrt{d_δ}}N^{−2α/d_δ}L^{−2α/d_δ})$ for ReLU FNNs to approximate $f$ in the $ε$-neighborhood, where $d_δ=\mathcal{O}(d_\mathcal{M}\frac{\rm{ln}(d/δ)}{δ^2})$ for any $δ∈(0,1)$ as a relative error for a projection to approximate an isometry when projecting $\mathcal{M}$ to a $d_δ$-dimensional domain.

We propose a data fusion method based on multi-fidelity Gaussian process regression (GPR) framework. This method combines available data of the quantity of interest (QoI) and its gradients with different fidelity levels, namely, it is a Gradient-enhanced Cokriging method (GE-Cokriging). It provides the approximations of both the QoI and its gradients simultaneously with uncertainty estimates. We compare this method with the conventional multi-fidelity Cokriging method that does not use gradients information, and the result suggests that GE-Cokriging has a better performance in predicting both QoI and its gradients. Moreover, GE-Cokriging even shows better generalization result in some cases where Cokriging performs poorly due to the singularity of the covariance matrix. We demonstrate the application of GE-Cokriging in several practical cases including reconstructing the trajectories and velocity of an underdamped oscillator with respect to time simultaneously, and investigating the sensitivity of power factor of a load bus with respect to varying power inputs of a generator bus in a large scale power system. Although GE-Cokriging requires slightly higher computational cost than Cokriging in some cases, the comparison of the accuracy shows that this cost is worthwhile.

Deep networks, especially convolutional neural networks (CNNs), have been successfully applied in various areas of machine learning as well as to challenging problems in other scientific and engineering fields. This paper introduces Butterfly-net, a low-complexity CNN with structured and sparse cross-channel connections, together with a Butterfly initialization strategy for a family of networks. Theoretical analysis of the approximation power of Butterfly-net to the Fourier representation of input data shows that the error decays exponentially as the depth increases. Combining Butterfly-net with a fully connected neural network, a large class of problems are proved to be well approximated with network complexity depending on the effective frequency bandwidth instead of the input dimension. Regular CNN is covered as a special case in our analysis. Numerical experiments validate the analytical results on the approximation of Fourier kernels and energy functionals of Poisson's equations. Moreover, all experiments support that training from Butterfly initialization outperforms training from random initialization. Also, adding the remaining cross-channel connections, although significantly increases the parameter number, does not much improve the post-training accuracy and is more sensitive to data distribution.

Algorithms based on deep neural networks (DNNs) have attracted increasing attention from the scientific computing community. DNN based algorithms are easy to implement, natural for nonlinear problems, and have shown great potential to overcome the curse of dimensionality. In this work, we utilize the multi-scale DNN-based algorithm (MscaleDNN) proposed by Liu, Cai and Xu (2020) to solve multi-scale elliptic problems with possible nonlinearity, for example, the p-Laplacian problem. We improve the MscaleDNN algorithm by a smooth and localized activation function. Several numerical examples of multi-scale elliptic problems with separable or non-separable scales in low-dimensional and high-dimensional Euclidean spaces are used to demonstrate the effectiveness and accuracy of the MscaleDNN numerical scheme.

Random batch algorithms are constructed for quantum Monte Carlo simulations. The main objective is to alleviate the computational cost associated with the calculations of two-body interactions, including the pairwise interactions in the potential energy, and the two-body terms in the Jastrow factor. In the framework of variational Monte Carlo methods, the random batch algorithm is constructed based on the over-damped Langevin dynamics, so that updating the position of each particle in an $N$-particle system only requires $\mathcal{O}(1)$ operations, thus for each time step the computational cost for $N$ particles is reduced from $\mathcal{O}(N^2)$ to $\mathcal{O}(N)$. For diffusion Monte Carlo methods, the random batch algorithm uses an energy decomposition to avoid the computation of the total energy in the branching step. The effectiveness of the random batch method is demonstrated using a system of liquid $^4$He atoms interacting with a graphite surface.

In scientific and engineering applications, often sufficient low-cost low-fidelity data is available while only a small fractional of high-fidelity data is accessible. The multi-fidelity model integrates a large set of low-cost but biased low-fidelity datasets with a small set of precise but high-cost high-fidelity data to make an accurate inference of quantities of interest. Under many circumstances, the number of model input dimensions is often high in real applications. To simplify the model, dimension reduction is often used. The gradient-free active subspace is employed in this research for dimension reduction. In this work, a novel predictive model for high-dimensional nonlinear problems by integrating the nonlinear multi-fidelity Gaussian progress regression and the gradient-free active subspace method is put forward. Numerical results demonstrated that the proposed approach can not only perform effective dimension reduction on the original data but also obtain accurate prediction results thanks to the effective dimension reduction procedure.

In this paper, we propose multi-scale deep neural networks (MscaleDNNs) using the idea of radial scaling in frequency domain and activation functions with compact support. The radial scaling converts the problem of approximation of high frequency contents of PDEs' solutions to a problem of learning about lower frequency functions, and the compact support activation functions facilitate the separation of frequency contents of the target function to be approximated by corresponding DNNs. As a result, the MscaleDNNs achieve fast uniform convergence over multiple scales. The proposed MscaleDNNs are shown to be superior to traditional fully connected DNNs and be an effective mesh-less numerical method for Poisson-Boltzmann equations with ample frequency contents over complex and singular domains.

Physics informed neural networks (PINNs) are deep learning based techniques for solving partial differential equations (PDEs) encountered in computational science and engineering. Guided by data and physical laws, PINNs find a neural network that approximates the solution to a system of PDEs. Such a neural network is obtained by minimizing a loss function in which any prior knowledge of PDEs and data are encoded. Despite its remarkable empirical success in one, two or three dimensional problems, there is little theoretical justification for PINNs.
As the number of data grows, PINNs generate a sequence of minimizers which correspond to a sequence of neural networks. We want to answer the question: Does the sequence of minimizers converge to the solution to the PDE? We consider two classes of PDEs: linear second-order elliptic and parabolic. By adapting the Schauder approach and the maximum principle, we show that the sequence of minimizers strongly converges to the PDE solution in $C^0$. Furthermore, we show that if each minimizer satisfies the initial/boundary conditions, the convergence mode becomes $H^1$. Computational examples are provided to illustrate our theoretical findings. To the best of our knowledge, this is the first theoretical work that shows the consistency of PINNs.

We propose a universal discontinuity detector using convolution neural network (CNN) and apply it in conjunction of solving nonlinear conservation laws in both 1D and 2D. The CNN detector is trained offline with synthetic data. The training data are generated using randomly constructed piecewise functions, which are then processed using randomized linear advection solver to count for the cases of numerical errors in practice. The detector is then paired with high-order numerical solvers. In particular, we combined high-order WENO in troubled cells with high-order central difference in smooth region. Extensive numerical examples are presented. We observe that the proposed method produces notably sharper and cleaner signals near the discontinuities, when compared to other well known troubled cell detector methods.

One of the most challenging issues in applied mathematics is to develop and analyze algorithms which are able to approximately compute solutions of high-dimensional nonlinear partial differential equations (PDEs). In particular, it is very hard to develop approximation algorithms which do not suffer under the curse of dimensionality in the sense that the number of computational operations needed by the algorithm to compute an approximation of accuracy $ε$>0 grows at most polynomially in both the reciprocal 1/$ε$ of the required accuracy and the dimension $d∈\mathbb{N}$ of the PDE. Recently, a new approximation method, the so-called full history recursive multilevel Picard (MLP) approximation method, has been introduced and, until today, this approximation scheme is the only approximation method in the scientific literature which has been proven to overcome the curse of dimensionality in the numerical approximation of semilinear PDEs with general time horizons. It is a key contribution of this article to extend the MLP approximation method to systems of semilinear PDEs and to numerically test it on several example PDEs. More specifically, we apply the proposed MLP approximation method in the case of Allen-Cahn PDEs, Sine-Gordon-type PDEs, systems of coupled semilinear heat PDEs, and semilinear Black-Scholes PDEs in up to 1000 dimensions. We also compare the performance of the proposed MLP approximation algorithm with a deep learning based approximation method from the scientific literature.

In this paper, we study a multi-scale deep neural network (MscaleDNN) as a meshless numerical method for computing oscillatory Stokes flows in complex domains. The MscaleDNN employs a multi-scale structure in the design of its DNN using radial scalings to convert the approximation of high frequency components of the highly oscillatory Stokes solution to one of lower frequencies. The MscaleDNN solution to the Stokes problem is obtained by minimizing a loss function in terms of $L^2$ norm of the residual of the Stokes equation. Three forms of loss functions are investigated based on vorticity-velocity-pressure, velocity-stress-pressure, and velocity-gradient of velocity-pressure formulations of the Stokes equation. We first conduct a systematic study of the MscaleDNN methods with various loss functions on the Kovasznay flow in comparison with normal fully connected DNNs. Then, Stokes flows with highly oscillatory solutions in a 2-D domain with six randomly placed holes are simulated by the MscaleDNN. The results show that MscaleDNN has faster convergence and consistent error decays in the simulation of Kovasznay flow for all three tested loss functions. More importantly, the MscaleDNN is capable of learning highly oscillatory solutions when the normal DNNs fail to converge.

Conservation laws are considered to be fundamental laws of nature. It has broad applications in many fields, including physics, chemistry, biology, geology, and engineering. Solving the differential equations associated with conservation laws is a major branch in computational mathematics. The recent success of machine learning, especially deep learning in areas such as computer vision and natural language processing, has attracted a lot of attention from the community of computational mathematics and inspired many intriguing works in combining machine learning with traditional methods. In this paper, we are the first to view numerical PDE solvers as an MDP and to use (deep) RL to learn new solvers. As proof of concept, we focus on 1-dimensional scalar conservation laws. We deploy the machinery of deep reinforcement learning to train a policy network that can decide on how the numerical solutions should be approximated in a sequential and spatial-temporal adaptive manner. We will show that the problem of solving conservation laws can be naturally viewed as a sequential decision-making process, and the numerical schemes learned in such a way can easily enforce long-term accuracy. Furthermore, the learned policy network is carefully designed to determine a good local discrete approximation based on the current state of the solution, which essentially makes the proposed method a meta-learning approach. In other words, the proposed method is capable of learning how to discretize for a given situation mimicking human experts. Finally, we will provide details on how the policy network is trained, how well it performs compared with some state-of-the-art numerical solvers such as WENO schemes, and supervised learning based approach L3D and PINN, and how well it generalizes.

In Bayesian inverse problems, surrogate models are often constructed to speed up the computational procedure, as the parameter-to-data map can be very expensive to evaluate. However, due to the curse of dimensionality and the nonlinear concentration of the posterior, traditional surrogate approaches (such us the polynomial-based surrogates) are still not feasible for large scale problems. To this end, we present in this work an adaptive multi-fidelity surrogate modeling framework based on deep neural networks (DNNs), motivated by the facts that the DNNs can potentially handle functions with limited regularity and are powerful tools for high dimensional approximations. More precisely, we first construct offline a DNN-based surrogate according to the prior distribution, and then, this prior-based DNN-surrogate will be adaptively & locally refined online using only a few high-fidelity simulations. In particular, in the refine procedure, we construct a new shallow neural network that views the previous constructed surrogate as an input variable – yielding a composite multi-fidelity neural network approach. This makes the online computational procedure rather efficient. Numerical examples are presented to confirm that the proposed approach can obtain accurate posterior information with a limited number of forward simulations.

Neural network-based machine learning is capable of approximating functions in very high dimension with unprecedented efficiency and accuracy. This has opened up many exciting new possibilities, not just in traditional areas of artificial intelligence, but also in scientific computing and computational science. At the same time, machine learning has also acquired the reputation of being a set of "black box" type of tricks, without fundamental principles. This has been a real obstacle for making further progress in machine learning.
In this article, we try to address the following two very important questions: (1) How machine learning has already impacted and will further impact computational mathematics, scientific computing and computational science? (2) How computational mathematics, particularly numerical analysis, can impact machine learning? We describe some of the most important progress that has been made on these issues. Our hope is to put things into a perspective that will help to integrate machine learning with computational mathematics.

We propose a generalized space-time domain decomposition approach for the physics-informed neural networks (PINNs) to solve nonlinear partial differential equations (PDEs) on arbitrary complex-geometry domains. The proposed framework, named eXtended PINNs ($XPINNs$), further pushes the boundaries of both PINNs as well as conservative PINNs (cPINNs), which is a recently proposed domain decomposition approach in the PINN framework tailored to conservation laws. Compared to PINN, the XPINN method has large representation and parallelization capacity due to the inherent property of deployment of multiple neural networks in the smaller subdomains. Unlike cPINN, XPINN can be extended to any type of PDEs. Moreover, the domain can be decomposed in any arbitrary way (in space and time), which is not possible in cPINN. Thus, XPINN offers both space and time parallelization, thereby reducing the training cost more effectively. In each subdomain, a separate neural network is employed with optimally selected hyperparameters, e.g., depth/width of the network, number and location of residual points, activation function, optimization method, etc. A deep network can be employed in a subdomain with complex solution, whereas a shallow neural network can be used in a subdomain with relatively simple and smooth solutions. We demonstrate the versatility of XPINN by solving both forward and inverse PDE problems, ranging from one-dimensional to three-dimensional problems, from time-dependent to time-independent problems, and from continuous to discontinuous problems, which clearly shows that the XPINN method is promising in many practical problems. The proposed XPINN method is the generalization of PINN and cPINN methods, both in terms of applicability as well as domain decomposition approach, which efficiently lends itself to parallelized computation. The XPINN code is available on $https://github.com/AmeyaJagtap/XPINNs$.