solutions of PDEs [31, 88] to complex mathematical modeling, prediction, and classification tasks in physics [13], biology [76,83,93], and engineering [59,82].
Simultaneously, the broader applied mathematics community has taken interest in the approximation capabilities of NNs [6,10,69,90]. The earliest results in this direction [19,51] established that even a single hidden layer fully-connected NN has universal approximation capability: so long as the number of nodes in the hidden layer are allowed to grow unbounded, such architectures are able to approximate any Borel measurable function on a compact domain to arbitrary uniform accuracy. More recent works have studied the connection between expressiveness of DNNs and their depth [57, 60, 91], while others have established connections between DNNs and other methods of approximation, e.g., sparse grids [62], splines [85], polynomials [22,72], and “h, p”-finite elements [67]. A plethora of results now exist concerning the approximation power of DNNs for different function spaces – e.g. Sobolev spaces [45], bandlimited functions [62], analytic functions [33,68], Barron functions [32], cartoon-like functions [44], H¨older spaces [73] – and tasks in scientific computing, such as approximation of high-dimensional functions [56, 72] and PDEs [11, 43, 55], dimensionality reduction [92], and methods for DEs [61, 88]. Theoretically, these works establish best-in-class approximation properties of DNNs for many problems. Concurrently, some works have sought to address the construction of DNN approximations directly, although typically without theoretical guarantees on trainability. For example [26, 36] employ ideas from greedy methods in order to do this, [27] derives a formula for integral representations of shallow ReLU networks, and [23] construct networks directly approximating the Legendre basis, using this as an initialization point for training. On the other hand, works such as [20] reinterpret DNN approximation as approximation in adaptive bases, which can be learned via training data.
1.1. Challenges. Yet despite the impressive empirical and theoretical results achieved in the broader DL community, there is concern that methods based on DNNs do not currently meet the usual rigorous standards for algorithms in computational science [9]. While the aforementioned theoretical results assert the expressibility of the class of DNNs – that is, the existence of a DNN of a given architecture that achieves a desired rate of convergence for a given problem – they say little about their practical performance when trained by modern approaches in DL. If such techniques are to achieve widespread adoption in scientific computing, it is vital they be understood through the lens of numerical analysis, namely, (i) stability, (ii) accuracy, (iii) sample complexity, (iv) curse of dimensionality and (v) computational cost.
(i) Stability. Recently, researchers have begun to question the stability properties of DNNs [35, 63, 81]. A series of works have demonstrated that DNNs trained on tasks such as image classification are vulnerable to misclassification when provided images with small “adversarial” perturbations [64] and can even completely fail on image reconstruction tasks in the presence of small structural changes in the data [4, 42]. As deep learning is increasingly being applied towards critical problems in healthcare, e.g., DeepMind’s recent work on machine-assisted diagnostic imaging in retinal disease [24], many have questioned the ethics of applying tools whose stability properties are not fully understood to such problems.
(ii) Accuracy. Over the past 5 years, many works have been published on the classes of functions, e.g., analytic or piece-wise continuous, that can be approximated by DNNs of a given size with a certain rate of convergence. These results are constructive, often showing
the existence of a DNN emulating another approximation scheme, e.g. polynomials, for which convergence rates have already been established. While such results provide a useful benchmark for DNN expressivity, they do not suggest methods for training DNNs that reliably achieve the tolerances required in computational science applications.
(iii) Sample complexity. Areas in which DL has seen the greatest success include problems in supervised learning such as image classification. In such settings, DNNs are trained on large sets of labeled images, yielding a model capable of predicting labels for unseen images. Popular datasets for DL competitions include the ImageNet database which contains 14 million hand-annotated images of more than 20,000 categories of subjects [25]. In contrast, problems in computational science are often relatively data-starved, e.g. applications in uncertainty quantification (UQ) which involve computing a quantity of interest from sampled solutions of a parameterized PDE [46]. As each sample involves the discretization and solution of a PDE, which may require thousands of degrees of freedom to accurately resolve, there is great attention paid in such problems to minimizing the required number of samples [3,29].
(iv) Curse of dimensionality. Many modern problems in scientific computing involve high dimensionality. High-dimensional PDEs occur in numerous applications, and parametrized PDEs in UQ applications often involve tens to hundreds of variables. Recent works have shown that certain DNNs have the expressive capabilities to mitigate the curse of dimensionality to the same extent as current best-in-class schemes [11,33,43,55,62,68,72]. Yet, as noted, this does not assert these rates can be achieved via training. Moreover, the curse of dimensionality is an important consideration in the sample complexity, as the cost of obtaining samples can often dominate the overall cost. A recent numerical study has shown that approximation quality degrades with increasing dimension for approximating solutions of high-dimesional parameterized PDEs [38]. However the observed scaling is not exponential in the dimension d, but dependent on the complexity of the underlying problem. Understanding this scaling with respect to the sample complexity is crucial to applying these methods in computational science applications.
(v) Computational cost. By far, the largest barrier to entry for DL research is the cost of training. DNNs are typically trained on graphics processing units (GPUs), and a single GPU can cost thousands of US dollars. In many industry applications, models are trained on hundreds of these specialized cards. In addition, the training process itself is very energyintensive, and can produce a large amount of excess CO. Even a small reduction in computational cost can yield large cost savings and greater access to resources for researchers.
In the near term, it seems likely that any DL implementation will pay a price in computational cost. Hence there needs to be a clear understanding of the benefits vis-a-vis properties (i)–(iv) above. The study of these concerns is the broad purpose of this paper.
1.2. Contributions. Our main objective is to examine practical DNN approximation on problems motivated by scientific computing. In many applications in computational science, the core task involves approximating a function , with domain
(often
1). Hence our main aim is to examine the performance of DL on practical function
approximation through the five considerations (i)–(v). Our main contributions are:
1. We develop a computational framework for examining the practical capabilities of DNNs in scientific computing, based on the rigorous testing principles of numerical analysis. We provide clear practical guidance on training DNNs for function approximation problems.
2. We conduct perhaps the first comprehensive empirical study of the performance of training fully-connected feedforward ReLU DNNs on standard function classes considered in numerical analysis, namely, (piecewise) smooth functions on bounded domains. We compare performance over a range of dimensions, examining the capability of DL for mitigating the curse of dimensionality. We examine the effect of network architecture (depth and width) on both ease of training and approximation performance of the trained DNNs. We also make a clear empirical comparison between DL and current best-in-class approximation schemes for smooth function approximation. The latter is based on polynomial approximation via compressed sensing (CS) [3], which (as we also show) achieves exponential rates of convergence for analytic functions in arbitrarily-many dimensions, see Section 5.1 for more details.
3. We present theoretical analysis that compares the performance of DL to that of CS for smooth function approximation. In particular, we establish a novel practical existence theorem, which asserts the existence of a DNN architecture and training procedure (based on minimizing a certain cost function) that attains the same exponential rate of convergence as CS for analytic function approximation with the same sample complexity.
1.3. Conclusions. The primary conclusion of this work is the following. While it is in- creasingly well understood that DNNs have substantial expressive power for problems relevant to scientific computing, there remains a large gap between expressivity and practical performance achieved with standard methods of training. Surprisingly, trained DNNs can perform very badly on functions for which there are strong expressivity results, such as smooth functions in high dimensions and piecewise smooth functions. Yet, on other examples, they are competitive with current best-in-class schemes based on CS. We also draw the following conclusions based on our experimental results training fully-connected feedforward ReLU DNNs:
1. The accuracy of trained DNNs is limited by the precision used, but is typically nowhere near machine epsilon despite training to such tolerances in the loss. In this work, we perform experiments in both single and double precision. Yet, in both cases, it is typically impossible to get beyond four digits of accuracy even when approximating extremely smooth functions. In contrast, a combination of Legendre polynomial approximation in the hyperbolic cross subspace with weighted -minimization and a specific choice of weights motivated by smooth function approximation can reliably achieve six or seven digits of accuracy on such problems. Since our theoretical contribution shows that DNNs should be capable of obtaining these results (albeit with a different initialization and training procedure), this fact suggests a limitation of standard training methods in obtaining more accurate results.
2. The training process itself is also highly sensitive to the parameterization of the solvers and initialization of the weights and biases of the DNNs. After extensive testing, we chose the Adam optimizer with exponentially decaying learning rate over a variety of optimizers including standard SGD, finding empirically that this strategy can help to mitigate some of the challenges of solving the non-convex optimization problem of training, e.g., non-monotonic decrease in the loss, and slow convergence to minimizers. We also initialize our networks with
symmetric uniform or normal distributions with small variance, finding larger variances can lead to failure in training or slow convergence. These choices combined lead to a training process that is largely stable and minimizes the probability of failure in training over a wide range of architectures. However, failures can still occur and are often a consequence of choosing a network that is either too shallow and narrow (where failure occurs immediately and the error stagnates), or too wide and deep (where the resulting network achieves order machine epsilon tolerance in the loss, but massively overfits exhibiting numerical artifacts). 3. Generally speaking, deeper architectures (which are sufficiently wide) are both easier
to train and perform better than shallower architectures. However, this is clearly highly-
dependent on the regularity and dimensionality of the target function for smooth-function approximation problems. The width of the network also plays an important role; we find networks with width 5-10 times larger than the depth perform better. While this trend appears to be general, some of our results on less-smooth and piecewise continuous functions indicate that it is not universal. 4. While much of the success of DL has been in the field of classification, surprisingly
performance of trained DNNs on simple piecewise constant functions is relatively poor in
comparison to smooth functions, and adversely affected by the curse of dimensionality. Yet, DNNs do approximate such functions to some accuracy, unlike polynomial-based CS tech-
niques. This highlights the flexibility of the DL approach.
We also draw the following theoretical conclusions: 5. For analytic function approximation, a certain DNN strategy based on emulating poly-
nomials via a suitable DNN can achieve the same guaranteed performance (up to constants)
as best-in-class schemes based on CS. 6. However, such scheme will not typically offer any superior performance. In particular,
it is not flexible. It only perform well on the function classes on which CS performs well, and
perform correspondingly poorly on other, e.g. piecewise smooth functions. This is in contrast to the experimental results in this paper, which show that standard feedforward architectures and training via the -loss can approximate both smooth and nonsmooth functions, and expressibility results, which show that DNN architectures can approximate functions from a range of different function classes.
1.4. Outlook. This paper raises and seeks to answer the following question: is DL a useful tool for problems in scientific computing? Some of the above conclusions may appear rather negative in this regard, certainly in comparison to the positive impression given by the plethora of expressivity results on DNNs. Let us raise several caveats. First, this study considers one particular setup: namely, fully-connected, ReLU networks of constant hidden layer widths, in combination with the -loss function. There are almost countless variations on this setup, some of which will undoubtedly perform better. These variations include different activation functions (e.g. sigmoid, hyperbolic tangent), different architectures (e.g. ResNets, sparsely-connected layers, convolutional layers) and different loss functions (e.g. those incorporating regularization). We elected to use this setup based on standards in the literature; for instance, most expressibility results consider ReLU activations. Given difficulties and intense computational resources required for training, it is beyond the scope of this first work to methodically compare all possible setups. Second, the practical existence theorem we prove
shows that a careful choice of architecture and cost function can allow DNNs to offer similar performance to state-of-the-art techniques. While it comes at the price of flexibility, this result nonetheless provides a theoretical benchmark for DNN approximations. It also demonstrates the potential of DNN strategies for eventually outperforming such methods, and indicates that theory-inspired architecture and training strategy design as a way to achieve such improvements. The extent to which this can (provably) be done while maintaining flexibility – i.e. the ability to approximate different function classes with the same DNN procedure – is an enticing challenge for future work.
1.5. Outline. The outline of the remainder of this paper is as follows. In we introduce the approximation problem, DNNs, DL and CS. In
we describe the experimental setup, including details of the training procedure used. Our numerical results are found in
in
we present our theoretical results. Additional information for this paper is contained in the Supplementary Materials. Code accompanying the computational framework is available at https://github.com/ndexter/MLFA.
2. Framework. In this section, we first describe the function approximation problem, and then introduce DNNs, DL, polynomial approximation and CS.
2.1. Problem formulation. Throughout this paper, we consider the unit cube in d dimensions, , equipped with the uniform probability measure d
is the Lebesgue measure and
-dimensional variable. Let
the space of real-valued square-integrable functions on U with respect to
. Our objective is to approximate an unknown function
) from samples. These samples are generated by simple Monte Carlo sampling: we draw
randomly and independently from the measure
. Hence the approximation problem we aim to solve is
(AP) Given the measurements , approximate f.
We note in passing that much of what follows in this paper can be extended to more general domains, sampling strategies and to functions taking values in other vector spaces (for instance, complex-valued functions, vector-valued functions, or even Hilbert-valued functions, as arise commonly in UQ applications [29]). We assume the above setup for ease of presentation.
We require several further pieces of notation. We write -norm with respect to
. The space of essentially bounded functions on U is denoted by
) and its norm by
) to denote a (multi)index of length
is a finite or countable (multi)index set, we write
) for the space of
sequences
, i.e. those satisfying
we define
in the usual way.
2.2. Deep Learning. We now introduce DL. First, we recall the definition of a DNN:
with affine linear maps , and the activation function
acting component-wise (i.e.,
) is called a neural network (NN). The map
corresponding to layer l is given by
th weight matrix and
th bias vector. We refer to L as the depth of the network and max
as its width.
Informally, we consider a DNN as any NN with 1 hidden layers. Definition 2.1 pertains to feedforward DNNs. We do not consider more exotic constructions such as recurrent networks or ResNets in this paper. We also consider so-called fully connected networks, meaning that the weights and biases can take arbitrary real values. The layers l = 1, . . . , L are referred to as hidden layers. In our experiments later, we set their widths to be equal,
. Note that
are specified by the problem. In our case,
and
There are numerous choices for the activation function , and moreover, one may also choose different activation functions in different layers. Since it is popular both in theory and in practice, we use the rectified linear unit (ReLU), defined by
also refer to a DNN architecture having ReLU activation function and L hidden layers with N nodes per layer as a ReLU
The architecture of a network is the specific choice of activation and parameters L and
. We denote the set of neural networks of a given architecture by N. Note that this family is parametrized by the weight matrices and biases. Selecting the right architecture for a given problem is a significant challenge. We discuss this topic further in
Given an unknown function is the process of computing a neural network Φ that approximates f from the data
. This is normally achieved by minimizing a
, i.e. we solve
where N is the family of neural networks of the chosen architecture. Note that this is equivalent to a minimization problem for the weights and biases
. A typical choice is the
-loss (also known as empirical risk, mean squared loss):
We primarily use this loss function in this paper. However, many other choices are possible. For instance, it is common to add a regularization term to the loss function, e.g.
Here is chosen to promotes some desirable features of the network. For instance, J may be a norm of the weight matrices, thus promoting small and/or sparse weights.
2.3. Polynomial approximation of smooth functions. We now introduce the polynomial approximation schemes against which we compare DL for function approximation. Polynomial
approximation is a vast and classical topic. Yet, it has received renewed attention in the last several decades, motivated by applications in UQ where one seeks to approximate a smooth quantity of interest of a parametric PDE [18]. The particular scheme we consider is based on orthogonal expansions in orthonormal polynomials in ), i.e. multivariate Legendre polynomials. The univariate, orthonormal Legendre polynomials on the [
1] are defined by
where is the classical Legendre polynomial with normalization
form an orthonormal basis of
1, we define the tensor orthonormal Legendre polynomials as
The set forms an orthonormal basis of
). Hence any function
convergent expansion
(2.3)
where ) is the coefficient of f with respect to Ψ
. Note that the sequence
is an element of
), the space of square-summable sequences with indices in
. By Parseval’s identity,
When d = 1, approximating a smooth function in the Legendre basis is typically achieved by truncating the expansion (2.3) after its first s terms, then using, for instance, least-squares to approximately recover the coefficients from the measurements
In
2 dimensions, the situation becomes more complicated, since there are many different choices of index set S of cardinality s one might employ to truncate the expansion (2.3):
By Parseval’s identity, the error
on the coefficients outside S. Hence, an a priori choice of S may have limited effectiveness, since it may fail to capture any anisotropic behaviour of f. This naturally motivates the concept of best s-term approximation [18]. In best s-term approximation, the index set S is chosen so that it contains the multi-indices corresponding to the large s coefficients
absolute value. This is a type of nonlinear approximation scheme [28]. If ˜
denotes the best s-term approximation, the error satisfies
Under appropriate conditions (e.g. f is analytic), this approximation converges exponentially fast in s (see Theorem 5.2). In high dimensions this significantly improves over any linear approximation scheme based on a fixed, isotropic choice of S. We discuss this further in
2.4. Polynomial approximation with compressed sensing. Computing the best s-term approximation is on the face of it a daunting task. In theory, it involves computing all infinitely-many of the coefficients c, then selecting the largest s. This is of course intractable, and generally still computationally infeasible even if one limits oneself to computing a large, but finite number of coefficients.
A solution is to use compressed sensing. Here one first selects a large, but finite multiindex set Λ . This set is generally assumed to contain the coefficients of some quasi-best s-term approximation, if not the coefficients of the true best s-term approximation itself. For reasons that will be made clear in
(2.4) Λ
is the hyperbolic cross index set of degree s. Having chosen Λ, the finite vector can now be assumed to be approximately sparse. Next, one formulates the normalized measurement matrix and vector of measurements
Then one searches for an approximately sparse solution of the linear system Az = f. A standard means to do this is to solve the quadratically-constrained basis pursuit problem
for suitably chosen
minimize
for appropriately chosen 0. A solution ˆ
of either problem yields an approximation ˆ
Unfortunately, simply promoting the sparsity of the polynomial coefficients via the
norm is not sufficient to achieve favourable sample complexity bounds. Bounds on
norm based approaches can be exponential in the dimension d [3]. Fortunately, as considered in [1, 15, 70], this issue can be overcome by replacing the
-norm with a certain weighted
-norm. For instance, instead of (2.6) one now solves
minimize
Here is a vector of weights and
. As shown in [1], an appropriate choice of weights is
(2.9)
For the remainder of this paper, we consider the CS polynomial approximation scheme ˆf = is either a solution of (2.8), or for the sake of comparison, (2.6).
3. Testing setup. We now describe the testing setup for our experiments. Training DNNs requires careful choices of the optimization solver, initialization and optimization parameters. This section describes the choices we made to deliver consistent performance across a range of DNN architectures. In summary, we find the Adam (adaptive moments) optimizer with exponentially-decaying learning rate and a specific random initialization to deliver this performance, whereas other solvers such as SGD and other learning rate schedules perform less consistently.
3.1. Setup. We first summarize the main methodology:
(i) Implementation. Our framework has been implemented in a package called MLFA (Machine Learning Function Approximation) in version 1.13 of Google’s TensorFlow software library https://www.tensorflow.org/, and is available on GitHub at https://github.com/ ndexter/MLFA. Details about the set of features supported by MLFA and data recorded by the code can also be found on the GitHub page.
(ii) Hardware. In the course of testing our DNN models, we observed improved accuracy on some of our test problems by initializing and training the networks in double precision. Modern GPUs often support half, single, and double precision arithmetic, though many commonly-available GPUs are optimized to perform single precision computations much more quickly (see ). The majority of our computations were performed in single precision using the Tesla P100 GPUs on Compute Canada’s Cedar compute cluster at Simon Fraser University
, though for some of our test problems we provide double precision results for a subset of the architectures considered for comparison.
(iii) Choice of architectures and initialization. We consider fully-connected ReLU networks. This choice is inspired by the many theoretical existence results on such networks (see ). There is a vast literature suggesting various strategies for designing architectures and initializing neural networks. For an introduction to these topics see, e.g., [41,
5.2 & 8.4]. We recall the work [47] focusing on which DNN architectures result in exploding and vanishing gradients, a common problem in backpropagation which can result in failure during training. Our empirical results confirm that choosing DNN architectures with a fixed number of nodes per layer N and depth L such that the ratio
is small, e.g.,
effective choice for training. In
, we study several popular strategies for initializing DNNs. Our results show that initializing the weights and biases to be normal random variables with mean 0 and variance 0.01 is an effective choice for many of the architectures and problems considered herein. We note that for the range of architectures studied, this choice results in weights and biases with variance smaller than many other popular choices, e.g., the strategies from [40,49], and is sufficiently small to avoid the failure modes analyzed in [48]. We also note that setting the variances smaller than 0.01 can result in networks which are more difficult to train, which can be attributed to the exponentially decaying length scales described in [48, Theorem 1]. The seed used for the random number generators can also have an important effect on the initialization. In this study, we initialize all of our networks from the same seed 0 for both TensorFlow and NumPy, in order to reduce the complexity of our experiments. A more comprehensive study of the average performance of DL would require averaging over a large set of seeds used in initialization. We leave such a study to a future work.
compares the performance of a variety of solvers and learning rate schedules on the function approximation problem. In practice, we find the Adam optimizer [52] with an exponentially decaying learning rate yields the most accurate results of the solvers tested in the least amount of training time. See
for implementation details. In single precision, the DNNs are trained for 50,000 epochs or to a tolerance of
, while in double precision DNNs are trained for 200,000 epochs or to a tolerance of
. Due to the non-monotonic convergence of minimizing the non-convex loss (2.2) with respect to the weights and biases, we checkpoint our partially trained networks once the training loss has been reduced to 1/8th of the previous checkpoint’s training loss, saving the best result at the end of training. We then average our testing results over the final trained networks.
(v) Training data and design of experiments. To understand the average performance of DL on a variety of reconstruction tasks, we run 20 trials of solving (AP) with each of our DNN architectures and CS over a range of data sets of increasing size. We then average the testing error and run statistics over all trials for each data set in plotting. We define a trial as one complete run of training a DNN, initialized as above, or solving a CS problem on a set of training data consisting of the values . To generate each data set of size
, we sample 20 i.i.d. sets of points
uniform distribution on (
and evaluate our target function f at these points to form our training data. Since we are interested in the sample complexity of the DL problem, for most of our examples we choose
to be relatively small, on the order of 1000 points.
(vi) Testing data and error metric. To study the generalization capabilities of DL, we use a common error metric for numerical analysis, the relative
where ˜f is an approximation obtained using either DL or CS. We purposefully choose this over other norms (e.g. the -norms) so that we can use the same error metric across all function classes, including both smooth and nonsmooth functions. As we run multiple trials of our experiments, we compute the average of (3.1) over all of our trials in testing. In contrast to the training data, we use deterministically generated points and function data for testing the relative
errors. More specifically, we compute an approximation to the
integrals using a high order isotropic Clenshaw Curtis sparse grid quadrature [39] rule, see, e.g., [66] for more details. As opposed to the training data, we use a large set of testing points, e.g., (d = 1) 65,537, (d = 2) 311,297, (d = 4) 643,073, and (d = 8) 1,863,937 points, to ensure a good covering of our parameter space U in computing these statistics. For such moderate dimensional problems an isotropic rule is sufficient to study the generalization performance of both CS and DL. For higher-dimensional instances of (AP), the points and weights can be pre-computed and re-used in testing multiple functions. We rely on the TASMANIAN sparse grid toolkit [77,78,79] for the generation of these rules.
(vii) Compressed sensing. We solve the problems (2.6) and (2.8) (with weights (2.9)) using the MATLAB solver SPGL1 [86,87] in double precision. The parameter is chosen as
(3.2)
To compute we use the same sparse grid rule described above in evaluating each coefficient
for further discussion. We set Λ as in (2.4) to be the hyperbolic cross index set of degree s chosen so that #(Λ
3.2. Solvers. The choice of solver and parameterization of the solver are important factors in training DNN models to a desired tolerance in a reasonable amount of time. A common choice in many computational science applications is the Adam optimizer, a variant of SGD (stochastic gradient descent) incorporating moment estimates of the gradient. In the course of testing our implementation of the MLFA package, we studied the effect of the solvers on the generalization error and time to train given a fixed budget of 50,000 epochs in single precision. Fig. 1 displays results for a variety of solvers and learning rate schedules. There we observe comparable performance for the Adam and RMSProp algorithms in terms of accuracy, with Adagrad, SGD, and PGD (proximal gradient descent) performing the worst in both accuracy and computational cost. When comparing run times and accuracy for all methods tested, the Adam optimizer achieves the best accuracy with the least computational cost. We also include results obtained with the AdamW optimizer [58], which implements a decoupled weight decay (a form of -regularization) on the weights and biases. In Fig. 1 we observe that smaller values of the weight decay parameter
optimizer to achieve identical performance to Adam with minimal overhead, but do not outperform standard Adam, while larger values of
both decrease accuracy and increase run time.
Figure 1. Comparison of (left) average relative average time in training a set of 20 ReLU
DNNs to approximate f(x) = log(sin(10x) + 2) + sin(x).
For certain solvers, it is often the case that some of the trials of a given test may fail to achieve the desired loss tolerance before arriving at the final epoch. An example of this can be seen in the left plot of Fig. 2 where none of the trials of SGD with a constant learning rate of 10were able to achieve the desired tolerance in 50,000 epochs of training. The middle plot of Fig. 2 displays the effect of using an exponentially decaying learning rate with SGD, though we also observe there that none of trials trained with this learning rate schedule were able to achieve the tolerance in 50,000 epochs of training. On the other hand the right plot shows that, on the same problem, all 20 trials trained with the Adam optimizer with the exponentially decaying learning rate were able to converge to the 5
loss tolerance in under 14,000 epochs of training. We also compared learning rate schedules for the Adam optimizer, but
found the exponentially decaying learning rate to give consistently good results. See
Figure 2. Examples of typical convergence in training a set of 20 ReLU with a constant learning rate
SGD with an exponentially decaying learning rate and (right) Adam with an exponentially decaying learning rate.
Previous studies on training DNNs suggest that the batch size can affect the convergence of the algorithms. Due to the small data set sizes considered herein over those in, e.g. computer vision, we have found that setting the batch size too small can result in longer training times (due to the increased transfer time between the CPU and GPU) and only marginal performance benefit. See . We leave a more detailed study of the affects of batch size on the performance of trained DNNs in higher-dimensional problems to a future work.
3.3. Initialization and Precision. The strategy used to initialize the weights and biases of a DNN can have a large effect on the training process, with some initializations resulting in failure to train at all. To explore the impact of this choice, we tested three different initialization strategies, all of which based on symmetric uniform or normal distributions with small variance. Empirically we find small variance to be an important factor in the success or failure of training. The strategies tested are (i) normal with mean 0 and variance 0.01, (ii) normal with mean 0 and variance 2/N, and (iii) uniform on (
Since we use constant layer widths in this work, i.e., for all hidden layers
+ 1, we can compare these strategies to the popular Xavier and He initializations, see [45], as follows. The Xavier initialization sets the weights as mean zero normal random variables with variance 1
is the number of nodes on layer
, while the He initialization modifies the variance to 2
. Therefore strategy (ii) is similar to the He initialization, with the exception of the input and output layers for which strategy (ii) prescribes much smaller variance for problems with input and output dimension in the ranges considered in this work. Also, for values of
less than 100, the Xavier initialization results in variances larger than the variance of 0.01 used in strategy (i), and similarly for values of
less than 200 for the He initialization.
Fig. 3 presents the results of these different initialization strategies for training three different architectures with Adam on the piecewise continuous function approximation problem (4.4). There we see that strategy (i) performs near the best for the smaller ReLU 2 3
30 architectures, and about as well for the 5
50 architecture as strategy (ii).
For the majority of our results, we focused on architectures with a small ratio (with L the number of layers and N the number of nodes per layer). However, we also note that setting
too small can result in DNNs which exhibit numerical instabilities after training, see Fig. 5, despite achieving the training loss tolerance. We also remark that this can occur for networks with relatively small weights, see the right plot of Fig. 14 which shows the
weights and biases for the 20 800 network remain on average bounded by 2 as we increase the samples.
We also observed that for some problems, initializing and training our DNN models in double precision can lead to improved accuracy in testing the generalization performance, and that the improvement is more pronounced for deeper networks, see Fig. 4. Due to the massive increase in computation time associated with training in double precision, a result of needing to train for more epochs to achieve the tolerance 5 and the increased time of double precision arithmetic, see item (ii) on Hardware in
, we limit our comparisons of single vs. double precision results to a handful of problems and architectures.
Figure 3. Comparison of strategies for initializing ReLU DNNs with L hidden layers and N nodes per hidden layer with (left) L = 2 and N = 20, (middle) L = 3 and N = 30, and (right) L = 5 and N = 50. Results were obtained by training with the Adam optimizer and an exponentially decaying learning rate schedule on function data generated from f(x) from (4.4) with d = 2.
Figure 4. Comparison of (a) & (b) average relative average training time for DNNs initialized and trained in single vs. double precision on a (a) & (c) one dimensional smooth function f(x) = log(sin(10x) + 2) + sin(x) and (b) & (d) eight-dimensional smooth function f(x) from (4.3).
3.4. Convergence of Adam on function approximation problems. We now discuss the convergence of the Adam optimizer. Our stopping criteria for testing the convergence depends only on the loss with respect to the training data and the current number of epochs of training. During the course of training, the process outlined above may result in networks which exhibit numerical instabilities or provide poor performance, see Fig. 5. In
that some trials may not achieve the desired accuracy tolerance within our budget of 50,000 epochs in single precision and 200,000 epochs in double precision. This can also occur for some of the trials with the Adam optimizer, see Fig. 6. Despite these issues, when viewing the average performance of the Adam optimizer on such problems we often find that the majority of networks are numerically stable and approximate the function well, see Fig. 6 which displays boxplots of the average relative error of trained ReLU DNNs with 10 hidden layers and 100 nodes per layer on a discontinuous two-dimensional function. Note that the spread shown in the boxplots is somewhat large due to the discontinuity in the underlying function. On smoother functions (not shown), this spread is smaller. We also remark that the observed performance improves with depth, see Fig. 7 which displays the outputs of all 20 trials of a variety of ReLU DNN architectures trained on a smooth, one-dimensional function.
Figure 5. Unstable ReLU networks with 20 hidden layers, trained in single precision to a loss tolerance of on noiseless data. The (left) network has 200 nodes per layer and exhibits a numerical instability (a spike) near x = 0.2334, while the (right) network has 800 nodes per layer and exhibits multiple instabilities throughout the domain, providing poor pointwise estimation of the target function.
Figure 6. Training a 10 layer DNN with 100 nodes per layer with Adam on the function (4.4) with d = 2. (left) An example where one out of 20 trials (each plotted with a different color) is unable to complete within 50,000 epochs. (middle) Boxplot of the average relative errors of 20 trials. For the boxplot, the tops and bottoms of the boxes represent the 25th and 75th percentiles, with the whiskers covering the most extreme datapoints and outliers (red plusses) plotted individually, see boxplot from MATLAB for more details.
4. Numerical experiments. We now present our numerical experiments. We test on smooth functions of one or more variables, as well as piecewise continuous functions. The results presented in this section elaborate on the main conclusions in
Figure 7. Comparison of outputs of a variety of ReLU DNN architectures trained with Adam on 750 samples of the function f(x) = log(sin(100x) + 2) + sin(10x) (plotted in black).
4.1. Smooth one-dimensional functions. We first consider problems where the target function has analytic dependence on only one variable, with regularity controllable by a single parameter. Specifically we consider
f(x) = log(sin(10Kx) + 2) + sin(Kx),(4.1)
for values of K = 1 and K = 10. When K = 1, f is smooth and has rapidly decaying Legendre coefficients, making it ideal for reconstruction with CS techniques. When K = 10, f oscillates rapidly and has slowly decaying coefficients, leading to less efficient approximation by polynomial-based methods. Fig. 8 displays results in the case of K = 1, and the right plot compares the average relative errors of the CS and DL approaches with respect to the number of samples m used in training. We observe shallower networks fail to achieve a good approximation, while increasing depth and width leads to networks which are competitive with the unweighted and weighted CS approaches. However, comparing the results obtained with the 10 and 20 hidden layer deep networks, we note diminishing returns in performance with increasing size.
Any number of factors can contribute to this phenomenon, including overfitting of the data, increased difficulty of training larger networks from our choice of random initialization, or development of numerical instabilities due to the accumulation of errors from standard sources, e.g., roundoff, overflow, and underflow. We remark that improved results have been obtained on this function by training some of our networks in double precision, see Fig. 4. However we also observe that training the 20 hidden layer DNN in double precision did not improve, but actually decreased its accuracy, suggesting room for improvement in training and initializing larger DNNs in double precision.
Fig. 9 displays the results of approximating f when K = 10, where the more rapid oscillation now leads to a degradation in performance for both methods. We again observe improved accuracy by increasing depth and width as in the previous example. We also again observe the diminishing returns increasing the size of the network architectures from 10 to 20
200, although the 20 hidden layer network provides the best accuracy of all tested DNNs after approximately 600 samples and achieves the same accuracy as the weighted CS approximations around 700 samples.
Fig. 10 displays the effect of choosing wider networks in approximating f with K = 10. There we see wider counterparts of the network architectures generally outperform narrower architectures with the same number of hidden layers. However, we also observe diminishing returns going from 10 to 20 hidden layers, and in the right plot we see the ReLU 20
Figure 8. Legendre coefficients of f(x) = log(sin(10x) + 2) + sin(x) sorted (left) lexicographically and (center) by decreasing magnitude. (right) Average relative error v.s. number of samples used in training. CS approximations were computed with the Legendre basis of cardinality n = 3, 000.
Figure 9. Legendre coefficients of f(x) = log(sin(100x) + 2) + sin(10x) sorted (left) lexicographically and (center) by decreasing magnitude. (right) Average relative error v.s. number of samples used in training. CS approximations were computed with the Legendre basis of cardinality n = 3, 000.
DNNs diverge, with the resulting trained networks exhibiting numerical artifacts as in the right plot of Fig. 5.
Fig. 11 compares the average absolute maximum of the weights and biases for ReLU DNNs trained on (4.1). We observe in the case K = 10, corresponding to the more oscillatory version of this function, on average the trained DNN architectures have larger weights and biases in magnitude when compared to those trained on the same function with K = 1. Also, when comparing the trained weights of the architectures when the ratio the function with K = 10, we note the similarity in the magnitudes of the weights despite the presence of severe artifacts for the wider networks that are not present for their narrower counterparts, see Fig. 5.
4.2. Smooth higher-dimensional functions. In Fig. 12, we present results for the function
studied in [15]. This function has smooth, isotropic dependence on its parameters and rapidly decaying Legendre coefficients with increasing polynomial order, making it ideal for approxi-
Figure 10. Average relative error v.s. number of samples of f(x) = log(sin(100x) + 2) + sin(10x) used in training DNNs with L hidden layers and N nodes per layer when
(right) 0.025. CS approximations were computed with the Legendre basis of cardinality n = 3, 000.
Figure 11. Average absolute maximum of weights and biases v.s. number of samples of f(x) =
mation with CS techniques. Despite the analyticity of f, the DNNs are unable to obtain an approximation accurate beyond 3 digits while both CS approaches achieve nearly 8 digits of accuracy. This result highlights the perceived gap between theory and practice in this work: theoretical results suggest DNNs are capable of achieving the same performance as CS on this problem, but practically we never achieve these accuracies.
Fig. 13 displays the results of approximating f from (4.2) with both narrower and wider networks. For narrower DNN architectures, corresponding to , we observe deeper networks perform better. For wider architectures corresponding to
we observe the shallower networks perform better. In all cases, the best performing networks had between 10
total trainable parameters. Also, in contrast to the diminishing returns with increasing depth observed in the one-dimensional smooth function examples, on this smooth higher-dimensional problem we now observe divergence of wider and deeper networks. These observations suggest the existence of fundamental stability barriers in current methods of training DNNs.
Fig. 14 displays the average absolute maximum of the weights and biases for each archi-
Figure 12. Legendre coefficients of f from (4.2) with d = 8 sorted (left) lexicographically and (center) by decreasing magnitude. (right) Average relative number of samples used in training. CS approximations were computed with the Legendre basis of cardinality n = 3, 023.
tecture. There we observe that the weights and biases of narrower network architectures, e.g. , remain relatively small as we increase the number of samples. However, as we increase the width of the networks to values corresponding to
begin to see the weights growing with depth as we train on larger sample sets.
Fig. 23 displays the average run time of the Adam optimizer in training our DNNs as we increase the number of training samples. For the narrower and shallower DNN architectures, we see the transfer time bottleneck between the CPU and GPU is yielding an effective linear scaling with respect to samples in the average training time. This effect is also present in the deeper and wider networks, though at a diminished amount due to the fast convergence of the Adam optimizer on larger networks.
Next we present results on the function
also from [15]. This function has anisotropic dependence on its parameters and is less smooth than the previous example. Consequently its Legendre coefficients decay more slowly, as seen in the left and middle plots of Fig. 15, and the right plot shows the CS approximations only achieve an error of roughly 0.02. The best performing ReLU DNNs achieve an error that approximately five times smaller.
Fig. 16 displays the effect of increasing the width of the DNNs as before. There, as in the previous example, we observe best performance is achieved by networks that are both wide and deep, e.g. the ReLU 5 20 networks. We also again observe that for wider networks, architectures with fewer hidden layers perform the best. Nonetheless, comparing the results between the ReLU 1
40 DNNs and those achieved by the 10
5
20 DNNs, we see that far better performance is achieved with fewer samples by the deeper and narrower DNN architectures.
Figure 17 compares the average absolute maximum of the weights and biases of the trained DNNs. There, as in the one-dimensional examples, we observe that the weights and biases
Figure 13. Comparison of average relative errors w.r.t. number of samples of f from (4.2) with d = 8 used in training ReLU architectures parameterized with
(hidden layers/nodes per hidden layer) for values
Table of best-performing architectures for each choice of
and number of parameters. CS approximations were computed with the Legendre basis of cardinality n = 3, 023.
of DNNs trained on function data from the less smooth function (4.3) are on average larger than those obtained after training on the smoother function (4.2).
Fig. 24 again displays the average run time of the Adam optimizer in training the DNNs as we increase the number of samples. There we observe similar patterns to the timing results obtained on function (4.2), e.g., linear scaling in the runtimes for narrower architectures due to the CPU-GPU transfer bottleneck. However, we also observe longer training time in general in approximating the function (4.3) compared to Fig. 23, suggesting the difficulty of approximating less smooth functions can impact training time.
4.3. Piecewise continuous functions. In this section, we present results on approximating piecewise continuous functions. We consider the function
f(x, . . . , x
, . . . , x
where and 0 otherwise. In one dimension
). We note this simple discontinuous function equally splits the data between values in {0, 1} in arbitrary dimension d, and that the separating hyperplane between the sets
is not aligned to any particular axis. In d > 1
Figure 14. Comparison of average absolute maximum of weights and biases w.r.t. number of samples of f from (4.2) with d = 8 used in training ReLU architectures parameterized with (hidden layers/nodes per hidden layer) for values
dimensions, the CS approximation fails to converge since the coefficients are not sufficiently summable. On the other hand, the work [69] shows the Heaviside function in d dimensions given by ) can be approximated to
by a 2-layer ReLU network with five nonzero weights taking value in
. As the function (4.4) is a rotation of the d-dimensional Heaviside function, the result holds in this case as well.
Fig. 18 displays the results of approximating this function in 1, 2, and 4 dimensions. There we observe the lack of convergence of the CS approximations when d > 1. While the DNNs perform better than the CS approximations for this problem, none of the achieved results obtain more than 2 digits of accuracy. In the right plot of Fig. 18, we see the double descent behavior that has been observed in other works, see, e.g., [65,
Fig. 19 displays the average absolute maximum of the weights and biases in training the DNNs. There we observe that while on average the maximum weights are larger than those found for the smooth functions of the previous section, they remain bounded for the trained DNNs. Comparing to the aforementioned existence results, we note the weights never grow large enough in our experiments to obtain such high accuracy approximations in practice. Due to the initialization strategy chosen in this work to prevent failure during training, all of our networks start from an initial point corresponding to weights and biases close to 0. Combined
Figure 15. Legendre coefficients of f from (4.3) with d = 8, sorted (left) lexicographically and (center) by decreasing magnitude. (right) Average relative error w.r.t. number of samples used in training. CS approximations were computed with the Legendre basis of cardinality n = 3, 023.
with the training process which uses an initial learning rate of 10which decays exponentially with the number of epochs, the algorithms may never reach a minimizer corresponding to higher-accuracy approximations of the halfspace function with larger weights.
5. Theoretical insights. We conclude this paper with some theoretical insights. Specifi- cally, we first establish convergence rates for polynomial approximation via CS (Theorem 5.4) and then show a practical existence theorem for DL which demonstrates the existence of an architecture and training procedure that achieves the same rates (Theorem 5.5). Proofs in this section can be found in the Supplementary Materials.
we introduced the best s-term Legendre polynomial approximation. We first establish exponential rates of convergence of such approximations. To this end, let
Bernstein ellipse with parameter
1 and define the Bernstein polyellipse
where
. This inequality is understood componentwise, i.e.
only if
. We now make the following assumption:
, the function
admits a holomorphic extension to an open set O containing
Note that this is a reasonable assumption in practice: it is known that functions that arise as quantities of interest for a range of parametric PDEs admit holomorphic extensions [17,18]. Under Assumption 5.1, it is known that the Legendre polynomial coefficients satisfy
where (see, for instance, the proof of Theorem 3.5 in [68]). When d = 1 this is a classical result in polynomial approximation, and guarantees convergence of the truncated expansion at an exponential rate
1, it also clarifies why best s-term approximation is a suitable strategy; namely, unless the parameter
is known, it is difficult to make an a priori choice of coefficients which lead to a fast rate of exponential convergence.
Figure 16. Comparison of average relative errors w.r.t. number of samples of f from (4.3) with d = 8 used in training ReLU architectures parameterized with
(hidden layers/nodes per hidden layer) for values
Table of best-performing architectures for each choice of
and number of parameters. CS approximations were computed with the Legendre basis of cardinality n = 3, 023.
On the other hand, the following theorem shows favourable exponential rates of convergence for the best s-term approximation. It also reveals another important property; namely, that the prescribed rate can be achieved using an index set that is lower. We recall that a set Λ of multi-indices is lower if, for any Λ whenever
(this inequality is understood componentwise, i.e.
). As discussed later, the lower set property is crucial for obtaining approximations using CS that attain the same rates.
For convenience, from now on we consider the s-term approximation error in the on U. This is slightly more convenient for the subsequent analysis, although
-norm bounds could also be proved with some additional technicalities.
satisfy Assumption 5.1 for some
. Then there exists a lower set Λ
of size at most s for which
Figure 17. Comparison of average absolute maximum of weights and biases w.r.t. number of samples of f from (4.3) with d = 8 used in training ReLU architectures parameterized with (hidden layers/nodes per hidden layer) for values
where depends on
only, for any
satisfying
(5.2) 0
This result asserts exponential convergence of the best s-term approximation at a rate depending on the product of log(), as opposed to log(
, which would arise if some fixed isotropic index set were used. Importantly, it also does not require any a priori knowledge of the parameters (
, and adapts to the underlying anisotropy of the function. In the next subsection we show this rate can be achieved using CS with a sampling complexity m that scales favourably with s. We refer to CS as a best-in-class scheme since there are, currently, no other methods which provably achieve this property.
Remark 5.3. There are numerous results on the convergence rate of the best s-term polynomial approximation of holomorphic functions. Algebraic rates of convergence can be found in, for instance, [16, 17], for not only Legendre, but also Chebyshev and Taylor expansions. Notably, these results also extend to the case , which theoretically permits the approximation of functions of infinitely-many variables. However, the constants in the error bounds
Figure 18. Average relative error w.r.t. number of samples used in training ReLU DNNs and solving the CS problem with Legendre basis of cardinality n in d dimensions with (left) d = 1 and n = 3, 000, (center) d = 2 and n = 3, 001, and (right) d = 4 and n = 3, 079.
Figure 19. Comparison of average absolute maximum of weights and biases vs. number of samples of the halfspace function from (4.4) used in training in (left) d = 1, (center) d = 2, and (right) d = 4 dimensions.
may be large in practice [84]. The rates shown above – which are based on [18, Sec. 3.9] and [68] – possess the advantage of always being attained in lower sets. However, they may also exhibit a poor scaling with d in the constant C. For quasi-optimal error bounds and rates, see [7,8,84]. The results shown in [84] are asymptotically sharp as , hence they provide the optimal rate in the finite-dimensional setting.
5.2. Exponential convergence via compressed sensing. Theorem 5.2 prescribes expo- nential rates of convergence for the best s-term approximation. We now show that these rates can be achieved via CS. To do so, following [2], we now suppose that the CS approximation ˆis computed by solving the so-called weighted square-root LASSO problem
minimize
rather than (2.8). The reasons for doing this are explained in
pose that 0
are drawn independently from the uniformly measure on U,
and . Then the following holds with probability at least 1
satisfy Assumption 5.1 for some
is any minimizer of (5.3) with
and weights u given by (2.9). Then
for all satisfying (5.2), where C > 0 depends on
This theorem asserts that CS achieves the same rates as those guaranteed in Theorem 5.2, with a sample complexity that is (up to the logarithmic factor) quadratic in s. In particular, scaling implied by (5.4) depends only on logarithmically on the dimension d. Hence, this estimate scales particularly well in higher dimensions.
It is important to note the key role the lower set property plays in this result. First, the union of all lower sets S of cardinality at most hyperbolic cross set of degree s. This is the rationale for choosing Λ in this way. Second, much like how sparse sets are promoted by the
-norm, sparse and lower sets are promoted by the weighted
-norm, with the weights taken to be u (these weights penalize high indices). Had one considered the unweighted
-norm in the above result, the sample complexity would have scaled with a higher algebraic power of s [3].
be as in Theorem 5.4 and suppose that
be drawn independently from the uniform measure on U. Then the following holds with probability at least 1
, satisfy Assumption 5.1 for some
. Then there exists a family of neural networks
trainable parameters and of size and depth
where is a universal constant, and a regularization functional
any minimizer ˆΦ of the regularized loss function
satisfies
for all satisfying (5.2), where C > 0 depends on
only. The functional J is a weighted
-norm penalty of the weights in the output layer.
Note that in this result the size of an architecture refers to the number of nonzero weights and biases. The number of trainable parameters in an architecture refers to the number of weights and biases that are trainable, as opposed to those that are fixed. The proof of this theorem uses a result of [68] (see Proposition C.6), which states that a finite set of Legendre polynomials can be uniformly approximated by a neural network of a given depth and size. Theorem 5.5 is obtained by rewriting the CS recovery using the weighted square-root LASSO as a neural network training problem. See for the details.
Theorem 5.5 asserts the existence of a DNN architecture, with relatively few training parameters, and a training procedure for which the resulting DNN approximations are guaranteed to perform as well as CS, up to a constant, in terms of sample complexity and convergence rates. Of course, the specific procedure suggested by the proof would not be expected to lead to any superior performance over CS. We make no claim that this approach is either practical or numerically stable (indeed, its proof relies on monomials). This is analogous to how standard existence theorems in DNN approximation theory (see ), while constructive, do not lead to superior approximations over the classical approximations on which they are based. However, it does indicate that with sufficiently careful architecture design and training, one may achieve superior performance with DNNs over CS. The extent to which this can be done is a largely open problem, requiring further theoretical and empirical investigation.
6. Conclusions. In this work we have presented results highlighting the key issues of accuracy, stability, sample complexity and computational efficiency of practical function approximation with DNNs. Our theoretical contribution on the existence of a DL procedure which performs as well as CS suggests that DL can, in theory, enjoy the same accuracy and sample complexity properties as CS. However, our numerical results comparing standard DNN architectures and common methods of training with CS techniques suggest that current methods in DL are generally unable to achieve these theoretical convergence rates on smooth function approximation problems. On the other hand, while DNNs perform relatively badly on highly smooth functions, their performance on more challenging problems, such as less smooth functions or functions with jump discontinuities is rather more promising. Certainly the fact that the same DNN architectures can be used on quite different problems sets them apart from traditional methods in scientific computing, e.g. polynomial-based methods, which are usually tied to a specific class of function (e.g. smooth functions). Hence there is ample scope and need for future work along the lines of investigation initiated in this paper. This includes both empirical investigations into architecture and cost function design, as well as algorithms for training, and further novel theoretical insights into practical existence theory; that is, the existence of not only effective DNN architectures, but also training procedure and sampling strategies which realize them efficiently. The hope is that, with these further efforts, DNNs may develop into effective tools for scientific computing that can consistently outperform current best-in-class approaches across a range of challenging problems.
7. Acknowledgements. The authors thank Clayton Webster for important discussions motivating our study of practical DNN approximation. We thank Simone Brugiapaglia for useful discussions related to various aspects of this work and Sebastian Moraga for his assistance with several technical steps in the proofs. We also acknowledge and thank the anonymous referees for proofreading and providing valuable comments on the manuscript. Finally, Boris
Hanin and Philipp Petersen are gratefully acknowledged for their helpful and inspiring comments on the approximation theory of DNNs and practical issues related to initialization and trainability.
In this section, we give further details of the testing and training setup for the DNNs and numerical experiments with CS. Section A.1 describes the hardware used in training and testing our DNN models in single and double precision. Section A.2 describes further aspects of training our DNNs with the Adam solver. Section A.3 presents results on training DNNs with Adam under a variety of learning rate schedules. Section A.4 discusses the importance of batch size on the convergence of the Adam optimizer. Section B.1 describes the selection of truncation parameters in our CS experiments.
A.1. Single versus double precision. Double precision calculations using GPUs a gen- erally more computationally demanding and hence we conducted a limited number of such studies. For example, the NVIDIA Tesla P100 GPUs operate at 4.7 TFLOPS (trillion floating point operations per second) in double precision vs. 9.3 TFLOPS in single precision, implying a 1:2 ratio for single vs. double precision computation time. On the other hand, common off-the-shelf consumer GPUs such as the NVIDIA GeForce GTX 1080 Ti operate at 0.355 TFLOPS in double precision vs. 11.5 TFLOPS in single precision, implying a 1:32 ratio for single vs. double precision computation time.
A.2. The Adam optimizer. We use the Adam optimizer [52] as follows. We use the default values , and set the initial learning rate
at epoch
we set the learning rate
where is the update frequency and
is the base, so that at the final epoch
the desired final learning rate. We allow the learning rate
to vary continuously between
, using the constants above only as scaling factors. In single precision, the DNNs are trained for
of
, while in double precision DNNs are trained for
to a tolerance of
. We set the final stepsize
, so that these choices ensure that the base b is approximately 0.85, implying the learning rate is reduced by 15% every
iterations.
A.3. Learning rate schedules for Adam. In Fig. 20 we compare learning rate schedules for the Adam optimizer on a discontinuous 2-dimensional function. The schedules compared are the exponentially decaying learning rate from (A.1), a constant learning rate of 10, and a linearly decaying learning rate, i.e., at each iteration setting
over a range of DNN architectures. We observe that shallower networks tend to perform marginally better under the exponentially-decaying learning rate schedule, while the effect is insignificant for deeper networks. In testing a range of values of the base parameter b, we found values smaller than b = 0.85 result in a learning rate that can decrease too quickly to achieve the desired tolerance with the error stagnating, while larger values of b can result in a learning rate not decreasing fast enough to quickly train larger models.
Figure 20. Comparison of learning rate schedules for the Adam optimizer for a variety of ReLU DNN architectures on function (4.4) with d = 2.
A.4. Batch size comparison. In many standard computer vision tasks, networks are trained on batches of images consisting of a small subset of the overall training set. Due to the lower dimension of the problems and smaller data set sizes considered herein, we find empirically that setting the batch size too small can result in the transfer time between the CPU and GPU contributing to a large portion of overall run time. Fig. 21 displays the effect of training four different architectures with batch sizes ranging from full-batch to quarter-batch on the generalization performance of the DNNs on a smooth 8-dimensional function. There we observe that decreasing the batch size leads to a marginal improvement in the overall error. Nonetheless for such moderate dimensional problems we find the increase in training time associated with smaller batches is not proportional to the improvement in the overall error. We leave a more detailed study of the affects of batch size on generalization performance of trained ReLU DNNs in higher-dimensional problems to a future work.
Figure 21. Comparison of batch sizes for the Adam optimizer on a smooth function approximation problem in d = 8 dimensions. For each of the plots, we fix an architecture and then train to 50,000 epochs in single precision or a tolerance of when the batch size is set to (circles) full batch (triangles) half batch (diamonds) quarter batch on function (4.3) with d = 8.
A.5. Stability of training with respect to noise. For many problems in computational science we are provided measurements that are corrupted by some amount of noise. Methods such as compressed sensing come with theoretical recovery guarantees that assert stable and robust recovery of sparse signal vectors from such measurements. Therefore, an important question for the application of DL techniques to problems in science, engineering, and medicine is how robust these methods are under realistic noise scenarios.
In Fig. 22, we include a selection of our experiments where the measurements given to both
CS and the best-performing DNNs are corrupted by additive Gaussian noise with various noise levels, i.e., given ) for various values
There we observe that both solving the CS weighted and unweighted
-minimization problems and training the DNNs are stable to noise, i.e., the
error is proportional to the noiseless approximation error plus the noise standard deviation.
Figure 22. Average relative errors of ReLU
DNNs in approximating (left) (4.1) with K = 10, and (middle) (4.2) and (right) (4.3) in d = 8 dimensions with and without
noise for various
Appendix B. Additional numerical results. In this section, we include additional observations on some of our experiments. Figure
23 displays the average run times for training ReLU DNNs of various sizes on function (4.2) with d = 8. Comparing errors from Fig. 13 with these timing results, we note that the best performing architectures in accuracy often require less time to train than some of their shallower or narrower counterparts. We also note that some of the larger architectures for which we observe divergence also required longer training times, suggesting a careful choice of architecture is key for efficient representation. Figure 24 plots the average run times for training DNNs on data from function (4.3). Comparing these results with the timings on the previous function, we observe increased training times for this less smooth function.
B.1. Truncation parameters in compressed sensing. In our numerical experiments, we solve the weighted (2.8) or unweighted (2.6) quadratically-constrained basis problems with truncation parameter chosen as in (3.2).
It is well-known that the accuracy of the quadratically-constrained basis pursuit is affected by the choice of , with it declining if
is too large or too small [12]. In (3.2), we chose
as small as possible such that the exact coefficients
are feasible for (2.6), which is in some sense an optimal choice. See [2] for further information. In practice, when the coefficients
cannot be feasibly computed, an alternative is to use cross validation, see, e.g., [30].
In our theoretical analysis in we consider the weighted square-root LASSO problem (5.3). As discussed above, the choice of
) depends on the unknown expansion coefficients
, and more precisely, the expansion tail since
. Hence, previous theoretical error estimates for polynomial approx- imations via CS often involve unrealistic assumptions on the size of this term [3, 15]. Much the same is true of the unconstrained LASSO (2.7), or its weighted variant. This problem was
Figure 23. Comparison of training time vs. number of samples of function (4.2) with d = 8 used in training for ReLU architectures parameterized with (hidden layers/nodes per hidden layer) for values (top-
studied in [12], and in [2] the weighted square-root LASSO (5.3) was proposed as a solution. While far less well known than the LASSO, the square-root LASSO has the beneficial property that the optimal choice of its parameter is independent of the noise term (i.e. the expansion tail), and depends only on the parameter s (this can be seen in Theorem 5.4). Hence, as shown in Theorem 5.4, it allows for explicit convergence rate estimates for CS, without unrealistic assumptions being imposed on the expansion tail.
On the other hand, in our experiments we continue to use (weighted) quadratically-constrained basis pursuit because it is the standard approach in the literature.
C.1. Proof of Theorem 5.2. The proof is based on techniques from [68, Sec. 3].
Proof of Theorem 5.2. It is clear from (5.1) that ). Notice the following straight- forward inequality
Figure 24. Comparison of training time vs. number of samples of function (4.3) with d = 8 used in training for ReLU architectures parameterized with (hidden layers/nodes per hidden layer) for values (top-
Hence it suffices to consider the right-hand side. Without loss of generality 0
1 such that
. In the proof of Theorem 3.5 in [68] it is shown that
for any satisfying
where C > 0 depends on only. We now derive upper and lower bounds for
in terms of s. First, observe that
and therefore
Hence
and
which gives . Therefore
Returning to (C.1), we deduce that
This completes the proof.
require some elements of compressed sensing theory, which we now introduce. We note that many of the constructions developed below apply more generally (for instance, to other measures and other orthonormal systems). However, for simplicity, we focus only on the case of Legendre polynomials. What follows is based primarily on [2,3,15,70].
Since our focus is on lower set recovery, the setup differs to the standard compressed sensing framework (see, e.g., [37]) for the recovery of arbitrary s-sparse vectors. Let recall that the union of all lower sets of cardinality at most s is the hyperbolic cross index set
, defined by (2.4). Write
. Throughout this section, we consider vectors in
indexed over Λ.
(C.2)
where
is the weighted cardinality of a subset S with respect to the weights u [70]. Note that K(s) is bounded, and satisfies
(C.4)
for weights u as in (2.9). See, for example, [2, Lem. 2.2]. Given this, we define the best s-term and lower approximation error as
Note that here and henceforth, we use to denote either the vector
equal to
and zero otherwise, or the vector
. The precise meaning will be clear from the context.
We now require the following (see [15] or [2, Defn. 5.3]):
lower robust null space property (lower rNSP) of order s if
for any is defined as in (C.2).
The lower rNSP is sufficient to provide a recovery guarantee for the weighted square-root LASSO decoder (5.3). In fact, although we shall not do it, this property also provides recovery guarantees for the decoders (2.6) and (2.7); see [2].
Theorem C.2. Suppose that satisfies the lower rNSP of order s with constants 0
and consider the the weighted square-root LASSO problem (5.3) with parameter
Then
Proof. By [2, Thm. 5.6], we have
Since is a minimizer, we obtain
and by assumption on we deduce that
as requred.
We note in passing one can also provide recovery guarantees in the 2-norm. See, for instance, [2]. In practice, it is difficult to work directly with the lower rNSP. Hence we consider the following (see [15] or [2, Defn. 5.3]):
is said to have the lower restricted isometry property of order s if there exists a constant 0
where supp(is its weighted cardinality defined as in (C.3). The smallest constant such that this property holds is called the
lower restricted isometry constant of A and it is denoted as
The following result, see [15] or [2, Lem. 5.4], asserts that the lower restricted isometry property is a sufficient condition for the lower rNSP:
and suppose
satisfies the lower restricted isometry property (with K(s) as in (C.2) for weights (2.9)) of order 2s with constant
Then A has the lower rNSP of order s with constants
Finally, we also need a result asserting the lower restricted isometry property for the measurement matrix (2.5). The following result was first proved in [15]. See, for example, [2, Thm. 5.5]:
where K(s) is as in (C.2), C > 0 is a universal constant,
and be drawn independently according to the uniform measure on U and consider the matrix
is the orthonormal tensor Legendre polynomial basis. Then, with probability at least 1
satisfies the lower RIP of order s with constant
C.3. Proof of Theorem 5.4. We now give the proof of Theorem 5.4. First, we recall the following inequality for the cardinality of the hyperbolic cross index set:
The first and third bounds are due to Theorems 3.7 and 3.5 of [14] respectively. The second bound is due to [54, Thm. 4.9]
Proof of Theorem 5.4. We claim that, with probability at least 1 , the matrix A has the lower restricted isometry property of order 2s with constant
6. Suppose this claim is true. Then Lemma C.4 gives that it satisfies the lower rNSP with constants
and
5. Also, by this and (C.4),
Recall by definition that
∥
Since this holds for all satisfying (5.2), and since the exponential term dominates as
we deduce (after possible change of C) that
To complete the proof, it remains to prove the claim. Let ) be the log factor in Theorem C.5. Then, since
Note that log(2). Using this and the estimate (C.5) for n =
Observe that for some universal constant
0 and therefore
universal constant
0. It follows that
for possibly different constant c > 0. Hence
The claim now follows immediately from Theorem C.5.
C.4. Proof of Theorem 5.5. We make use of the following result, which can be found in [68, Prop. 2.13]:
Proposition C.6. For every finite subset Λ there exists a ReLU neural network Φ
such that, if Φ
The depth and size of this network satisfy
where is a universal constant.
The general idea of the proof is to use Proposition C.6 to approximately express matrixvector multiplication Az as a neural network Φ evaluated at the sample points use the compressed sensing results to establish an error bound. Since this process commits an error, we first require the following result, which shows that the lower rNSP is robust to small matrix perturbations:
Lemma C.7. Suppose that satisfies the lower rNSP of order s with constants 0
Then satisfies the lower rNSP of order s with constants 0
Hence
The result now follows immediately.
be as in Proposition C.6. For the moment, we do not specify the choice of
1). This will be done at the end of the proof. Now define the family of neural networks
Notice that this family has size(Φ
. Then observe that
and therefore
Hence ˆis a minimizer of L over N if and only if
is a minimizer of
(C.6) minimize
We seek to use Lemma C.7. Observe that
Hence, we require to satisfy
As in the proof of Theorem 5.4, the conditions on m assert that A has the lower rNSP of order s with constants 5. Therefore, by Lemma C.7,
has the lower rNSP of order s with constants
5 with probability at least 1
Now we derive an error bound for Φ. To do so, we write ˆ
corresponding minimizer of (C.6), and set
Then following similar arguments as in proof of Theorem 5.4, we obtain
where is such that
. Moreover, we have
and now using the fact that is a minimizer of (C.6), we deduce that
Here, in the last step, we use that fact that 1 by assumption. Combining these estimates, noticing that
1 and using the same bound for
as in Theorem 5.4, we deduce that
Recall that . We now set
to deduce that
after possible change of C. It remains to establish the size and depth bounds for the network. First, from the definition of , notice that
from which, we easily obtain
Now, notice that Λ, and therefore
Proposition C.6 therefore gives
for some universal constant 1, we deduce that
which completes the proof.
[1] B. Adcock, Infinite-dimensional compressed sensing and function interpolation, Found. Comput. Math., 18 (2018), pp. 661–701.
[2] B. Adcock, A. Bao, and S. Brugiapaglia, Correcting for unknown errors in sparse high-dimensional function approximation, Numer. Math., 142 (2019), pp. 667–711.
[3] B. Adcock, S. Brugiapaglia, and C. G. Webster, Compressed sensing approaches for polynomial approximation of high-dimensional functions2017.
[4] V. Antun, F. Renna, C. Poon, B. Adcock, and A. C. Hansen, On instabilities of deep learning in image reconstruction – Does AI come at a cost?, arXiv:1902.05300, (2019).
[5] S. Arridge, P. Maass, O. ¨Oktem, and C.-B. Sch¨onlieb, Solving inverse problems using data-driven models, Acta Numer., 28 (2019), pp. 1–174.
[6] F. Bach, Breaking the Curse of Dimensionality with Convex Neural Networks, J. Mach. Learn. Res., 18 (2017), pp. 1–53.
Convergence of quasi-optimal stochastic Galerkin methods for a class of PDEs with random coefficients, Comput. Math. Appl., 67 (2014), pp. 732–751.
On the optimal polynomial approximation of stochastic PDEs by Galerkin and collocation methods, Math. Models and Methods Appl. Sci., 22 (2012).
[9] N. Baker, F. Alexander, T. Bremer, A. Hagberg, Y. Y. Kevrekidis, H. Najm, M. Parashar, A. Patra, J. Sethian, S. Wild, and K. Willcox, Workshop report on basic research needs for scientific machine learning: Core technologies for artificial intelligence, U.S. Department of Energy Advanced Scientific Computing Research, (2019).
[10] C. Beck, A. Jentzen, and B. Kuckuck, Full error analysis for the training of deep neural networks, (2019).
[11] J. Berner, P. Grohs, and A. Jentzen, Analysis of the generalization error: Empirical risk minimization over deep artificial neural networks overcomes the curse of dimensionality in the numerical approximation of Black-Scholes partial differential equations, arXiv:1809.03062, (2018).
[12] S. Brugiapaglia and B. Adcock, Robustness to unknown error in sparse regularization, IEEE Trans. Inform. Theory, 64 (2018), pp. 6638–6661.
[13] J. Carrasquilla and R. Melko, Machine learning phases of matter, Nature Physics, 13 (2017), pp. 431– 434.
New explicit-in-dimension estimates for the cardinality of high-dimensional hyperbolic crosses and approximation of functions having mixed smoothness, J. Complexity, 32 (2016), pp. 92–121.
THE GAP BETWEEN THEORY AND PRACTICE 41
[15] A. Chkifa, N. Dexter, H. Tran, and C. G. Webster, Polynomial approximation via compressed sensing of high-dimensional functions on lower sets, Math. Comp., 87 (2018), pp. 1415–1450.
[16] A. Cohen, R. DeVore, and C. Schwab, Convergence rates of best N-term Galerkin approximations for a class of elliptic sPDEs, Found. Comput. Math., 10 (2010), pp. 615–646.
[17] A. Cohen, R. DeVore, and C. Schwab, Analytic regularity and polynomial approximation of parametric and stochastic elliptic PDEs, Anal. Appl., 9 (2011), pp. 11–47.
[18] A. Cohen and R. A. DeVore, Approximation of high-dimensional parametric PDEs, Acta Numer., 24 (2015), pp. 1–159.
[19] G. Cybenko, Approximation by Superpositions of a Sigmoidal Function, Math. Control Signals Systems, 2 (1989), pp. 303–314.
[20] E. C. Cyr, M. A. Gulian, R. G. Patel, M. Perego, and N. A. Trask, Robust training and initialization of deep neural networks: An adaptive basis viewpoint, vol. 107 of Proceedings of Machine Learning Research, Princeton University, Princeton, NJ, USA, 20–24 Jul 2020, PMLR, pp. 512–536, http://proceedings.mlr.press/v107/cyr20a.html.
[21] G. E. Dahl, D. Yu, L. Deng, and A. Acero, Context-Dependent Pre-Trained Deep Neural Networks for Large-Vocabulary Speech Recognition, IEEE Trans. Audio, Speech, Language Process., 20 (2012), pp. 30–42.
[22] J. Daws and C. Webster, Analysis of Deep Neural Networks with Quasi-optimal polynomial approximation rates, 107 (2019), pp. 1–13.
[23] J. Daws and C. G. Webster, A Polynomial-Based Approach for Architectural Design and Learning with Deep Neural Networks, (2019).
[24] J. De Fauw, J. R. Ledsam, B. Romera-Paredes, S. Nikolov, N. Tomasev, S. Blackwell, H. Askham, X. Glorot, B. O’Donoghue, D. Visentin, G. van den Driessche, B. Lakshminarayanan, C. Meyer, F. Mackinder, S. Bouton, K. Ayoub, R. Chopra, D. King, A. Karthikesalingam, C. O. Hughes, R. Raine, J. Hughes, D. A. Sim, C. Egan, A. Tufail, H. Montgomery, D. Hassabis, G. Rees, T. Back, P. T. Khaw, M. Suleyman, J. Cornebise, P. A. Keane, and O. Ronneberger, Clinically applicable deep learning for diagnosis and referral in retinal disease, Nature Medicine, 24 (2018), pp. 1342–1350.
[25] J. Deng, W. Dong, R. Socher, L.-J. Li, K. Li, and L. Fei-Fei, ImageNet: A Large-Scale Hierarchical Image Database, in 2009 IEEE Conference on Computer Vision and Pattern Recognition, 2009, pp. 248–255.
[26] A. Dereventsov, A. Petrosyan, and C. Webster, Greedy Shallow Networks: A New Approach for Constructing and Training Neural Networks, (2019).
[27] A. Dereventsov, A. Petrosyan, and C. Webster, Neural network integral representations with the ReLU activation function, (2019), pp. 1–15.
[28] R. A. DeVore, Nonlinear approximation, Acta Numer., 7 (1998), pp. 51–150.
regularization approach for sparse simultaneous approximation of parameterized PDEs, ESAIM Math. Model. Numer. Anal., 53 (2019), pp. 2025–2045.
[30] A. Doostan and H. Owhadi, A non-adapted sparse approximation of PDEs with stochastic inputs, J. Comput. Phys., 230 (2011), pp. 3015–3034.
[31] W. E, J. Han, and A. Jentzen, Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations, Commun. Math. Stat., 5 (2017), pp. 349–380.
[32] W. E, C. Ma, and L. Wu, Barron spaces and the compositional function spaces for neural network models, 2019, https://arxiv.org/abs/1906.08039.
[33] W. E and Q. Wang, Exponential convergence of the deep neural network approximation for analytic functions, 2018, https://arxiv.org/abs/1807.00297.
[34] C. Farabet, C. Couprie, L. Najman, and Y. LeCun, Scene parsing with multiscale feature learning, purity trees, and optimal covers, (2012).
[35] A. Fawzi, S.-M. Moosavi-Dezfooli, and P. Frossard, The robustness of deep networks: A geometrical perspective, IEEE Signal Process. Mag., 34 (2017), pp. 50–62.
[36] D. Fokina and I. Oseledets, Growing axons: greedy learning of neural networks with application to function approximation, (2019), pp. 1–17.
[37] S. Foucart and H. Rauhut, A Mathematical Introduction to Compressive Sensing, Applied and Nu-
42 B. ADCOCK AND N. DEXTER
merical Harmonic Analysis, Birkh¨auser, 2013.
[38] M. Geist, P. Petersen, M. Raslan, R. Schneider, and G. Kutyniok, Numerical solution of the parametric diffusion equation by deep neural networks, 2020, https://arxiv.org/abs/2004.12131.
[39] T. Gerstner and M. Griebel, Numerical integration using sparse grids, Numer. Algorithms, 18 (1998), pp. 209–232.
[40] X. Glorot and Y. Bengio, Understanding the difficulty of training deep feedforward neural networks, Journal of Machine Learning Research, 9 (2010), pp. 249–256.
[41] I. Goodfellow, Y. Bengio, and A. Courville, Deep Learning, The MIT Press, 2016.
[42] N. M. Gottschling, V. Antun, B. Adcock, and A. C. Hansen, The troublesome kernel: why deep learning for inverse problems is typically unstable, arXiv:2001.01258, (2020).
[43] P. Grohs, F. Hornung, A. Jentzen, and P. Von Wurstemberger, A proof that artificial neural networks overcome the curse of dimensionality in the numerical approximation of Black-Scholes partial differential equations, arXiv:1809.02362, (2018).
[44] P. Grohs, T. Wiatowski, and H. Bolcskei, Deep convolutional neural networks on cartoon functions, in IEEE International Symposium on Information Theory, ISIT 2016, Barcelona, Spain, July 10-15, 2016, 2016, pp. 1163–1167.
Error bounds for approximations with deep ReLU neural networks in
[46] M. D. Gunzburger, C. G. Webster, and G. Zhang, Stochastic finite element methods for partial differential equations with random input data, Acta Numer., 23 (2014), pp. 521–650.
[47] B. Hanin, Which neural net architectures give rise to exploding and vanishing gradients?, Advances in Neural Information Processing Systems, (2018), pp. 582–591.
[48] B. Hanin and D. Rolnick, How to start training: The effect of initialization and architecture, in Advances in Neural Information Processing Systems, 2018, pp. 571–581.
[49] K. He, X. Zhang, S. Ren, and J. Sun, Delving deep into rectifiers: Surpassing human-level performance on imagenet classification, Proceedings of the IEEE International Conference on Computer Vision, 2015 International Conference on Computer Vision, ICCV 2015 (2015), pp. 1026–1034.
[50] G. Hinton, L. Deng, D. Yu, G. E. Dahl, A. Mohamed, N. Jaitly, A. Senior, V. Vanhoucke, P. Nguyen, T. N. Sainath, and B. Kingsbury, Deep neural networks for acoustic modeling in speech recognition: The shared views of four research groups, IEEE Signal Process. Mag., 29 (2012), pp. 82–97.
[51] K. Hornik, M. Stinchcombe, and H. White, Multilayer feedforward networks are universal approximators (Book review), Neural Networks, 2 (1989), pp. 359–366.
[52] D. P. Kingma and J. Ba, Adam: A method for stochastic optimization, arXiv:1412.6980, (2014).
[53] A. Krizhevsky, I. Sutskever, and G. E. Hinton, ImageNet Classification with Deep Convolutional Neural Networks, in Advances in Neural Information Processing Systems, 2012, pp. 1097–1105.
Approximation of mixed order Sobolev functions on the d-torus: asymptotics, preasymptotics, and d-dependence, Constr. Approx., 42 (2015), pp. 353–398.
[55] G. Kutyniok, P. Petersen, M. Raslan, and R. Schneider, A theoretical analysis of deep neural networks and parametric pdes, 2019, https://arxiv.org/abs/1904.00377.
[56] B. Li, S. Tang, and H. Yu, Better approximations of high dimensional smooth functions by deep neural networks with rectified power units, 1903.05858, (2019).
[57] S. Liang and R. Srikant, Why deep neural networks for function approximation?, arXiv:1610.04161, (2016).
[58] I. Loshchilov and F. Hutter, Decoupled weight decay regularization, in 7th International Conference on Learning Representations, ICLR 2019, 2019.
[59] J.-L. Loyer, E. Henriques, M. Fontul, and S. Wiseall, Comparison of machine learning methods applied to the estimation of manufacturing cost of jet engine components, Int. J. Prod. Econ., 178 (2016), pp. 109–119.
[60] J. Lu, Z. Shen, H. Yang, and S. Zhang, Deep Network Approximation for Smooth Functions, (2020), pp. 1–33.
[61] Y. Lu, A. Zhong, Q. Li, and B. Dong, Beyond finite layer neural networks: bridging deep architectures and numerical differential equations, arXiv:1710.10121, (2017).
[62] H. Montanelli, H. Yang, and Q. Du, Deep ReLU networks overcome the curse of dimensionality for
THE GAP BETWEEN THEORY AND PRACTICE 43
bandlimited functions, 1903.00735, (2019).
[63] S.-M. Moosavi-Dezfooli, A. Fawzi, O. Fawzi, and P. Frossard, Universal adversarial perturbations, (2016).
[64] S.-M. Moosavi-Dezfooli, A. Fawzi, and P. Frossard, Deepfool: a simple and accurate method to fool deep neural networks, (2015).
[65] P. Nakkiran, G. Kaplun, Y. Bansal, T. Yang, B. Barak, and I. Sutskever, Deep Double Descent: Where Bigger Models and More Data Hurt, arXiv:1912.02292, (2019).
[66] F. Nobile, R. Tempone, and C. Webster, A sparse grid stochastic collocation method for elliptic partial differential equations with random input data, SIAM J. Numer. Anal., 46 (2008), pp. 2309– 2345.
[67] J. A. A. Opschoor, P. C. Petersen, and C. Schwab, Deep ReLU networks and high-order finite element methods
[68] J. A. A. Opschoor, C. Schwab, and J. Zech, Exponential ReLU DNN expression of holomorphic maps in high dimension
[69] P. Petersen and F. Voigtlaender, Optimal approximation of piecewise smooth functions using deep ReLU neural networks, Neural Netw., 108 (2018), pp. 296–330.
[70] , Appl. Comput. Harm. Anal., 40 (2016), pp. 321–351.
[71] S. H. Rudy, S. L. Brunton, J. L. Proctor, and J. N. Kutz, Data-driven discovery of partial differential equations, Sci. Adv., 3 (2017).
[72] C. Schwab and J. Zech, Deep learning in high dimension: Neural network expression rates for generalized polynomial chaos expansions in UQ, Anal. Appl., 17 (2019), pp. 19–55.
[73] Z. Shen, H. Yang, and S. Zhang, Deep network approximation characterized by number of neurons, arXiv:1906.05497, (2019).
[74] D. Silver, J. Schrittwieser, K. Simonyan, I. Antonoglou, A. Huang, A. Guez, T. Hubert, L. Baker, M. Lai, A. Bolton, Y. Chen, T. Lillicrap, F. Hui, L. Sifre, G. Van Den Driessche, T. Graepel, and D. Hassabis, Mastering the game of Go without human knowledge, Nature, 550 (2017), pp. 354–359.
[75] K. Simonyan and A. Zisserman, Very deep convolutional networks for large-scale image recognition, in International Conference on Learning Representations (ICLR), 2015.
[76] C. Sommer and D. W. Gerlich, Machine learning in cell biology - teaching computers to recognize phenotypes, Journal of cell science, 126 (2013).
[77] M. Stoyanov, User manual: Tasmanian sparse grids, Tech. Report ORNL/TM-2015/596, Oak Ridge National Laboratory, One Bethel Valley Road, Oak Ridge, TN, 2015.
[78] M. Stoyanov, Adaptive sparse grid construction in a context of local anisotropy and multiple hierarchical parents, in Sparse Grids and Applications-Miami 2016, Springer, 2018, pp. 175–199.
[79] M. K. Stoyanov and C. G. Webster, A dynamically adaptive sparse grids method for quasi-optimal interpolation of multidimensional functions, Computers Math. Applic., 71 (2016), pp. 2449–2465.
[80] E. Strubell, A. Ganesh, and A. McCallum, Energy and policy considerations for deep learning in NLP, arXiv:1906.02243, (2019).
[81] C. Szegedy, W. Zaremba, I. Sutskever, J. Bruna, D. Erhan, I. Goodfellow, and R. Fergus, Intriguing properties of neural networks, (2013).
[82] W. Z. Taffese and E. Sistonen, Machine learning for durability and service-life assessment of reinforced concrete structures: Recent advances and future directions, Automation in Construction, 77 (2017), pp. 1 – 14.
[83] A. L. Tarca, V. J. Carey, X.-w. Chen, R. Romero, and S. Dr˘aghici, Machine learning and its applications to biology, PLoS Computational Biology, 3 (2007).
[84] H. Tran, C. G. Webster, and G. Zhang, Analysis of quasi-optimal polynomial approximations for parameterized PDEs with deterministic and stochastic coefficients, Numer. Math., (2017), pp. 1–43.
[85] M. Unser, A representer theorem for deep neural networks, J. Mach. Learn. Res., 20 (2019), pp. 1–28.
[86] E. van den Berg and M. P. Friedlander, SPGL1: A solver for large-scale sparse reconstruction, June 2007. http://www.cs.ubc.ca/labs/scl/spgl1.
[87] E. van den Berg and M. P. Friedlander, Probing the pareto frontier for basis pursuit solutions, SIAM J. Sci. Comput., 31 (2008), pp. 890–912.
44 B. ADCOCK AND N. DEXTER
[88] E. Weinan and B. Yu, The Deep Ritz Method: A Deep Learning-Based Numerical Algorithm for Solving Variational Problems, Commun. Math. Stat., 6 (2018), pp. 1–14.
[89] C. Wu, P. Karanasou, M. J. Gales, and K. C. Sim, Stimulated deep neural network for speech recognition, in Interspeech 2016, 2016, pp. 400–404.
[90] D. Yarotsky, Error bounds for approximations with deep ReLU networks, Neural Netw., 94 (2017), pp. 103–114.
[91] D. Yarotsky, Optimal approximation of continuous functions by very deep ReLU networks, arXiv:1802.03620, (2018).
[92] G. Zhang, J. Zhang, and J. Hinkle, Learning nonlinear level sets for dimensionality reduction in function approximation, in Advances in Neural Information Processing Systems 32, Curran Associates, Inc., 2019, pp. 13199–13208.
[93] B. Zieli´nski, A. Plichta, K. Misztal, P. Spurek, M. Brzychczy-W�loch, and D. Ocho´nska, Deep learning approach to bacterial colony classification, PLOS ONE, 12 (2017), p. e0184554.