Where can I find a paper  explaining the motivation behind using exponential moving averaged weights in a neural network?
https://arxiv.org/abs/1803.05407
Averaging Weights Leads to Wider Optima and Better Generalization
2018-03-14
Izmailov; Pavel; Podoprikhin; Dmitrii; Garipov; Timur; Vetrov; Dmitry; Wilson; Andrew Gordon
<div><div>
 
 
 
 
 
 
 <p><a href="/pdf/1803.05407">Download PDF</a></p><blockquote>
 Abstract: Deep neural networks are typically trained by optimizing a loss function with
an SGD variant, in conjunction with a decaying learning rate, until
convergence. We show that simple averaging of multiple points along the
trajectory of SGD, with a cyclical or constant learning rate, leads to better
generalization than conventional training. We also show that this Stochastic
Weight Averaging (SWA) procedure finds much flatter solutions than SGD, and
approximates the recent Fast Geometric Ensembling (FGE) approach with a single
model. Using SWA we achieve notable improvement in test accuracy over
conventional SGD training on a range of state-of-the-art residual networks,
PyramidNets, DenseNets, and Shake-Shake networks on CIFAR-10, CIFAR-100, and
ImageNet. In short, SWA is extremely easy to implement, improves
generalization, and has almost no computational overhead.
 </blockquote>
 
 
 </div><div>
 <h2>Submission history</h2><p> From: Andrew Wilson [<a href="/show-email/e92e5f48/1803.05407">view email</a>]
 <br />
 <strong><a href="/abs/1803.05407v1">[v1]</a></strong>
 Wed, 14 Mar 2018 17:09:27 UTC (1,402 KB)<br />
 <strong><a href="/abs/1803.05407v2">[v2]</a></strong>
 Wed, 8 Aug 2018 08:49:15 UTC (1,404 KB)<br /><strong>[v3]</strong>
Mon, 25 Feb 2019 14:18:11 UTC (1,404 KB)<br /></p></div>||||I|||| Skip to main content
 We gratefully acknowledge support from
 the Simons Foundation and member institutions.
 > cs > arXiv:1803.05407

 Help | Advanced Search

 All fields Title Author Abstract Comments Journal reference ACM classification MSC classification Report number arXiv identifier DOI ORCID arXiv author ID Help pages Full text
 Search
 GO

 quick links

 * Login
 * Help Pages
 * About

 Computer Science > Machine Learning

 arXiv:1803.05407 (cs)
 [Submitted on 14 Mar 2018 (v1), last revised 25 Feb 2019 (this version, v3)]

 Title: Averaging Weights Leads to Wider Optima and Better Generalization

 Authors: Pavel Izmailov, Dmitrii Podoprikhin, Timur Garipov, Dmitry Vetrov, Andrew Gordon Wilson
 Download PDF
 Abstract: Deep neural networks are typically trained by optimizing a loss function with an SGD variant, in conjunction with a decaying learning rate, until convergence. We show that simple averaging of multiple points along the trajectory of SGD, with a cyclical or constant learning rate, leads to better generalization than conventional training. We also show that this Stochastic Weight Averaging (SWA) procedure finds much flatter solutions than SGD, and approximates the recent Fast Geometric Ensembling (FGE) approach with a single model. Using SWA we achieve notable improvement in test accuracy over conventional SGD training on a range of state-of-the-art residual networks, PyramidNets, DenseNets, and Shake-Shake networks on CIFAR-10, CIFAR-100, and ImageNet. In short, SWA is extremely easy to implement, improves generalization, and has almost no computational overhead.
 Comments: Appears at the Conference on Uncertainty in Artificial Intelligence (UAI), 2018 
 Subjects: Machine Learning (cs.LG) ; Artificial Intelligence (cs.AI); Computer Vision and Pattern Recognition (cs.CV); Machine Learning (stat.ML)
 Cite as: arXiv:1803.05407 [cs.LG]  
 (or arXiv:1803.05407v3 [cs.LG] for this version)  
 https://doi.org/10.48550/arXiv.1803.05407  
 Focus to learn more  
 arXiv-issued DOI via DataCite  
 

 Submission history

 From: Andrew Wilson [view email]
 [v1] Wed, 14 Mar 2018 17:09:27 UTC (1,402 KB)
 [v2] Wed, 8 Aug 2018 08:49:15 UTC (1,404 KB)
 [v3] Mon, 25 Feb 2019 14:18:11 UTC (1,404 KB)
 Full-text links:

 Download:

 * PDF
 * Other formats
 (license)
 Current browse context:
 cs.LG
 < prev | next >
 new | recent | 1803
 Change to browse by:
 cs
 cs.AI
 cs.CV
 stat
 stat.ML

 References & Citations

 * NASA ADS
 * Google Scholar
 * Semantic Scholar

 5 blog links

 (what is this?)

 DBLP - CS Bibliography

 listing | bibtex
 Pavel Izmailov
 Dmitrii Podoprikhin
 Timur Garipov
 Dmitry P. Vetrov
 Andrew Gordon Wilson
 a export bibtex citation Loading...

 Bibtex formatted citation

 ×
 loading...
 Data provided by:

 Bookmark

 Bibliographic Tools

 Bibliographic and Citation Tools

 Bibliographic Explorer Toggle
 Bibliographic Explorer (What is the Explorer?)
 Litmaps Toggle
 Litmaps (What is Litmaps?)
 scite.ai Toggle
 scite Smart Citations (What are Smart Citations?)
 Code, Data, Media

 Code, Data and Media Associated with this Article

 Links to Code Toggle
 Papers with Code (What is Papers with Code?)
 ScienceCast Toggle
 ScienceCast (What is ScienceCast?)
 Demos

 Demos

 Replicate Toggle
 Replicate (What is Replicate?)
 Spaces Toggle
 Hugging Face Spaces (What is Spaces?)
 Related Papers

 Recommenders and Search Tools

 Connected Papers Toggle
 Connected Papers (What is Connected Papers?)
 Core recommender toggle
 CORE Recommender (What is CORE?)
 IArxiv recommender toggle
 IArxiv Recommender (What is IArxiv?)
 About arXivLabs

 arXivLabs: experimental projects with community collaborators

 arXivLabs is a framework that allows collaborators to develop and share new arXiv features directly on our website.

 Both individuals and organizations that work with arXivLabs have embraced and accepted our values of openness, community, excellence, and user data privacy. arXiv is committed to these values and only works with partners that adhere to them.

 Have an idea for a project that will add value for arXiv's community? Learn more about arXivLabs.

 Which authors of this paper are endorsers? | Disable MathJax (What is MathJax?)
 * About
 * Help
 * Click here to contact arXiv Contact
 * Click here to subscribe Subscribe
 * Copyright
 * Privacy Policy
 * Web Accessibility Assistance

 * arXiv Operational Status
 Get status notifications via email or slack

Averaging Weights Leads to Wider Optima and Better Generalization
Pavel Izmailov 
Cornell University
Dmitrii Podoprikhin 
Higher School of Economics
Samsung-HSE Laboratory
Timur Garipov 
Samsung AI Center in Moscow
Lomonosov Moscow State University
Dmitry Vetrov 
Higher School of Economics
Samsung-HSE Laboratory
Andrew Gordon Wilson 
Cornell University
Averaging Weights Leads to Wider Optima and Better Generalization
Deep neural networks are typically trained by optimizing a loss function with an SGD variant, in conjunction with a decaying learning rate, until convergence. We show that simple averaging of multiple points along the trajectory of SGD, with a cyclical or constant learning rate, leads to better generalization than conventional training. We also show that this Stochastic Weight Averaging (SWA) procedure finds much broader optima than SGD, and approximates the recent Fast Geometric Ensembling (FGE) approach with a single model. Using SWA we achieve notable improvement in test accuracy over conventional SGD training on a range of state-of-the-art residual networks, PyramidNets, DenseNets, and Shake-Shake networks on CIFAR-10, CIFAR-100, and ImageNet. In short, SWA is extremely easy to implement, improves generalization, and has almost no computational overhead.
INTRODUCTION
With a better understanding of the loss surfaces for multilayer networks, we can accelerate the convergence, stability, and accuracy of training procedures in deep learning. Recent work [Garipov et al., 2018, Draxler et al., 2018 shows that local optima found by SGD can be connected by simple curves of near constant loss. Building upon this insight, Garipov et al. [2018] also developed Fast Geometric Ensembling (FGE) to sample multiple nearby points in weight space to create high performing ensembles in the time required to train a single DNN.
FGE uses a high frequency cyclical learning rate with SGD to select networks to ensemble. In Figure 1 (left) * Equal contribution. we see that the weights of the networks ensembled by FGE are on the periphery of the most desirable solutions. This observation suggests it is promising to average these points in weight space, and use a network with these averaged weights, instead of forming an ensemble by averaging the outputs of networks in model space. Although the general idea of maintaining a running average of weights traversed by SGD dates back to Ruppert [1988], this procedure is not typically used to train neural networks. It is sometimes applied as an exponentially decaying running average in combination with a decaying learning rate (where it is called an exponential moving average), which smooths the trajectory of conventional SGD but does not perform very differently. However, we show that an equally weighted average of the points traversed by SGD with a cyclical or constant learning rate, which we refer to as Stochastic Weight Averaging (SWA), has many surprising and promising features for training deep neural networks, leading to a better understanding of the geometry of their loss surfaces. Indeed, SWA with cyclical or constant learning rates can be used as a dropin replacement for standard SGD training of multilayer networks -but with improved generalization and essentially no overhead. In particular:
• We show that SGD with cyclical [e.g., Loshchilov and Hutter, 2017] and constant learning rates traverses regions of weight space corresponding to high-performing networks. We find that while these models are moving around this optimal set they never reach its central points. We show that we can move into this more desirable space of points by averaging the weights proposed over SGD iterations.
• While FGE ensembles [Garipov et al., 2018] can be trained in the same time as a single model, test predictions for an ensemble of k models requires k times more computation. We show that SWA can be interpreted as an approximation to FGE ensembles but with the test-time, convenience, and inter- pretability of a single model.
• We demonstrate that SWA leads to solutions corresponding to wider optima than SGD. Keskar et al.
[2017] and Hochreiter and Schmidhuber [1997] conjecture that the width of the optima is critically related to generalization. We illustrate that the loss on the train is shifted with respect to the test error ( Figure 1, middle and right panels, and sections 3, 4). We show that SGD generally converges to a point near the boundary of the wide flat region of optimal points. SWA on the other hand is able to find a point centered in this region, often with slightly worse train loss but with substantially better test error.
• SWA achieves notable improvement for training a broad range of architectures over several consequential benchmarks. In particular, running SWA for just 10 epochs on ImageNet we are able to achieve 0.8% improvement for ResNet-50 and DenseNet-161, and 0.6% improvement for ResNet-150. We achieve improvement of over 1.3% on CIFAR-100 and of over 0.4% on CIFAR-10 with Preactivation ResNet-164, VGG-16 and Wide ResNet-28-10. We also achieve substantial improvement for the recent Shake-Shake Networks and PyramidNets.
• SWA is extremely easy to implement and has virtually no computational overhead compared to the conventional training schemes.
• We provide an implementation of SWA at https://github.com/timgaripov/swa.
RELATED WORK
This paper is fundamentally about better understanding the geometry of loss surfaces and generalization in deep learning. We follow the trajectory of weights traversed by SGD, leading to new geometric insights and the intuition that SWA will lead to better results than standard training. Empirically, we make the discovery that SWA notably improves training of many state-of-the-art deep neural networks over a range of consequential benchmarks, with essentially no overhead.
The procedures for training neural networks are constantly being improved. New methods are being proposed for architecture design, regularization and optimization. The SWA approach is related to work in both optimization and regularization.
In optimization, there is great interest in how different types of local optima affect generalization in deep learning. Keskar et al. [2017] claim that SGD is more likely to converge to broad local optima than batch gradient methods, which tend to converge to sharp optima. Moreover, they argue that the broad optima found by SGD are more likely to have good test performance, even if the training loss is worse than for the sharp optima. On the other hand Dinh et al. [2017] argue that all the known definitions of sharpness are unsatisfactory and cannot on their own explain generalization. Chaudhari et al. [2017] propose the Entropy-SGD method that explicitly forces optimization towards wide valleys. They report that although the optima found by Entropy-SGD are wider than those found by conventional SGD, the generalization performance is still comparable.
The SWA method is based on averaging multiple points along the trajectory of SGD with cyclical or constant learning rates. The general idea of maintaining a running average of weights proposed by SGD was first consid-1 Suppose we have three weight vectors w1, w2, w3. We set
u = (w2 −w1), v = (w3 −w1)− w3 − w1, w2 − w1 / w2 − w1 2 · (w2 − w1)
. Then the normalized vectorsû = u/ u , v = v/ v form an orthonormal basis in the plane containing w1, w2, w3. To visualize the loss in this plane, we define a Cartesian grid in the basisû,v and evaluate the networks corresponding to each of the points in the grid. A point P with coordinates (x, y) in the plane would then be given by P = w1 + x ·û + y ·v. ered in convex optimization by Ruppert [1988] and later by Polyak and Juditsky [1992]. However, this procedure is not typically used to train neural networks. Practitioners instead sometimes use an exponentially decaying running average of the weights found by SGD with a decaying learning rate, which smooths the trajectory of SGD but performs comparably.
SWA is making use of multiple samples gathered through exploration of the set of points corresponding to high performing networks. To enforce exploration we run SGD with constant or cyclical learning rates. Mandt et al. [2017] show that under several simplifying assumptions running SGD with a constant learning rate is equivalent to sampling from a Gaussian distribution centered at the minimum of the loss, and the covariance of this Gaussian is controlled by the learning rate. Following this explanation from [Mandt et al., 2017], we can interpret points proposed by SGD as being constrained to the surface of a sphere, since they come from a high dimensional Gaussian distribution. SWA effectively allows us to go inside the sphere to find higher density solutions.
In a procedure called Fast Geometric Ensembling (FGE), Garipov et al. [2018] showed that using a cyclical learning rate it is possible to gather models that are spatially close to each other but produce diverse predictions. They used the gathered models to train ensembles with no computational overhead compared to training a single DNN model. In recent work Neklyudov et al. [2018] also discuss an efficient approach for model averaging of Bayesian neural networks. SWA was inspired by following the trajectories of FGE proposals, in order to find a single model that would approximate an FGE ensemble, but provide greater interpretability, convenience, and test-time scalability.
Dropout [Srivastava et al., 2014] is an extremely popular approach to regularizing DNNs. Across each minibatch used for SGD, a different architecture is created by randomly dropping out neurons. The authors make analogies between dropout, ensembling, and Bayesian model averaging. At test time, an ensemble approach is proposed, but then approximated with similar results by multiplying each connection by the dropout rate. At a high level, SWA and Dropout are both at once regularizers and training procedures, motivated to approximate an ensemble. Each approach implements these high level ideas quite differently, and as we show in our experiments, can be combined for improved performance.
STOCHASTIC WEIGHT AVERAGING
We present Stochastic Weight Averaging (SWA) and analyze its properties. In section 3.1, we consider trajec-tories of SGD with a constant and cyclical learning rate, which helps understand the geometry of SGD training for neural networks, and motivates the SWA procedure. Then in section 3.2 we present the SWA algorithm in detail, in section 3.3 we derive its complexity, and in section 3.4 we analyze the width of optima found by SWA versus conventional SGD training. In section 3.5 we then examine the relationship between SWA and the recently proposed Fast Geometric Ensembling [Garipov et al., 2018]. Finally, in section 3.6 we consider SWA from the perspective of stochastic convex optimization.
We note the name SWA has two meanings: on the one hand, it is an average of SGD weights. On the other, with a cyclical or constant learning rate, SGD proposals are approximately sampling from the loss surface of the DNN, leading to stochastic weights.
ANALYSIS OF SGD TRAJECTORIES
SWA is based on averaging the samples proposed by SGD using a learning rate schedule that allows exploration of the region of weight space corresponding to high-performing networks. In particular we consider cyclical and constant learning rate schedules.
The cyclical learning rate schedule that we adopt is inspired by Garipov et al. [2018] and Smith and Topin [2017]. In each cycle we linearly decrease the learning rate from α 1 to α 2 . The formula for the learning rate at iteration i is given by
α(i) = (1 − t(i))α 1 + t(i)α 2 , t(i) = 1 c (mod(i − 1, c) + 1) .
The base learning rates α 1 ≥ α 2 and the cycle length c are the hyper-parameters of the method. Here by iteration we assume the processing of one batch of data. Figure 2 illustrates the cyclical learning rate schedule and the test error of the corresponding points. Note that unlike the cyclical learning rate schedule of Garipov et al. [2018] and Smith and Topin [2017], here we propose to use a discontinuous schedule that jumps directly from the minimum to maximum learning rates, and does not steadily increase the learning rate as part of the cycle. We use this more abrupt cycle because for our purposes exploration is more important than the accuracy of individual proposals. For even greater exploration, we also consider constant learning rates α(i) = α 1 .
We run SGD with cyclical and constant learning rate schedules starting from a pretrained point for a Preactivation ResNet-164 on CIFAR-100. We then use the first, middle and last point of each of the trajectories to define a 2-dimensional plane in the weight space containing all affine combinations of these points. In Figure 3  we plot the loss on train and error on test for points in these planes. We then project the other points of the trajectory to the plane of the plot. Note that the trajectories do not generally lie in the plane of the plot, except for the first, last and middle points, showed by black crosses in the figure. Therefore for other points of the trajectories it is not possible to tell the value of train loss and test error from the plots.
The key insight from Figure 3 is that both methods explore points close to the periphery of the set of highperforming networks. The visualizations suggest that both methods are doing exploration in the region of space corresponding to DNNs with high accuracy. The main difference between the two approaches is that the individual proposals of SGD with a cyclical learning rate schedule are in general much more accurate than the proposals of a fixed-learning rate SGD. After making a large step, SGD with a cyclical learning rate spends several epochs fine-tuning the resulting point with a decreasing learning rate. SGD with a fixed learning rate on the other hand is always making steps of relatively large sizes, exploring more efficiently than with a cyclical learning rate, but the individual proposals are worse.
Another important insight we can get from Figure 3 is that while the train loss and test error surfaces are qualitatively similar, they are not perfectly aligned. The shift between train and test suggests that more robust central points in the set of high-performing networks can lead to better generalization. Indeed, if we average several proposals from the optimization trajectories, we get a more robust point that has a substantially higher test performance than the individual proposals of SGD, and is es-sentially centered on the shifted mode for test error. We further discuss the reasons for this behaviour in sections 3.4, 3.5, 3.6.
SWA ALGORITHM
We now present the details of the Stochastic Weight Averaging algorithm, a simple but effective modification for training neural networks, motivated by our observations in section 3.1.
Following Garipov et al. [2018], we start with a pretrained modelŵ. We will refer to the number of epochs required to train a given DNN with the conventional training procedure as its training budget and will denote it by B. The pretrained modelŵ can be trained with the conventional training procedure for full training budget or reduced number of epochs (e.g. 0.75B). In the latter case we just stop the training early without modifying the learning rate schedule. Starting fromŵ we continue training, using a cyclical or constant learning rate schedule. When using a cyclical learning rate we capture the models w i that correspond to the minimum values of the learning rate (see Figure 2), following Garipov et al. [2018]. For constant learning rates we capture models at each epoch. Next, we average the weights of all the captured networks w i to get our final model w SWA .
Note that for cyclical learning rate schedule, the SWA algorithm is related to FGE [Garipov et al., 2018], except that instead of averaging the predictions of the models, we average their weights, and we use a different type of learning rate cycle. In section 3.5 we show how SWA can approximate FGE, but with a single model.
Batch normalization. If the DNN uses batch normalization [Ioffe and Szegedy, 2015], we run one additional pass over the data, as in Garipov et al. [2018], to compute the running mean and standard deviation of the activations for each layer of the network with w SWA weights after the training is finished, since these statistics are not collected during training. For most deep learning libraries, such as PyTorch or Tensorflow, one can typically collect these statistics by making a forward pass over the data in training mode.
The SWA procedure is summarized in Algorithm 1.
COMPUTATIONAL COMPLEXITY
The time and memory overhead of SWA compared to conventional training is negligible. During training, we need to maintain a copy of the running average of DNN weights. Note however that the memory consumption in storing a DNN is dominated by its activations rather than its weights, and thus is only slightly increased by the   
Algorithm 1 Stochastic Weight Averaging
Require: weightsŵ, LR bounds α 1 , α 2 , cycle length c (for constant learning rate c = 1), number of iterations n Ensure: w SWA w ←ŵ {Initialize weights withŵ}
w SWA ← w for i ← 1, 2, . . . , n do α ← α(i) {Calculate LR for the iteration} w ← w − α∇L i (w) {Stochastic gradient update} if mod(i, c) = 0 then n models ← i/c {Number of models} w SWA ← wSWA·nmodels+w nmodels+1
{Update average} end if end for {Compute BatchNorm statistics for w SWA weights} SWA procedure, even for large DNNs (e.g., on the order of 10%). After the training is complete we only need to store the model that aggregates the average, leading to the same memory requirements as standard training.
During training extra time is only spent to update the aggregated weight average. This operation is of the form w SWA ← w SWA · n models + w n models + 1 , and it only requires computing a weighted sum of the weights of two DNNs. As we apply this operation at most once per epoch, SWA and SGD require practically the same amount of computation. Indeed, a similar operation is performed as a part of each gradient step, and each epoch consists of hundreds of gradient steps.
OPTIMA WIDTH
Keskar et al. [2017] and Chaudhari et al. [2017] conjecture that the width of a local optimum is related to generalization. The general explanation for the importance of width is that the surfaces of train loss and test error are shifted with respect to each other and it is thus desirable to converge to the modes of broad optima, which stay approximately optimal under small perturbations. In this section we compare the solutions found by SWA and SGD and show that SWA generally leads to much wider optima.
Let w SWA and w SGD denote the weights of DNNs trained using SWA and conventional SGD, respectively. Consider the rays
w SWA (t, d) = w SWA + t · d, w SGD (t, d) = w SGD + t · d,
which follow a direction vector d on the unit sphere, starting at w SWA and w SGD , respectively. In Figure 4 we plot train loss and test error of w SWA (t, d i ) and w SGD (t, d i ) as a function of t for 10 random directions d i , i = 1, 2, . . . , 10 drawn from a uniform distribution on the unit sphere. For this visualization we use a Preactivation ResNet-164 on CIFAR-100.
First, while the loss values on train for w SGD and w SWA are quite similar (and in fact w SGD has a slightly lower train loss), the test error for w SGD is lower by 1.5% (at the converged value corresponding to t = 0). Further, the shapes of both train loss and test error curves are considerably wider for w SWA than for w SGD , suggesting that SWA indeed converges to a wider optimum: we have to step much further away from the solution found by w SWA to increase error by a given amount. We even see the error curve for SGD has an inflection point that is not present for these distances with SWA.
Notice that in Figure 4 any of the random directions from w SGD increase test error. However, we know that the direction from w SGD to w SWA would decrease test error, since w SWA has considerably lower test error than w SGD . In other words, the path from w SGD to w SWA is qualitatively different from all directions shown in Figure 4, because along this direction w SGD is far from optimal. We therefore consider the line segment connecting w SGD and w SWA : In Figure 5 we plot the train loss and test error of w(t) as a function of signed distance from w SWA for Preactivation ResNet-164 and VGG-16 on CIFAR-100.
w(t) = t · w SGD + (1 − t) · w SWA .
We can extract several key insights about w SWA and w SGD from Figure 5. First, the train loss and test error plots are indeed substantially shifted, and the point obtained by minimizing the train loss is far from optimal on test. Second, w SGD lies near the boundary of a wide flat region of the train loss. Further, the loss is very steep near w SGD .
Keskar et al. [2017] argue that the loss near sharp optima found by SGD with very large batches are actually flat in most directions, but there exist directions in which the optima are extremely steep. They conjecture that because of this sharpness the generalization performance of large batch optimization is substantially worse than that of solutions found by small batch SGD. Remarkably, in our experiments in this section we observe that there exist directions of steep ascent even for small batch optima, and that SWA provides even wider solutions (at least along random directions) with better generalization.
CONNECTION TO ENSEMBLING
Garipov et al. [2018] proposed the Fast Geometric Ensembling (FGE) procedure for training ensembles in the time required to train a single model. Using a cyclical learning rate, FGE generates a sequence of points that are close to each other in the weight space, but produce diverse predictions. In SWA instead of averaging the predictions of the models we average their weights. However, the predictions proposed by FGE ensembles and SWA models have similar properties.
Let f (·) denote the predictions of a neural network parametrized by weights w. We will assume that f is a scalar (e.g. the probability for a particular class) twice continuously differentiable function with respect to w.
Consider points w i proposed by FGE. These points are close in the weight space by design, and concentrated around their average w SWA = 1
n n i=1 w i . We denote ∆ i = w i − w SWA . Note n i=1 ∆ i = 0.
Ensembling the networks corresponds to averaging the function values
f = 1 n n i=1 f (w i ).
Consider the linearization of f at w SWA .
f (w j ) = f (w SWA ) + ∇f (w SWA ), ∆ j + O( ∆ j 2 ),
where ·, · denotes the dot product. Thus, the difference between averaging the weights and averaging the predic-
tions f − f (w SWA ) = 1 n n i=1 ∇f (w SWA ), ∆ i + O( ∆ i 2 ) = ∇f (w SWA ), 1 n n i=1 ∆ i + O(∆ 2 ) = O(∆ 2 ),
where ∆ = max n i=1 ∆ i . Note that the difference between the predictions of different perturbed networks is
f (w i ) − f (w j ) = ∇f (w SWA ), ∆ i − ∆ j + O(∆ 2 ),
and is thus of the first order of smallness, while the difference between averaging predictions and averaging weights is of the second order of smallness. Note that for the points proposed by FGE the distances between proposals are relatively small by design, which justifies the local analysis.
To analyze the difference between ensembling and averaging the weights of FGE proposals in practice, we run FGE for 20 epochs and compare the predictions of different models on the test dataset with a Preactivation ResNet-164 [He et al., 2016] on CIFAR-100. The norm of the difference between the class probabilities of consecutive FGE proposals averaged over the test dataset is 0.126. We then average the weights of the proposals and compute the class probabilities on the test dataset. The norm of difference of the probabilities for the SWA model and the FGE ensemble is 0.079, which is substantially smaller than the difference between the probabilities of consecutive FGE proposals. Further, the fraction of objects for which consecutive FGE proposals output the same labels is not greater than 87.33%. For FGE and SWA the fraction of identically labeled objects is 95.26%.
The theoretical considerations and empirical results presented in this section suggest that SWA can approximate the FGE ensemble with a single model.
CONNECTION TO CONVEX MINIMIZATION
Mandt et al. [2017] showed that under strong simplifying assumptions SGD with a fixed learning rate approximately samples from a Gaussian distribution centered at the minimum of the loss. Suppose this is the case when we run SGD with a fixed learning rate for training a DNN.
Let us denote the dimensionality of the weight space of the neural network by d. Denote the samples produced by SGD by w i , i = 1, 2, . . . , k. Assume the points w i are concentrated around the local optimumŵ. The SWA solution is given by w SWA = 1 n k i=1 w i . The points w i are samples from a multidimensional Gaussian N (ŵ, Σ) for some covariance matrix Σ defined by the curvature of the loss, batch size and the learning rate. Note that the samples from a multidimensional Gaussian are concentrated on the ellipsoid
z ∈ R d | Σ − 1 2 (z −ŵ) = √ d ,
and the probability mass for a sample to end up inside the ellipsoid nearŵ is negligible. On the other hand, w SWA is guaranteed to converge toŵ as k → ∞.
Moreover, Polyak and Juditsky [1992] showed that averaging SGD proposals achieves the best possible convergence rate among all stochastic gradient algorithms. The proof relies on the convexity of the underlying problem and in general there are no convergence guarantees if the loss function is non-convex [see e.g. Ghadimi and Lan, 2013]. While DNN loss functions are known to be nonconvex [e.g. Choromanska et al., 2015], over the trajectory of SGD these loss surfaces are approximately convex [e.g. Goodfellow et al., 2015]. However, even when the loss is locally non-convex, SWA can improve generalization. For example, in Figure 5 we see that SWA converges to a central point of the training loss.
EXPERIMENTS
We compare SWA against conventional SGD training on CIFAR-10, CIFAR-100 and ImageNet ILSVRC-2012[Russakovsky et al., 2012. We also compare to Fast Geometric Ensembling (FGE) [Garipov et al., 2018], but we note that FGE is an ensemble whereas SWA corresponds to a single model. Conventional SGD training uses a standard decaying learning rate schedule (details in the Appendix) until convergence. We found an exponentially decaying average of SGD to perform comparably to conventional SGD at convergence. We release the code for reproducing the results in this paper at https://github.com/timgaripov/swa.
CIFAR DATASETS
For the experiments on CIFAR datasets we use VGG-16 [Simonyan and Zisserman, 2014], a 164-layer Preactivation-ResNet [He et al., 2016] and Wide ResNet-28-10 [Zagoruyko and Komodakis, 2016] models. Additionally, we experiment with the recent Shake-Shake-2x64d [Gastaldi, 2017] on CIFAR-10 and PyramidNet-272 (bottleneck, α = 200) [Han et al., 2016] on CIFAR-100. All models are trained using L 2 -regularization, and VGG-16 also uses dropout.
For each model we define budget as the number of epochs required to train the model until convergence with conventional SGD training, such that we do not see improvement with SGD beyond this budget. We use the same budgets for VGG, Preactivation ResNet and Wide ResNet models as Garipov et al. [2018]. For Shake-Shake and PyramidNets we use the budgets indicated by the papers that proposed these models [Gastaldi, 2017, Han et al., 2016. We report the results of SWA training within 1, 1.25 and 1.5 budgets of epochs.
For VGG, Wide ResNet and Preactivation-ResNet models we first run standard SGD training for ≈ 75% of the training budget, and then use the weights at the last epoch as an initialization for SWA with a fixed learning rate schedule. We ran SWA for 0.25, 0.5 and 0.75 budget to complete the training within 1, 1.25 and 1.5 budgets respectively.
For Shake-Shake and PyramidNet architectures we do not report the results in one budget. For these models we use a full budget to get an initialization for the procedure, and then train with a cyclical learning rate schedule for 0.25 and 0.5 budgets. We used long cycles of small learning rates for Shake-Shake, because this architecture already involves many stochastic components.
We present the details of the learning rate schedules for each of these models in the Appendix.
For each model we also report the results of conventional SGD training, which we denote by SGD. For VGG, Preactivation ResNet and Wide ResNet we also provide the results of the FGE method with one budget reported in Garipov et al. [2018]. Note that for FGE we report the accuracy of an ensemble of 6 to 12 networks, while for SWA we report the accuracy of a single model.
We summarize the experimental results in Table 1. For all models we report the mean and standard deviation of test accuracy over 3 runs. In all conducted experi-ments SWA substantially outperforms SGD in one budget, and improves further, as we allow more training epochs. Across different architectures we see consistent improvement by ≈ 0.5% on CIFAR-10 (excluding Shake-Shake, for which SGD performance is already extremely high) and by 0.75-1.5% on CIFAR-100. Amazingly, SWA is able to achieve comparable or better performance than FGE ensembles with just one model. On CIFAR-100 SWA usually needs more than one budget to get results comparable with FGE ensembles, but on CIFAR-10 even with 1 budget SWA outperforms FGE.
IMAGENET
On ImageNet we experimented with ResNet-50, ResNet-152 [He et al., 2016] and DenseNet-161 [Huang et al., 2017]. For these architectures we used pretrained models from PyTorch.torchvision. For each of the models we ran SWA for 10 epochs with a cyclical learning rate schedule with the same parameters for all models (the details can be found in the Appendix), and report the mean and standard deviation of test error averaged over 3 runs. The results are shown in Table 2. For all 3 architectures SWA provides consistent improvement by 0.6-0.9% over the pretrained models. In this section we explore how the learning rate schedule affects the performance of SWA. We run experiments on Preactivation ResNet-164 on CIFAR-100. For all schedules we use the same initialization from a model trained for 125 epochs using the conventional SGD training. As a baseline we use a fully-trained model trained with conventional SGD for 150 epochs.
EFFECT OF THE LEARNING RATE SCHEDULE
We consider a range of constant and cyclical learning rate schedules. For cyclical learning rates we fix the cycle length to 5, and consider the pairs of base learning rate parameters (α 1 , α 2 ) ∈ {(10 −1 , 10 −3 ), (5 · 10 −2 , 5 · 10 −4 ), (10 −2 , 10 −4 ), (5 · 10 −3 , 5 · 10 −5 )}. Among the constant learning rates we consider α 1 ∈ {10 −1 , 5 · 10 −2 , 10 −2 , 10 −3 }.
We plot the test error of the SWA procedure for different learning rate schedules as a function of the number of training epochs in Figure 6.
We find that in general the more aggressive constant learning rate schedule leads to faster convergence of SWA. In our experiments we found that setting the learning rate to some intermediate value between the largest and the smallest learning rate used in the annealing scheme in conventional training usually gave us the best results. The approach is however universal and can work well with different learning rate schedules tailored for particular tasks.
DNN TRAINING WITH A FIXED LEARNING RATE
In this section we show that it is possible to train DNNs from scratch with a fixed learning rate using SWA. We run SGD with a fixed learning rate of 0.  Figure 7 illustrates the test error as a function of the number of training epochs for SWA and conventional training. The accuracy of the individual models with weights averaged by SWA stays at the level of ≈ 65% which is 16% less than the accuracy of the SWA model. These results correspond to our intuition presented in section 3.6 that SGD with a constant learning rate oscillates around the optimum, but SWA converges.
While being able to train a DNN with a fixed learning rate is a surprising property of SWA, for practical purposes we recommend initializing SWA from a model pretrained with conventional training (possibly for a reduced number of epochs), as it leads to faster and more stable convergence than running SWA from scratch.
DISCUSSION
We have presented Stochastic Weight Averaging (SWA) for training neural networks. SWA is extremely easy to implement, architecture-agnostic, and improves generalization performance at virtually no additional cost over conventional training.
There are so many exciting directions for future research. SWA does not require each weight in its average to correspond to a good solution, due to the geometry of weights traversed by the algorithm. It therefore may be possible to develop SWA for much faster convergence than standard SGD. One may also be able to combine SWA with large batch sizes while preserving generalization performance, since SWA discovers much broader optima than conventional SGD training. Furthermore, a cyclic learning rate enables SWA to explore regions of high posterior density over neural network weights. Such learning rate schedules could be developed in conjunction with stochastic MCMC approaches, to encourage exploration while still providing high quality samples. One could also develop SWA to average whole regions of good solutions, using the high-accuracy curves discovered in Garipov et al. [2018].
A better understanding of the loss surfaces for multilayer networks will help continue to unlock the potential of these rich models. We hope that SWA will inspire further progress in this area. 
A Appendix
A.1 EXPERIMENTAL DETAILS
For the experiments on CIFAR datasets (section 4.1) we used the following implementations (embedded links):
• Shake-Shake-2x64d
• PyramidNet-272
• VGG-16
• Preactivation-ResNet-164
• Wide ResNet-28-10
Models for ImageNet are from here. Pretrained networks can be found here.
SWA learning rates. For PyramidNet SWA uses a cyclic learning rate with α 1 = 0.05 and α 2 = 0.001 and cycle length 3. For VGG and Wide ResNet we used constant learning α 1 = 0.01. For ResNet we used constant learning rates α 1 = 0.01 on CIFAR-10 and 0.05 on CIFAR-100.
For Shake-Shake Net we used a custom cyclic learning rate based on the cosine annealing used when training Shake-Shake with SGD. Each of the cycles replicate the learning rates corresponding to epochs 1600 − 1700 of the standard training and the cycle length c = 100 epochs. The learning rate schedule is depicted in Figure  8 and follows the formula α(i) = 0.1 · 1 + cos π · 1600 + epoch(i) mod 100) 1800 ,
where epoch(i) is the number of data passes completed before iteration i.
For all experiments with ImageNet we used cyclic learning rate schedule with the same hyperparameters α 1 = 0.001, α 2 = 10 −5 and c = 1. SGD learning rates. For conventional SGD training we used SGD with momentum 0.9 and with an annealed learning rate schedule. For VGG, Wide ResNet and Preactivation ResNet we fixed the learning rate to α 1 for the first half of epochs (0B-0.5B), then linearly decreased the learning rate to 0.01α 1 for the next 40% of epochs (0.5B-0.9B), and then kept it constant for the last 10% of epochs (0.9B -1B). For VGG we set α 1 = 0.05, and for Preactivation ResNet and Wide ResNet we set α 1 = 0.1. For Shake-Shake Net and PyramidNets we used the cosine and piecewise-constant learning rate schedules described in Gastaldi [2017] and Han et al.
[2016] respectively.
A.2 TRAINING RESNET WITH A CONSTANT LEARNING RATE
In this section we present the experiment on training Preactivation ResNet-164 using a constant learning rate. The experimental setup is the same as in section 4.4. We set the learning rate to α 1 = 0.1 and start averaging after epoch 200. The results are presented in Figure 9. 
Figure 2 :
2Top: cyclical learning rate as a function of iteration. Bottom: test error as a function of iteration for cyclical learning rate schedule with Preactivation-ResNet-164 on CIFAR-100. Circles indicate iterations corresponding to the minimum learning rates.
Figure 3 :
3The L 2 -regularized cross-entropy train loss and test error surfaces of a Preactivation ResNet-164 on CIFAR-100 in the plane containing the first, middle and last points (indicated by black crosses) in the trajectories with (left two) cyclical and (right two) constant learning rate schedules.
Figure 4 :Figure 5 :
45(Left) Test error and (Right) L 2 -regularized cross-entropy train loss as a function of a point on a random ray starting at SWA (blue) and SGD (green) solutions for Preactivation ResNet-164 on CIFAR-100. Each line corresponds to a different random ray. L 2 -regularized cross-entropy train loss and test error as a function of a point on the line connecting SWA and SGD solutions on CIFAR-100. Left: Preactivation ResNet-164. Right: VGG-16.
Figure 6 :
6Test error as a function of training epoch for SWA with different learning rate schedules with a Preactivation ResNet-164 on CIFAR-100.
Figure 8 :
8Cyclical learning rate used for Shake-Shake as a function of iteration.
Figure 9 :
9Test error as a function of training epoch for constant (green) and decaying (blue) learning rate schedules for a Preactivation ResNet-164 on CIFAR-100. In red we average the points along the trajectory of SGD with constant learning rate starting at epoch 200.
Figure 1: Illustrations of SWA and SGD with a Preactivation ResNet-164 on CIFAR-100 1 . Left: test error surface for three FGE samples and the corresponding SWA solution (averaging in weight space). Middle and Right: test error and train loss surfaces showing the weights proposed by SGD (at convergence) and SWA, starting from the same initialization of SGD after 125 training epochs.arXiv:1803.05407v2 [cs.LG] 8 Aug 2018 
−10 
0 
10 
20 
30 
40 
50 
−10 
0 
10 
20 
30 
W1 
W2 
W3 
W SWA 
Test error (%) 
19.95 
20.64 
21.24 
22.38 
24.5 
28.49 
35.97 
50 
> 50 
−5 
0 
5 
10 
15 
20 
25 
0 
5 
10 
epoch 125 
W SGD 
W SWA 
Test error (%) 
19.62 
20.15 
20.67 
21.67 
23.65 
27.52 
35.11 
50 
> 50 
−5 
0 
5 
10 
15 
20 
25 
0 
5 
10 
epoch 125 
W SGD 
W SWA 
Train loss 
0.00903 
0.02142 
0.03422 
0.06024 
0.1131 
0.2206 
0.4391 
0.8832 
> 0.8832 
Table 1 :
1Accuracies (%) of SWA, SGD and FGE methods on CIFAR-100 and CIFAR-10 datasets for different training 
budgets. Accuracies for the FGE ensemble are from Garipov et al. [2018]. 
SWA 
DNN (Budget) 
SGD 
FGE (1 Budget) 
1 Budget 
1.25 Budgets 1.5 Budgets 
CIFAR-100 
VGG-16 (200) 
72.55 ± 0.10 
74.26 
73.91 ± 0.12 74.17 ± 0.15 74.27 ± 0.25 
ResNet-164 (150) 
78.49 ± 0.36 
79.84 
79.77 ± 0.17 80.18 ± 0.23 80.35 ± 0.16 
WRN-28-10 (200) 
80.82 ± 0.23 
82.27 
81.46 ± 0.23 81.91 ± 0.27 82.15 ± 0.27 
PyramidNet-272 (300) 
83.41 ± 0.21 
-
-
83.93 ± 0.18 84.16 ± 0.15 
CIFAR-10 
VGG-16 (200) 
93.25 ± 0.16 
93.52 
93.59 ± 0.16 93.70 ± 0.22 93.64 ± 0.18 
ResNet-164 (150) 
95.28 ± 0.10 
95.45 
95.56 ± 0.11 95.77 ± 0.04 95.83 ± 0.03 
WRN-28-10 (200) 
96.18 ± 0.11 
96.36 
96.45 ± 0.11 96.64 ± 0.08 96.79 ± 0.05 
ShakeShake-2x64d (1800) 96.93 ± 0.10 
-
-
97.16 ± 0.10 97.12 ± 0.06 
Table 2 :
2Accuracies (%) on ImageNet dataset for SWA 
and SGD with different architectures. 
SWA 
DNN 
SGD 
5 epochs 
10 epochs 
ResNet-50 
76.15 76.83 ± 0.01 76.97 ± 0.05 
ResNet-152 
78.31 78.82 ± 0.01 78.94 ± 0.07 
DenseNet-161 77.65 78.26 ± 0.09 78.44 ± 0.06 
05 on a Wide ResNet-28-10 [Zagoruyko and Komodakis, 2016] for Figure 7: Test error as a function of training epoch for constant (green) and decaying (blue) learning rate schedules for a Wide ResNet-28-10 on CIFAR-100. In red we average the points along the trajectory of SGD with constant learning rate starting at epoch 140. 300 epochs from a random initialization on CIFAR-100. We then averaged the weights at the end of each epoch from epoch 140 and until the end of training. The final test accuracy of this SWA model was 81.7.0 
50 
100 
150 
200 
250 
300 
Epochs 
15 
20 
25 
30 
35 
40 
45 
50 
Test error (%) 
SGD 
Const LR SGD 
Const LR SWA 
Karen Simonyan and Andrew Zisserman. Very deep convolutional networks for large-scale image recognition. arXiv preprint arXiv:1409.1556, 2014. Leslie N Smith and Nicholay Topin. Exploring loss function topology with cyclical learning rates. arXiv preprint arXiv:1702.04283, 2017. Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov. Dropout: A simple way to prevent neural networks from overfitting. The Journal of Machine Learning Research, 15 (1):1929-1958, 2014. Sergey Zagoruyko and Nikos Komodakis. Wide residual networks. arXiv preprint arXiv:1605.07146, 2016.
Acknowledgements. This work was supported by NSF IIS-1563887, Samsung Research, Samsung Electronics and Russian Science Foundation grant 17-11-01027. We also thank Vadim Bereznyuk for helpful comments.
Entropy-sgd: Biasing gradient descent into wide valleys. P Chaudhari, Anna Choromanska, S Soatto, Yann Le-Cun, C Baldassi, C Borgs, J Chayes, Levent Sagun, R Zecchina, International Conference on Learning Representations (ICLR. P. Chaudhari, Anna Choromanska, S. Soatto, Yann Le- Cun, C. Baldassi, C. Borgs, J. Chayes, Levent Sagun, and R. Zecchina. Entropy-sgd: Biasing gradient de- scent into wide valleys. In International Conference on Learning Representations (ICLR), 2017.
The loss surfaces of multilayer networks. Anna Choromanska, Mikael Henaff, Michael Mathieu, Gérard Ben Arous, Yann Lecun, Artificial Intelligence and Statistics. Anna Choromanska, Mikael Henaff, Michael Mathieu, Gérard Ben Arous, and Yann LeCun. The loss surfaces of multilayer networks. In Artificial Intelligence and Statistics, pages 192-204, 2015.
Sharp minima can generalize for deep nets. Laurent Dinh, Razvan Pascanu, Samy Bengio, Yoshua Bengio, International Conference on Machine Learning. Laurent Dinh, Razvan Pascanu, Samy Bengio, and Yoshua Bengio. Sharp minima can generalize for deep nets. In International Conference on Machine Learn- ing, pages 1019-1028, 2017.
Essentially no barriers in neural network energy landscape. Felix Draxler, Kambis Veschgini, Manfred Salmhofer, Fred Hamprecht, Proceedings of the 35th International Conference on Machine Learning. the 35th International Conference on Machine LearningFelix Draxler, Kambis Veschgini, Manfred Salmhofer, and Fred Hamprecht. Essentially no barriers in neu- ral network energy landscape. In Proceedings of the 35th International Conference on Machine Learning, pages 1308-1317, 2018.
Timur Garipov, Pavel Izmailov, Dmitrii Podoprikhin, P Dmitry, Andrew Gordon Vetrov, Wilson, arXiv:1802.10026Loss surfaces, mode connectivity, and fast ensembling of dnns. arXiv preprintTimur Garipov, Pavel Izmailov, Dmitrii Podoprikhin, Dmitry P Vetrov, and Andrew Gordon Wilson. Loss surfaces, mode connectivity, and fast ensembling of dnns. arXiv preprint arXiv:1802.10026, 2018.
Xavier Gastaldi, arXiv:1705.07485Shake-shake regularization. arXiv preprintXavier Gastaldi. Shake-shake regularization. arXiv preprint arXiv:1705.07485, 2017.
Stochastic first-and zeroth-order methods for nonconvex stochastic programming. Saeed Ghadimi, Guanghui Lan, SIAM Journal on Optimization. 234Saeed Ghadimi and Guanghui Lan. Stochastic first-and zeroth-order methods for nonconvex stochastic pro- gramming. SIAM Journal on Optimization, 23(4): 2341-2368, 2013.
Qualitatively characterizing neural network optimization problems. J Ian, Oriol Goodfellow, Andrew M Vinyals, Saxe, International Conference on Learning Representations. Ian J Goodfellow, Oriol Vinyals, and Andrew M Saxe. Qualitatively characterizing neural network optimiza- tion problems. International Conference on Learning Representations, 2015.
Dongyoon Han, Jiwhan Kim, Junmo Kim, arXiv:1610.02915Deep pyramidal residual networks. arXiv preprintDongyoon Han, Jiwhan Kim, and Junmo Kim. Deep pyramidal residual networks. arXiv preprint arXiv:1610.02915, 2016.
Deep residual learning for image recognition. Kaiming He, Xiangyu Zhang, Shaoqing Ren, Jian Sun, Proceedings of the IEEE conference on computer vision and pattern recognition. the IEEE conference on computer vision and pattern recognitionKaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recognition. In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770-778, 2016.
Flat minima. Sepp Hochreiter, Jürgen Schmidhuber, Neural Computation. 91Sepp Hochreiter and Jürgen Schmidhuber. Flat minima. Neural Computation, 9(1):1-42, 1997.
Densely connected convolutional networks. Gao Huang, Zhuang Liu, Q Kilian, Laurens Weinberger, Van Der Maaten, Proceedings of the IEEE conference on computer vision and pattern recognition. the IEEE conference on computer vision and pattern recognition13Gao Huang, Zhuang Liu, Kilian Q Weinberger, and Lau- rens van der Maaten. Densely connected convolutional networks. In Proceedings of the IEEE conference on computer vision and pattern recognition, volume 1, page 3, 2017.
Batch normalization: Accelerating deep network training by reducing internal covariate shift. Sergey Ioffe, Christian Szegedy, International Conference on Machine Learning. Sergey Ioffe and Christian Szegedy. Batch normaliza- tion: Accelerating deep network training by reducing internal covariate shift. In International Conference on Machine Learning, pages 448-456, 2015.
Jorge Nocedal, Mikhail Smelyanskiy, and Ping Tak Peter Tang. On large-batch training for deep learning: Generalization gap and sharp minima. Dheevatsa Nitish Shirish Keskar, Mudigere, International Conference on Learning Representations. Nitish Shirish Keskar, Dheevatsa Mudigere, Jorge No- cedal, Mikhail Smelyanskiy, and Ping Tak Peter Tang. On large-batch training for deep learning: Generaliza- tion gap and sharp minima. International Conference on Learning Representations, 2017.
Sgdr: stochastic gradient descent with restarts. Ilya Loshchilov, Frank Hutter, International Conference on Learning Representations. Ilya Loshchilov and Frank Hutter. Sgdr: stochastic gra- dient descent with restarts. International Conference on Learning Representations, 2017.
Stochastic gradient descent as approximate bayesian inference. Stephan Mandt, D Matthew, David M Hoffman, Blei, The Journal of Machine Learning Research. 181Stephan Mandt, Matthew D Hoffman, and David M Blei. Stochastic gradient descent as approximate bayesian inference. The Journal of Machine Learning Research, 18(1):4873-4907, 2017.
Variance networks: When expectation does not meet your expectations. Kirill Neklyudov, Dmitry Molchanov, Arsenii Ashukha, Dmitry Vetrov, arXiv:1803.03764arXiv preprintKirill Neklyudov, Dmitry Molchanov, Arsenii Ashukha, and Dmitry Vetrov. Variance networks: When expec- tation does not meet your expectations. arXiv preprint arXiv:1803.03764, 2018.
Acceleration of stochastic approximation by averaging. T Boris, Anatoli B Juditsky Polyak, SIAM Journal on Control and Optimization. 304Boris T Polyak and Anatoli B Juditsky. Acceleration of stochastic approximation by averaging. SIAM Journal on Control and Optimization, 30(4):838-855, 1992.
Efficient estimations from a slowly convergent robbins-monro process. David Ruppert, Cornell University Operations Research and Industrial EngineeringTechnical reportDavid Ruppert. Efficient estimations from a slowly convergent robbins-monro process. Technical report, Cornell University Operations Research and Industrial Engineering, 1988.
Averaging Weights Leads to Wider Optima and Better Generalization
https://arxiv.org/abs/1803.05407
https://arxiv.org/abs/1703.01780
Mean teachers are better role models: Weight-averaged consistency targets improve semi-supervised deep learning results
2017-03-06
Tarvainen; Antti; Valpola; Harri
<div><div>
 
 
 
 
 
 
 <p><a href="/pdf/1703.01780">Download PDF</a></p><blockquote>
 Abstract: The recently proposed Temporal Ensembling has achieved state-of-the-art
results in several semi-supervised learning benchmarks. It maintains an
exponential moving average of label predictions on each training example, and
penalizes predictions that are inconsistent with this target. However, because
the targets change only once per epoch, Temporal Ensembling becomes unwieldy
when learning large datasets. To overcome this problem, we propose Mean
Teacher, a method that averages model weights instead of label predictions. As
an additional benefit, Mean Teacher improves test accuracy and enables training
with fewer labels than Temporal Ensembling. Without changing the network
architecture, Mean Teacher achieves an error rate of 4.35% on SVHN with 250
labels, outperforming Temporal Ensembling trained with 1000 labels. We also
show that a good network architecture is crucial to performance. Combining Mean
Teacher and Residual Networks, we improve the state of the art on CIFAR-10 with
4000 labels from 10.55% to 6.28%, and on ImageNet 2012 with 10% of the labels
from 35.24% to 9.11%.
 </blockquote>
 
 
 </div><div>
 <h2>Submission history</h2><p> From: Antti Tarvainen [<a href="/show-email/84758636/1703.01780">view email</a>]
 <br />
 <strong><a href="/abs/1703.01780v1">[v1]</a></strong>
 Mon, 6 Mar 2017 09:34:56 UTC (177 KB)<br />
 <strong><a href="/abs/1703.01780v2">[v2]</a></strong>
 Mon, 12 Jun 2017 07:41:30 UTC (659 KB)<br />
 <strong><a href="/abs/1703.01780v3">[v3]</a></strong>
 Thu, 30 Nov 2017 11:14:43 UTC (522 KB)<br />
 <strong><a href="/abs/1703.01780v4">[v4]</a></strong>
 Mon, 18 Dec 2017 09:13:01 UTC (522 KB)<br />
 <strong><a href="/abs/1703.01780v5">[v5]</a></strong>
 Mon, 8 Jan 2018 08:10:09 UTC (522 KB)<br /><strong>[v6]</strong>
Mon, 16 Apr 2018 10:39:11 UTC (522 KB)<br /></p></div>||||I|||| Skip to main content
 We gratefully acknowledge support from
 the Simons Foundation and member institutions.
 > cs > arXiv:1703.01780

 Help | Advanced Search

 All fields Title Author Abstract Comments Journal reference ACM classification MSC classification Report number arXiv identifier DOI ORCID arXiv author ID Help pages Full text
 Search
 GO

 quick links

 * Login
 * Help Pages
 * About

 Computer Science > Neural and Evolutionary Computing

 arXiv:1703.01780 (cs)
 [Submitted on 6 Mar 2017 (v1), last revised 16 Apr 2018 (this version, v6)]

 Title: Mean teachers are better role models: Weight-averaged consistency targets improve semi-supervised deep learning results

 Authors: Antti Tarvainen, Harri Valpola
 Download PDF
 Abstract: The recently proposed Temporal Ensembling has achieved state-of-the-art results in several semi-supervised learning benchmarks. It maintains an exponential moving average of label predictions on each training example, and penalizes predictions that are inconsistent with this target. However, because the targets change only once per epoch, Temporal Ensembling becomes unwieldy when learning large datasets. To overcome this problem, we propose Mean Teacher, a method that averages model weights instead of label predictions. As an additional benefit, Mean Teacher improves test accuracy and enables training with fewer labels than Temporal Ensembling. Without changing the network architecture, Mean Teacher achieves an error rate of 4.35% on SVHN with 250 labels, outperforming Temporal Ensembling trained with 1000 labels. We also show that a good network architecture is crucial to performance. Combining Mean Teacher and Residual Networks, we improve the state of the art on CIFAR-10 with 4000 labels from 10.55% to 6.28%, and on ImageNet 2012 with 10% of the labels from 35.24% to 9.11%.
 Comments: In this version: Corrected hyperparameters of the 4000-label CIFAR-10 ResNet experiment. Changed Antti's contact info, Advances in Neural Information Processing Systems 30 (NIPS 2017) pre-proceedings
 Subjects: Neural and Evolutionary Computing (cs.NE) ; Machine Learning (cs.LG); Machine Learning (stat.ML)  
 Cite as: arXiv:1703.01780 [cs.NE]   
 (or arXiv:1703.01780v6 [cs.NE] for this version)   
 https://doi.org/10.48550/arXiv.1703.01780   
 Focus to learn more   
 arXiv-issued DOI via DataCite   
 

 Submission history

 From: Antti Tarvainen [view email]
 [v1] Mon, 6 Mar 2017 09:34:56 UTC (177 KB)
 [v2] Mon, 12 Jun 2017 07:41:30 UTC (659 KB)
 [v3] Thu, 30 Nov 2017 11:14:43 UTC (522 KB)
 [v4] Mon, 18 Dec 2017 09:13:01 UTC (522 KB)
 [v5] Mon, 8 Jan 2018 08:10:09 UTC (522 KB)
 [v6] Mon, 16 Apr 2018 10:39:11 UTC (522 KB)
 Full-text links:

 Download:

 * PDF
 * Other formats
 (license)
 Current browse context:
 cs.NE
 < prev | next >
 new | recent | 1703
 Change to browse by:
 cs
 cs.LG
 stat
 stat.ML

 References & Citations

 * NASA ADS
 * Google Scholar
 * Semantic Scholar

 2 blog links

 (what is this?)

 DBLP - CS Bibliography

 listing | bibtex
 Antti Tarvainen
 Harri Valpola
 a export bibtex citation Loading...

 Bibtex formatted citation

 ×
 loading...
 Data provided by:

 Bookmark

 Bibliographic Tools

 Bibliographic and Citation Tools

 Bibliographic Explorer Toggle
 Bibliographic Explorer (What is the Explorer?)
 Litmaps Toggle
 Litmaps (What is Litmaps?)
 scite.ai Toggle
 scite Smart Citations (What are Smart Citations?)
 Code, Data, Media

 Code, Data and Media Associated with this Article

 Links to Code Toggle
 Papers with Code (What is Papers with Code?)
 ScienceCast Toggle
 ScienceCast (What is ScienceCast?)
 Demos

 Demos

 Replicate Toggle
 Replicate (What is Replicate?)
 Spaces Toggle
 Hugging Face Spaces (What is Spaces?)
 Related Papers

 Recommenders and Search Tools

 Connected Papers Toggle
 Connected Papers (What is Connected Papers?)
 Core recommender toggle
 CORE Recommender (What is CORE?)
 About arXivLabs

 arXivLabs: experimental projects with community collaborators

 arXivLabs is a framework that allows collaborators to develop and share new arXiv features directly on our website.

 Both individuals and organizations that work with arXivLabs have embraced and accepted our values of openness, community, excellence, and user data privacy. arXiv is committed to these values and only works with partners that adhere to them.

 Have an idea for a project that will add value for arXiv's community? Learn more about arXivLabs and how to get involved.

 Which authors of this paper are endorsers? | Disable MathJax (What is MathJax?)
 * About
 * Help
 * Click here to contact arXiv Contact
 * Click here to subscribe Subscribe
 * Copyright
 * Privacy Policy
 * Web Accessibility Assistance

 * arXiv Operational Status
 Get status notifications via email or slack

Mean teachers are better role models: Weight-averaged consistency targets improve semi-supervised deep learning results
Authors:
Antti  Tarvainen, Aalto University, antti.tarvainen@aalto.fi
Harri  Valpola, The Curious AI Company
Abstract
The recently proposed Temporal Ensembling has achieved state-of-the-art results in several semi-supervised learning benchmarks. It maintains an exponential moving average of label predictions on each training example, and penalizes predictions that are inconsistent with this target. However, because the targets change only once per epoch, Temporal Ensembling becomes unwieldy when learning large datasets. To overcome this problem, we propose Mean Teacher, a method that averages model weights instead of label predictions. As an additional benefit, Mean Teacher improves test accuracy and enables training with fewer labels than Temporal Ensembling. Without changing the network architecture, Mean Teacher achieves an error rate of 4.35% on SVHN with 250 labels, outperforming Temporal Ensembling trained with 1000 labels. We also show that a good network architecture is crucial to performance. Combining Mean Teacher and Residual Networks, we improve the state of the art on CIFAR-10 with 4000 labels from 10.55% to 6.28%, and on ImageNet 2012 with 10% of the labels from 35.24% to 9.11%.
1 Introduction
Deep learning has seen tremendous success in areas such as image and speech recognition. In order to learn useful abstractions, deep learning models require a large number of parameters, thus making them prone to over-fitting (Figure 1a ). Moreover, adding high-quality labels to training data manually is often expensive. Therefore, it is desirable to use regularization methods that exploit unlabeled data effectively to reduce over-fitting in semi-supervised learning.
When a percept is changed slightly, a human typically still considers it to be the same object. Correspondingly, a classification model should favor functions that give consistent output for similar data points. One approach for achieving this is to add noise to the input of the model. To enable the model to learn more abstract invariances, the noise may be added to intermediate representations, an insight that has motivated many regularization techniques, such as Dropout [28] . Rather than minimizing the classification cost at the zero-dimensional data points of the input space, the regularized model minimizes the cost on a manifold around each data point, thus pushing decision boundaries away from the labeled data points (Figure 1b ).
Since the classification cost is undefined for unlabeled examples, the noise regularization by itself does not aid in semi-supervised learning. To overcome this, the Γ model [21] evaluates each data point with and without noise, and then applies a consistency cost between the two predictions. In this case, the model assumes a dual role as a teacher and a student. As a student, it learns as before; as a teacher, it generates targets, which are then used by itself as a student for learning. Since the model itself generates targets, they may very well be incorrect. If too much weight is given to the generated targets, the cost of inconsistency outweighs that of misclassification, preventing the learning of new (e) An ensemble of models gives an even better expected target. Both Temporal Ensembling and the Mean Teacher method use this approach.
information. In effect, the model suffers from confirmation bias (Figure 1c ), a hazard that can be mitigated by improving the quality of targets.
There are at least two ways to improve the target quality. One approach is to choose the perturbation of the representations carefully instead of barely applying additive or multiplicative noise. Another approach is to choose the teacher model carefully instead of barely replicating the student model. Concurrently to our research, Miyato et al. [16] have taken the first approach and shown that Virtual Adversarial Training can yield impressive results. We take the second approach and will show that it too provides significant benefits. To our understanding, these two approaches are compatible, and their combination may produce even better outcomes. However, the analysis of their combined effects is outside the scope of this paper.
Our goal, then, is to form a better teacher model from the student model without additional training.
As the first step, consider that the softmax output of a model does not usually provide accurate predictions outside training data. This can be partly alleviated by adding noise to the model at inference time [4] , and consequently a noisy teacher can yield more accurate targets (Figure 1d ). This approach was used in Pseudo-Ensemble Agreement [2] and has lately been shown to work well on semi-supervised image classification [13, 23] . Laine & Aila [13] named the method the Π model; we will use this name for it and their version of it as the basis of our experiments.
The Π model can be further improved by Temporal Ensembling [13] , which maintains an exponential moving average (EMA) 
2 Mean Teacher
To overcome the limitations of Temporal Ensembling, we propose averaging model weights instead of predictions. Since the teacher model is an average of consecutive student models, we call this the Mean Teacher method (Figure 2 ). Averaging model weights over training steps tends to produce a Both the student and the teacher model evaluate the input applying noise (η, η ) within their computation. The softmax output of the student model is compared with the one-hot label using classification cost and with the teacher output using consistency cost. After the weights of the student model have been updated with gradient descent, the teacher model weights are updated as an exponential moving average of the student weights. Both model outputs can be used for prediction, but at the end of the training the teacher prediction is more likely to be correct. A training step with an unlabeled example would be similar, except no classification cost would be applied.
more accurate model than using the final weights directly [19] . We can take advantage of this during training to construct better targets. Instead of sharing the weights with the student model, the teacher model uses the EMA weights of the student model. Now it can aggregate information after every step instead of every epoch. In addition, since the weight averages improve all layer outputs, not just the top output, the target model has better intermediate representations. These aspects lead to two practical advantages over Temporal Ensembling: First, the more accurate target labels lead to a faster feedback loop between the student and the teacher models, resulting in better test accuracy. Second, the approach scales to large datasets and on-line learning.
More formally, we define the consistency cost J as the expected distance between the prediction of the student model (with weights θ and noise η) and the prediction of the teacher model (with weights θ and noise η ).
J(θ) = E x,η ,η f (x, θ , η ) -f (x, θ, η) 2
The difference between the Π model, Temporal Ensembling, and Mean teacher is how the teacher predictions are generated. Whereas the Π model uses θ = θ, and Temporal Ensembling approximates f (x, θ , η ) with a weighted average of successive predictions, we define θ t at training step t as the EMA of successive θ weights:
θ t = αθ t-1 + (1 -α)θ t
where α is a smoothing coefficient hyperparameter. An additional difference between the three algorithms is that the Π model applies training to θ whereas Temporal Ensembling and Mean Teacher treat it as a constant with regards to optimization.
We can approximate the consistency cost function J by sampling noise η, η at each training step with stochastic gradient descent. Following Laine & Aila [13] , we use mean squared error (MSE) as the consistency cost in most of our experiments.
Table 1 : Error rate percentage on SVHN over 10 runs (4 runs when using all labels). We use exponential moving average weights in the evaluation of all our models. All the methods use a similar 13-layer ConvNet architecture. See Table 5 in the Appendix for results without input augmentation.
250 labels 73257 images 500 labels 73257 images 1000 labels 73257 images 73257 labels 73257 images GAN [25] 18.44 ± 4.8 8.11 ± 1.3 Π model [13] 6.65 ± 0.53 4.82 ± 0.17 2.54 ± 0.04 Temporal Ensembling [13] 5.12 ± 0. 
3 Experiments
To test our hypotheses, we first replicated the Π model [13] in TensorFlow [1] as our baseline. We then modified the baseline model to use weight-averaged consistency targets. The model architecture is a 13-layer convolutional neural network (ConvNet) with three types of noise: random translations and horizontal flips of the input images, Gaussian noise on the input layer, and dropout applied within the network. We use mean squared error as the consistency cost and ramp up its weight from 0 to its final value during the first 80 epochs. The details of the model and the training procedure are described in Appendix B.1.
3.1 Comparison to other methods on SVHN and CIFAR-10
We ran experiments using the Street View House Numbers (SVHN) and CIFAR-10 benchmarks [17] .
Both datasets contain 32x32 pixel RGB images belonging to ten different classes. In SVHN, each example is a close-up of a house number, and the class represents the identity of the digit at the center of the image. In CIFAR- [16] performs even better than Mean Teacher on the 1000-label SVHN and the 4000-label CIFAR-10. As discussed in the introduction, VAT and Mean Teacher are complimentary approaches. Their combination may yield better accuracy than either of them alone, but that investigation is beyond the scope of this paper. 
3.2 SVHN with extra unlabeled data
Above, we suggested that Mean Teacher scales well to large datasets and on-line learning. In addition, the SVHN and CIFAR-10 results indicate that it uses unlabeled examples efficiently. Therefore, we wanted to test whether we have reached the limits of our approach. Table 3 shows the results.
3.3 Analysis of the training curves
The training curves on Figure 3 help us understand the effects of using Mean Teacher. As expected, the EMA-weighted models (blue and dark gray curves in the bottom row) give more accurate predictions than the bare student models (orange and light gray) after an initial period.
Using the EMA-weighted model as the teacher improves results in the semi-supervised settings.
There appears to be a virtuous feedback cycle of the teacher (blue curve) improving the student (orange) via the consistency cost, and the student improving the teacher via exponential moving averaging. If this feedback cycle is detached, the learning is slower, and the model starts to overfit earlier (dark gray and light gray).
Mean Teacher helps when labels are scarce. When using 500 labels (middle column) Mean Teacher learns faster, and continues training after the Π model stops improving. On the other hand, in the all-labeled case (left column), Mean Teacher and the Π model behave virtually identically. 
3.4 Ablation experiments
To assess the importance of various aspects of the model, we ran experiments on SVHN with 250 labels, varying one or a few hyperparameters at a time while keeping the others fixed.
Removal of noise (Figures 4(a ) and 4(b)). In the introduction and Figure 1 , we presented the hypothesis that the Π model produces better predictions by adding noise to the model on both sides. But after the addition of Mean Teacher, is noise still needed? Yes. We can see that either input augmentation or dropout is necessary for passable performance. On the other hand, input noise does not help when augmentation is in use. Dropout on the teacher side provides only a marginal benefit over just having it on the student side, at least when input augmentation is in use.
Sensitivity to EMA decay and consistency weight (Figures 4(c) and 4(d))
. The essential hyperparameters of the Mean Teacher algorithm are the consistency cost weight and the EMA decay α. How sensitive is the algorithm to their values? We can see that in each case the good values span roughly an order of magnitude and outside these ranges the performance degrades quickly. Note that EMA decay α = 0 makes the model a variation of the Π model, although somewhat inefficient one because the gradients are propagated through only the student path. Note also that in the evaluation runs we used EMA decay α = 0.99 during the ramp-up phase, and α = 0.999 for the rest of the training. We chose this strategy because the student improves quickly early in the training, and thus the teacher should forget the old, inaccurate, student weights quickly. Later the student improvement slows, and the teacher benefits from a longer memory.
Decoupling classification and consistency (Figure 4 (e)). The consistency to teacher predictions may not necessarily be a good proxy for the classification task, especially early in the training. So far our model has strongly coupled these two tasks by using the same output for both. How would decoupling the tasks change the performance of the algorithm? To investigate, we changed the model to have two top layers and produce two outputs. We then trained one of the outputs for classification and the other for consistency. We also added a mean squared error cost between the output logits, and then varied the weight of this cost, allowing us to control the strength of the coupling. Looking at the results (reported using the EMA version of the classification output), we can see that the strongly coupled version performs well and the too loosely coupled versions do not. On the other hand, a moderate decoupling seems to have the benefit of making the consistency ramp-up redundant. Changing from MSE to KL-divergence (Figure 4 (f)) Following Laine & Aila [13] , we use mean squared error (MSE) as our consistency cost function, but KL-divergence would seem a more natural choice. Which one works better? We ran experiments with instances of a cost function family ranging from MSE (τ = 0 in the figure) to KL-divergence (τ = 1), and found out that in this setting MSE performs better than the other cost functions. See Appendix C for the details of the cost function family and for our intuition about why MSE performs so well.
3.5 Mean Teacher with residual networks on CIFAR-10 and ImageNet
In the experiments above, we used a traditional 13-layer convolutional architecture (ConvNet), which has the benefit of making comparisons to earlier work easy. In order to explore the effect of the model architecture, we ran experiments using a 12-block (26-layer) Residual Network [8] (ResNet) with Shake-Shake regularization [5] on CIFAR-10. The details of the model and the training procedure are described in Appendix B.2. As shown in Table 4 , the results improve remarkably with the better network architecture.
To test whether the methods scales to more natural images, we ran experiments on Imagenet 2012 dataset [22] using 10% of the labels. We used a 50-block (152-layer) ResNeXt architecture [33] , and saw a clear improvement over the state of the art. As the test set is not publicly available, we measured the results using the validation set.
4 Related work
Noise regularization of neural networks was proposed by Sietsma & Dow [26] . More recently, several types of perturbations have been shown to regularize intermediate representations effectively in deep learning. Adversarial Training [6] changes the input slightly to give predictions that are as different as possible from the original predictions. Dropout [28] zeroes random dimensions of layer outputs. Dropconnect [31] generalizes Dropout by zeroing individual weights instead of activations.
Stochastic Depth [11] drops entire layers of residual networks, and Swapout [27] generalizes Dropout and Stochastic Depth. Shake-shake regularization [5] duplicates residual paths and samples a linear combination of their outputs independently during forward and backward passes.
Several semi-supervised methods are based on training the model predictions to be consistent to perturbation. The Denoising Source Separation framework (DSS) [29] uses denoising of latent variables to learn their likelihood estimate. The Γ variant of Ladder Network [21] implements DSS with a deep learning model for classification tasks. It produces a noisy student predictions and clean teacher predictions, and applies a denoising layer to predict teacher predictions from the student predictions. The Π model [13] improves the Γ model by removing the explicit denoising layer and applying noise also to the teacher predictions. Similar methods had been proposed already earlier for linear models [30] and deep learning [2] . Virtual Adversarial Training [16] is similar to the Π model but uses adversarial perturbation instead of independent noise.
The idea of a teacher model training a student is related to model compression [3] and distillation [9] . The knowledge of a complicated model can be transferred to a simpler model by training the simpler model with the softmax outputs of the complicated model. The softmax outputs contain more information about the task than the one-hot outputs, and the requirement of representing this knowledge regularizes the simpler model. Besides its use in model compression, distillation can be used to harden trained models against adversarial attacks [18] . The difference between distillation and consistency regularization is that distillation is performed after training whereas consistency regularization is performed on training time.
Consistency regularization can be seen as a form of label propagation [34] . Training samples that resemble each other are more likely to belong to the same class. Label propagation takes advantage of this assumption by pushing label information from each example to examples that are near it according to some metric. Label propagation can also be applied to deep learning models [32] . However, ordinary label propagation requires a predefined distance metric in the input space. In contrast, consistency targets employ a learned distance metric implied by the abstract representations of the model. As the model learns new features, the distance metric changes to accommodate these features. Therefore, consistency targets guide learning in two ways. On the one hand they spread the labels according to the current distance metric, and on the other hand, they aid the network learn a better distance metric.
5 Conclusion
Temporal Ensembling, Virtual Adversarial Training and other forms of consistency regularization have recently shown their strength in semi-supervised learning. In this paper, we propose Mean Teacher, a method that averages model weights to form a target-generating teacher model. Unlike Temporal Ensembling, Mean Teacher works with large datasets and on-line learning. Our experiments suggest that it improves the speed of learning and the classification accuracy of the trained network.
In addition, it scales well to state-of-the-art architectures and large image sizes.
The success of consistency regularization depends on the quality of teacher-generated targets. If the targets can be improved, they should be. Mean Teacher and Virtual Adversarial Training represent two ways of exploiting this principle. Their combination may yield even better targets. There are probably additional methods to be uncovered that improve targets and trained models even further. a Not applied on SVHN experiments our baseline Π model and our Mean Teacher model is whether the teacher weights are identical to the student weights or an EMA of the student weights. In addition, the Π models (both the original and ours) backpropagate gradients to both sides of the model whereas Mean Teacher applies them only to the student side.
Table 6 describes the architecture of the convolutional network. We applied mean-only batch normalization and weight normalization [24] on convolutional and softmax layers. We used Leaky ReLu [15] with α = 0.1 as the nonlinearity on each of the convolutional layers.
We used cross-entropy between the student softmax output and the one-hot label as the classification cost, and the mean square error between the student and teacher softmax outputs as the consistency cost. The total cost was the weighted sum of these costs, where the weight of classification cost was the expected number of labeled examples per minibatch, subject to the ramp-ups described below.
We trained the network with minibatches of size 100. We used Adam Optimizer [12] for training with learning rate 0.003 and parameters β 1 = 0.9, β 2 = 0.999, and ε = 10 -8 . In our baseline Π model we applied gradients through both teacher and student sides of the network. In Mean teacher model, the teacher model parameters were updated after each training step using an EMA with α = 0.999. These hyperparameters were subject to the ramp-ups and ramp-downs described below.
We applied a ramp-up period of 40000 training steps at the beginning of training. The consistency cost coefficient and the learning rate were ramped up from 0 to their maximum values, using a sigmoid-shaped function e -5(1-x) 2 , where x ∈ [0, 1].
We used different training settings in different experiments. In the CIFAR-10 experiment, we matched the settings of Laine & Aila [13] as closely as possible. In the SVHN experiments, we diverged from Laine & Aila [13] to accommodate for the sparsity of labeled data. Table 7 summarizes the differences between our experiments.
B.1.1 ConvNet on CIFAR-10
We normalized the input images with ZCA based on training set statistics. We used a total cost function consisting of classification cost and three other costs: We used the dual output trick described in subsection 3.4 and Figure 4 (e) with MSE cost between logits with coefficient 0.01. This simplified other hyperparameter choices and improved the results. We used MSE consistency cost with coefficient ramping up from 0 to 100.0 during the first 5 epochs, using the same sigmoid ramp-up shape as in the experiments above. We also used an L2 weight decay with coefficient 2e-4. We used EMA decay value 0.97 (when 1000 labels) or 0.99 (when 4000 labels).
B.2.2 ResNet on ImageNet
On our ImageNet evaluation runs, we used a 152-layer ResNeXt architecture [33] consisting of 3+8+36+3 residual blocks, with 32 groups of 4 channels on the first block.
We trained the network on 10 GPUs using minibatches of 400 images, 200 of which were labeled. We sampled the images in the same way as described in the SVHN experiments above. Following [10] , we randomly augmented images using a 10 degree rotation, a crop with aspect ratio between 3/4 and 4/3 resized to 224x224 pixels, a random horizontal flip and a color jitter. We then normalized images to have channel-wise zero mean and unit variance over training data.
We trained the network using stochastic gradient descent with maximum learning rate 0.25 and Nesterov momentum 0.9. We ramped up the learning rate linearly during the first two epochs from 0.1 to 0.25. We trained for 60 epochs, decaying the learning rate with cosine annealing so that it would have reached zero after 75 epochs.
We used a total cost function consisting of classification cost and three other costs: We used the dual output trick described in subsection 3.4 and Figure 4 (e) with MSE cost between logits with coefficient 0.01. We used a KL-divergence consistency cost with coefficient ramping up from 0 to 10.0 during the first 5 epochs, using the same sigmoid ramp-up shape as in the experiments above. We also used an L2 weight decay with coefficient 5e-5. We used EMA decay value 0.9997.
MSE
0.01 0.1 0.9 0.99 KL-div cons. cost function 5% 10% 15% (f) 
B.3 Use of training, validation and test data
In the development phase of our work with CIFAR-10 and SVHN datasets, we separated 10% of training data into a validation set. We removed randomly most of the labels from the remaining training data, retaining an equal number of labels from each class. We used a different set of labels for each of the evaluation runs. We retained labels in the validation set to enable exploration of the results. In the final evaluation phase we used the entire training set, including the validation set but with labels removed.
On a real-world use case we would not possess a large fully-labeled validation set. However, this setup is useful in a research setting, since it enables a more thorough analysis of the results. To the best of our knowledge, this is the common practice when carrying out research on semi-supervised learning. By retaining the hyperparameters from previous work where possible we decreased the chance of over-fitting our results to validation labels.
In the ImageNet experiments we removed randomly most of the labels from the training set, retaining an equal number of labels from each class. For validation we used the given validation set without modifications. We used a different set of training labels for each of the evaluation runs and evaluated the results against the validation set.
C Varying between mean squared error and KL-divergence
As mentioned in subsection 3.4, we ran an experiment varying the consistency cost function between MSE and KL-divergence (reproduced in Figure 5 ). The exact consistency function we used was C τ (p, q) = Z τ D KL (p τ q τ ), where Z τ = 2 N 2 τ 2 , p τ = τ p + 1 -τ N , q τ = τ q + 1 -τ N , τ ∈ (0, 1] and N is the number of classes. Taking the Taylor expansion we get
D KL (p i q i ) = i 1 2 τ 2 N (p i -q i ) 2 + O N 2 τ 3
where the zeroth-and first-order terms vanish. Consequently,
C τ (p, q) → 1 N i (p i -q i ) 2 when τ → 0 C τ (p, q) = 2 N 2 D KL (p q) when τ = 1.
The results in Figure 5 show that MSE performs better than KL-divergence or C τ with any τ . We also tried other consistency cost weights with KL-divergence and did not reach the accuracy of MSE.
The exact reason why MSE performs better than KL-divergence remains unclear, but the form of C τ may help explain it. Modern neural network architectures tend to produce accurate but overly confident predictions [7] . We can assume that the true labels are accurate, but we should discount the confidence of the teacher predictions. We can do that by having τ = 1 for the classification cost and τ
Mean teachers are better role models: Weight-averaged consistency targets improve semi-supervised deep learning results
https://arxiv.org/abs/1703.01780
https://arxiv.org/abs/2010.12217
Exponential ReLU Neural Network Approximation Rates for Point and Edge Singularities
2020-10-23
Marcati; Carlo; Opschoor; Joost A A; Petersen; Philipp C; Schwab; Christoph
<div><div>
 
 
 
 
 
 
 <p><a href="/pdf/2010.12217">Download PDF</a></p><blockquote>
 Abstract: We prove exponential expressivity with stable ReLU Neural Networks (ReLU NNs)
in $H^1(\Omega)$ for weighted analytic function classes in certain polytopal
domains $\Omega$, in space dimension $d=2,3$. Functions in these classes are
locally analytic on open subdomains $D\subset \Omega$, but may exhibit isolated
point singularities in the interior of $\Omega$ or corner and edge
singularities at the boundary $\partial \Omega$. The exponential expression
rate bounds proved here imply uniform exponential expressivity by ReLU NNs of
solution families for several elliptic boundary and eigenvalue problems with
analytic data. The exponential approximation rates are shown to hold in space
dimension $d = 2$ on Lipschitz polygons with straight sides, and in space
dimension $d=3$ on Fichera-type polyhedral domains with plane faces. The
constructive proofs indicate in particular that NN depth and size increase
poly-logarithmically with respect to the target NN approximation accuracy
$\varepsilon&gt;0$ in $H^1(\Omega)$. The results cover in particular solution sets
of linear, second order elliptic PDEs with analytic data and certain nonlinear
elliptic eigenvalue problems with analytic nonlinearities and singular,
weighted analytic potentials as arise in electron structure models. In the
latter case, the functions correspond to electron densities that exhibit
isolated point singularities at the positions of the nuclei. Our findings
provide in particular mathematical foundation of recently reported, successful
uses of deep neural networks in variational electron structure algorithms.
 </blockquote>
 
 
 </div><div>
 <h2>Submission history</h2><p> From: Carlo Marcati [<a href="/show-email/e92242dd/2010.12217">view email</a>]
 <br /><strong>[v1]</strong>
Fri, 23 Oct 2020 07:44:32 UTC (192 KB)<br /></p></div>||||I|||| Skip to main content
 We gratefully acknowledge support from
 the Simons Foundation and member institutions.
 > math > arXiv:2010.12217

 Help | Advanced Search

 All fields Title Author Abstract Comments Journal reference ACM classification MSC classification Report number arXiv identifier DOI ORCID arXiv author ID Help pages Full text
 Search
 GO

 quick links

 * Login
 * Help Pages
 * About

 Mathematics > Numerical Analysis

 arXiv:2010.12217 (math)
 [Submitted on 23 Oct 2020]

 Title: Exponential ReLU Neural Network Approximation Rates for Point and Edge Singularities

 Authors: Carlo Marcati, Joost A. A. Opschoor, Philipp C. Petersen, Christoph Schwab
 Download PDF
 Abstract: We prove exponential expressivity with stable ReLU Neural Networks (ReLU NNs) in $H^1(\Omega)$ for weighted analytic function classes in certain polytopal domains $\Omega$, in space dimension $d=2,3$. Functions in these classes are locally analytic on open subdomains $D\subset \Omega$, but may exhibit isolated point singularities in the interior of $\Omega$ or corner and edge singularities at the boundary $\partial \Omega$. The exponential expression rate bounds proved here imply uniform exponential expressivity by ReLU NNs of solution families for several elliptic boundary and eigenvalue problems with analytic data. The exponential approximation rates are shown to hold in space dimension $d = 2$ on Lipschitz polygons with straight sides, and in space dimension $d=3$ on Fichera-type polyhedral domains with plane faces. The constructive proofs indicate in particular that NN depth and size increase poly-logarithmically with respect to the target NN approximation accuracy $\varepsilon>0$ in $H^1(\Omega)$. The results cover in particular solution sets of linear, second order elliptic PDEs with analytic data and certain nonlinear elliptic eigenvalue problems with analytic nonlinearities and singular, weighted analytic potentials as arise in electron structure models. In the latter case, the functions correspond to electron densities that exhibit isolated point singularities at the positions of the nuclei. Our findings provide in particular mathematical foundation of recently reported, successful uses of deep neural networks in variational electron structure algorithms.
 Comments: Found Comput Math (2022) 
 Subjects: Numerical Analysis (math.NA) ; Machine Learning (cs.LG)
 MSC classes: 35Q40, 41A25, 41A46, 65N30 
 Cite as: arXiv:2010.12217 [math.NA] 
 (or arXiv:2010.12217v1 [math.NA] for this version) 
 https://doi.org/10.48550/arXiv.2010.12217 
 Focus to learn more 
 arXiv-issued DOI via DataCite 
 https://doi.org/10.1007/s10208-022-09565-9 
 Related DOI: Focus to learn more 
 DOI(s) linking to related resources 
 

 Submission history

 From: Carlo Marcati [view email]
 [v1] Fri, 23 Oct 2020 07:44:32 UTC (192 KB)
 Full-text links:

 Download:

 * PDF
 * Other formats
 (license)
 Current browse context:
 math.NA
 < prev | next >
 new | recent | 2010
 Change to browse by:
 cs
 cs.LG
 cs.NA
 math

 References & Citations

 * NASA ADS
 * Google Scholar
 * Semantic Scholar
 a export bibtex citation Loading...

 Bibtex formatted citation

 ×
 loading...
 Data provided by:

 Bookmark

 Bibliographic Tools

 Bibliographic and Citation Tools

 Bibliographic Explorer Toggle
 Bibliographic Explorer (What is the Explorer?)
 Litmaps Toggle
 Litmaps (What is Litmaps?)
 scite.ai Toggle
 scite Smart Citations (What are Smart Citations?)
 Code, Data, Media

 Code, Data and Media Associated with this Article

 Links to Code Toggle
 Papers with Code (What is Papers with Code?)
 ScienceCast Toggle
 ScienceCast (What is ScienceCast?)
 Demos

 Demos

 Replicate Toggle
 Replicate (What is Replicate?)
 Spaces Toggle
 Hugging Face Spaces (What is Spaces?)
 Related Papers

 Recommenders and Search Tools

 Connected Papers Toggle
 Connected Papers (What is Connected Papers?)
 Core recommender toggle
 CORE Recommender (What is CORE?)
 About arXivLabs

 arXivLabs: experimental projects with community collaborators

 arXivLabs is a framework that allows collaborators to develop and share new arXiv features directly on our website.

 Both individuals and organizations that work with arXivLabs have embraced and accepted our values of openness, community, excellence, and user data privacy. arXiv is committed to these values and only works with partners that adhere to them.

 Have an idea for a project that will add value for arXiv's community? Learn more about arXivLabs.

 Which authors of this paper are endorsers? | Disable MathJax (What is MathJax?)
 * About
 * Help
 * Click here to contact arXiv Contact
 * Click here to subscribe Subscribe
 * Copyright
 * Privacy Policy
 * Web Accessibility Assistance

 * arXiv Operational Status
 Get status notifications via email or slack

Exponential ReLU Neural Network Approximation Rates for Point and Edge Singularities
January 11, 2022
Carlo Marcati 
Joost A A Opschoor 
Philipp C Petersen 
Christoph Schwab 
Exponential ReLU Neural Network Approximation Rates for Point and Edge Singularities
January 11, 2022Neural networksfinite element methodsexponential convergenceanalytic regularitysingularities Subject Classification: 35Q4041A2541A4665N30
We prove exponential expressivity with stable ReLU Neural Networks (ReLU NNs) in H 1 (Ω) for weighted analytic function classes in certain polytopal domains Ω, in space dimension d = 2, 3. Functions in these classes are locally analytic on open subdomains D ⊂ Ω, but may exhibit isolated point singularities in the interior of Ω or corner and edge singularities at the boundary ∂Ω. The exponential expression rate bounds proved here imply uniform exponential expressivity by ReLU NNs of solution families for several elliptic boundary and eigenvalue problems with analytic data. The exponential approximation rates are shown to hold in space dimension d = 2 on Lipschitz polygons with straight sides, and in space dimension d = 3 on Fichera-type polyhedral domains with plane faces. The constructive proofs indicate in particular that NN depth and size increase poly-logarithmically with respect to the target NN approximation accuracy ε > 0 in H 1 (Ω). The results cover in particular solution sets of linear, second order elliptic PDEs with analytic data and certain nonlinear elliptic eigenvalue problems with analytic nonlinearities and singular, weighted analytic potentials as arise in electron structure models. In the latter case, the functions correspond to electron densities that exhibit isolated point singularities at the positions of the nuclei. Our findings provide in particular mathematical foundation of recently reported, successful uses of deep neural networks in variational electron structure algorithms.
Introduction
The application of neural networks (NNs) as approximation architecture in numerical solution methods of partial differential equations (PDEs), possibly on high-dimensional parameter-and state-spaces, has attracted significant and increasing attention in recent years. We mention only [49,8,40,41,47] and the references therein. In these works, the solution of elliptic and parabolic boundary value problems is approximated by NNs which are found by minimization of a residual of the NN in the PDE.
A necessary condition for the performance of the mentioned NN-based numerical approximation methods is a high rate of approximation which is to hold uniformly over the solution set associated with the PDE under consideration. For elliptic boundary and eigenvalue problems, the function classes that weak solutions of the problems belong to are well known. Moreover, in many cases, representation systems such as splines or polynomials that achieve optimal linear or nonlinear approximation rates for functions belonging to these function classes have been identified. For functions belonging to a Sobolev or Besov type smoothness space of finite differentiation order such as continuously differentiable or Sobolev-regular functions on a bounded domain, upper bounds for the approximation rate by NNs were established for example in [52,13,53,27,50]. Here, we only mentioned results that use the ReLU activation function. Besides, for PDEs, the solutions of which have a particular structure, approximation rates of the solution that go beyond classical smoothness-based results can be achieved, such as in [9,46,24,4,21]. Again, we confine the list to publications with approximation rates for NNs with the ReLU activation function (referred to as ReLU NNs below).
In the present paper, we analyze approximation rates provided by ReLU NNs for solution classes of linear and nonlinear elliptic source and eigenvalue problems on polygonal and polyhedral domains. Mathematical results on weighted analytic regularity [15,16,14,1,5,17,34,7,30,33] imply that these classes consist of functions that are analytic with possible corner, edge, and corner-edge singularities.
Our analysis provides, for the aforementioned functions, approximation errors in Sobolev norms that decay exponentially in terms of the number of parameters M of the ReLU NNs.
Contribution
The principal contribution of this work is threefold:
1. We prove, in Theorem 4.3, a general result on the approximation by ReLU NNs of weighted analytic function classes on Q := (0, 1) d , where d = 2, 3. The analytic regularity of solutions is quantified via countably normed, analytic classes, based on weighted Sobolev spaces of Kondrat'ev type in Q, which admit corner and, in space dimension d = 3, also edge singularities. Such classes were introduced, e.g., in [5,1,14,15,16,7] and in the references there. We prove exponential expression rates by ReLU NNs in the sense that for a number M of free parameters of the NNs, the approximation error is bounded, in the H 1 -norm, by C exp(−bM 1/(2d+1) ) for constants b, C > 0.
2. Based on the ReLU NN approximation rate bound of Theorem 4.3, we establish, in Section 5, approximation results for solutions of different types of PDEs by NNs with ReLU activation. Concretely, in Section 5.1.1, we study the reapproximation of solutions of nonlinear Schrödinger equations with singular potentials in space dimension d = 2, 3. We prove that for solutions which are contained in weighted, analytic classes in Ω, ReLU NNs (whose realizations are continuous, piecewise affine) with at most M free parameters yield an approximation with accuracy of the order exp(−bM 1/(2d+1) ) for some b > 0. Importantly, this convergence is in the H 1 (Ω)-norm. In Section 5.1.2, we establish the same exponential approximation rates for the eigenstates of the Hartree-Fock model with singular potential in R 3 . This result constitutes the first, to our knowledge, mathematical underpinning of the recently reported, high efficiency of various NNbased approaches in variational electron structure computations, e.g., [39,20,18]. In Section 5.2, we demonstrate the same approximation rates also for elliptic boundary value problems with analytic coefficients and analytic right-hand sides, in plane, polygonal domains Ω. The approximation error of the NNs is, again, bound in the H 1 (Ω)-norm. We also infer an exponential NN expression rate bound for corresponding traces, in H 1/2 (∂Ω) and for viscous, incompressible flow. Finally, in Section 5.3, we obtain the same asymptotic exponential rates for the approximation of solutions to elliptic boundary value problems, with analytic data, on so-called Fichera-type domains Ω ⊂ R 3 (being, roughly speaking, finite unions of tensorized hexahedra). These solutions exhibit corner, edge and corner-edge singularities.
3. The exponential approximation rates of the ReLU NNs established here are based on emulating corresponding variable grid and degree ("hp") piecewise polynomial approximations. In particular, our construction comprises tensor product hp-approximations on Cartesian products of geometric partitions of intervals. In Theorem A.25, we establish novel tensor product hp-approximation results for weighted analytic functions on Q of the form u − u hp H 1 (Q) ≤ C exp(−b 2d √ N ) for d = 1, 2, 3, where N is the number of degrees of freedom in the representation of u hp and C, b > 0 are independent of N (but depend on u). The geometric partitions employed in Q and the architectures of the constructed ReLU NNs are of tensor product structure, and generalize to space dimension d > 3. We note that hp-approximations based on non-tensor-product, geometric partitions of Q into hexahedra have been studied before e.g. in [42,43] in space dimension d = 3. There, the rate of u − u hp H 1 (Q) exp(−b 5 √ N ) was found. Being based on tensorization, the present construction of exponentially convergent, tensorized hp-approximations in Appendix A does not invoke the rather involved polynomial trace liftings in [42,43], and is interesting in its own right: the geometric and mathematical simplification comes at the expense of a slightly smaller (still exponential) rate of approximation. Moreover, we expect that this construction of u hp will allow a rather direct derivation of rank bounds for tensor structured function approximation of u in Q, generalizing results in [22,23] and extending [32] from point to edge and corner-edge singularities.
Neural network approximation of weighted analytic function classes
The proof strategy that yields the main result, Theorem 4.3, is as follows. We first establish exponential approximation rates in the Sobolev space H 1 for tensor-product, so-called hp-finite elements for weighted analytic functions. Then, we re-approximate the corresponding quasi-interpolants by ReLU NNs.
The emulation of hp-finite element approximation by ReLU NNs is fundamentally based on the approximate multiplication network formalized in [52]. Based on the approximate multiplication operation and an extension thereof to errors measured with respect to W 1,q -norms, for q ∈ [1, ∞], we established already in [36] a reapproximation theorem of univariate splines of order p ∈ N on arbitrary meshes with N ∈ N cells. There, we observed that there exists a NN that reapproximates a variable-order, free-knot spline u in the H 1 -norm up to an error of ε > 0 with a number of free parameters that scales logarithmically in ε and |u| H 1 , linearly in N and quadratically in p. We recall this result in Proposition 3.7 below.
From this, it is apparent by the triangle inequality that, in univariate approximation problems where hp-finite elements yield exponential approximation rates, also ReLU NNs achieve exponential approximation rates, (albeit with a possibly smaller exponent, because of the quadratic dependence on p, see [36,Theorem 5.12]).
The extension of this result to higher dimensions for high-order finite elements that are built from univariate finite elements by tensorization is based on the underlying compositionality of NNs. Because of that, it is possible to compose a NN implementing a multiplication of d inputs with d approximations of univariate finite elements. We introduce a formal framework describing these operations in Section 3.
We remark that for high-dimensional functions with a radial structure, of which the univariate radial profile allows an exponentially convergent hp-approximation, exponential convergence was obtained in [36,Section 6] by composing ReLU approximations of univariate splines with an exponentially convergent approximation of the Euclidean norm, obtaining exponential convergence without the curse of dimensionality.
Outline
The manuscript is structured as follows: in Section 2, in particular Section 2.2, we review the weighted function spaces which will be used to describe the weighted analytic function classes in polytopes Ω that underlie our approximation results. In Section 2.3, we present an approximation result by tensor-product hp-finite elements for functions from the weighted analytic class. A proof of this result is provided in Appendix A. In Section 3 we review definitions of NNs and a "ReLU calculus" from [9,38] whose operations will be required in the ensuing NN approximation results.
In Section 4, we state and prove the key results of the present paper. In Section 5, we illustrate our results by deducing novel NN expression rate bounds for solution classes of several concrete examples of elliptic boundary-value and eigenvalue problems where solutions belong to the weighted analytic function classes introduced in Section 2. Some of the more technical proofs of Section 5 are deferred to Appendix B. In Section 6, we briefly recapitulate the principal mathematical results of this paper and indicate possible consequences and further directions.
Setting and functional spaces
We start by recalling some general notation that will be used throughout the paper. We also introduce some tools that are required to describe two and three dimensional domains as well as the associated weighted Sobolev spaces.
Notation
For α ∈ N d 0 , define |α| := α1 + · · · + α d and |α|∞ := max{α1, . . . , α d }.
When we indicate a relation on |α| or |α|∞ in the subscript of a sum, we mean the sum over all multi-indices that fulfill that relation: e.g., for a k ∈ N0
|α|≤k = α∈N d 0 :|α|≤k .
For a domain Ω ⊂ R d , k ∈ N0 and for 1 ≤ p ≤ ∞, we indicate by W k,p (Ω) the classical L p (Ω)-based Sobolev space of order k. We write H k (Ω) = W k,2 (Ω). We introduce the norms · W 1,p mix (Ω) as For Ω = I1 × · · · × I d , with bounded intervals Ij ⊂ R, j = 1, . . . , d, H 1 mix (Ω) = H 1 (I1) ⊗ · · · ⊗ H 1 (I d ) with Hilbertian tensor products. Throughout, C will denote a generic positive constant whose value may change at each appearance, even within an equation.
The p norm, 1 ≤ p ≤ ∞, on R n is denoted by x p . The number of nonzero entries of a vector or matrix x is denoted by x 0.
Three dimensional domain.
Let Ω ⊂ R 3 be a bounded, polygonal/polyhedral domain. Let C denote a set of isolated points, situated either at the corners of Ω or in the interior of Ω (that we refer to as the singular corners in either case, for simplicity), and let E be a subset of the edges of Ω (the singular edges). Furthermore, denote by Ec ⊂ E the set of singular edges abutting at a corner c. For each c ∈ C and each e ∈ E, we introduce the following weights:
rc(x) := |x − c| = dist(x, c), re(x) := dist(x, e), ρce(x) := re(x) rc(x) for x ∈ Ω.
For ε > 0, we define edge-, corner-, and corner-edge neighborhoods:
Ω ε e := x ∈ Ω : re(x)  ε, ∀c ∈ C , Ω ε c := x ∈ Ω : rc(x)  ε, ∀e ∈ E , Ω ε ce := x ∈ Ω : rc(x)  0 small enough so that Ω ε c ∩ Ω ε c = ∅ for all c = c ∈ C and Q ε ce ∩ Ω ε ce = Ω ε e ∩ Ω ε e = ∅ for all singular edges e = e . In the sequel, we omit the dependence of the subdomains on ε. Let also Ωce, and Ω0 := Ω \ (ΩC ∪ ΩE ∪ ΩCE ). In each subdomain Ωce and Ωe, for any multi-index α ∈ N 3 0 , we denote by α the multi-index whose component in the coordinate direction parallel to e is equal to the component of α in the same direction, and which is zero in every other component. Moreover, we set α ⊥ := α − α .
Two dimensional domain.
Let Ω ⊂ R 2 be a polygon. We adopt the convention that E := ∅. For
c ∈ C, we define Q ε c := x ∈ Ω : rc(x)  0 so that Ω ε c ∩ Ω ε c = ∅ for c = c ∈ C and write Ωc = Ω ε c . Furthermore, ΩC is defined as for d = 3, and Ω0 := Ω \ ΩC.
Weighted spaces with nonhomogeneous norms
We introduce classes of weighted, analytic functions in space dimension d = 3, as arise in analytic regularity theory for linear, elliptic boundary value problems in polyhedra, in the particular form introduced in [7]. There, the structure of the weights is in terms of Cartesian coordinates which is particularly suited for the presently adopted, tensorized approximation architectures. The definition of the corresponding classes when d = 2 is analogous. For a weight exponent vector γ = {γc, γe, c ∈ C, e ∈ E}, we introduce the nonhomogeneous, weighted Sobolev norms
v J k γ (Ω) := |α|≤k ∂ α v L 2 (Ω 0 ) + c∈C |α|≤k r (|α|−γc) + c ∂ α v L 2 (Ωc) + e∈E |α|≤k r (|α ⊥ |−γe) + e ∂ α v L 2 (Ωe) + c∈C e∈Ec |α|≤k r (|α|−γc) + c ρ (|α ⊥ |−γe) + ce ∂ α v L 2 (Ωce)
where (x)+ = max{0, x}. Moreover, we define the associated function space by
J k γ (Ω; C, E) := v ∈ L 2 (Ω) : v J k γ (Ω)  0, we define the space of weighted analytic functions with nonhomogeneous norm as
J γ (Ω; C, E; C, A) := v ∈ J ∞ γ (Ω; C, E) : |α|=k ∂ α v L 2 (Ω 0 ) ≤ CA k k!, |α|=k r (|α|−γc) + c ∂ α v L 2 (Ωc) ≤ CA k k! ∀c ∈ C, |α|=k r (|α ⊥ |−γe) + e ∂ α v L 2 (Ωe) ≤ CA k k! ∀e ∈ E, |α|=k r (|α|−γc) + c ρ (|α ⊥ |−γe) + ce ∂ α v L 2 (Ωce) ≤ CA k k!
∀c ∈ C and ∀e ∈ Ec, ∀k ∈ N0 . 
Approximation of weighted analytic functions on tensor product geometric meshes
The approximation result of weighted analytic functions via NNs that we present below is based on emulating an approximation strategy of tensor product hp-finite elements. In this section, we present this hp-finite element approximation. Let I ⊂ R be an interval. A partition of I into N ∈ N intervals is a set G such that |G| = N , all elements of G are disjoint, connected, and open subsets of I, and
U ∈G U = I.
We denote, for all p ∈ N0, by Qp(G) the piecewise polynomials of degree p on G.
One specific partition of I = [0, 1] is given by the one dimensional geometrically graded grid, which for σ ∈ (0, 1/2] and ∈ N, is given by
G 1 := J k , k = 0, . . . , , where J 0 := (0, σ ) and J k := (σ −k+1 , σ −k ), k = 1, . . . , . (2.2) Theorem 2.1. Let d ∈ {2, 3} and Q := (0, 1) d . Let C = {c}
where c is one of the corners of Q and let E = Ec contain the edges adjacent to c when d = 3, E = ∅ when d = 2. Further assume given constants C f , A f > 0, and
γ = {γc : c ∈ C}, with γc > 1, for all c ∈ C if d = 2, γ = {γc, γe : c ∈ C, e ∈ E}, with γc > 3/2 and γe > 1, for all c ∈ C and e ∈ E if d = 3.
Then, there exist Cp > 0, CL > 0 such that, for every 0  0 independent of p, L, ε.
We present the proof in Subsection A.9.3 after developing an appropriate framework of hp-approximation in Section A.
Basic ReLU neural network calculus
In the sequel, we distinguish between a neural network, as a collection of weights, and the associated realization of the NN. This is a function that is determined through the weights and an activation function. In this paper, we only consider the so-called ReLU activation: where N0 := d and N1, . . . , NL ∈ N, and where A ∈ R N ×N −1 and b ∈ R N for = 1, ..., L.
: R → R : x → max{0, x}.
For a NN Φ, we define the associated realization of the NN Φ as
R(Φ) : R d → R N L : x → xL =: R(Φ)(x),
where the output xL ∈ R N L results from
x0 := x, x := (A x −1 + b ) for = 1, . . . , L − 1, xL := AL xL−1 + bL. (3.1)
Here is understood to act component-wise on vector-valued inputs, i.e., for y = (y 1 , . . . , y m ) ∈ R m , (y) := ( (y 1 ), . . . , (y m )). We call N (Φ) := d + L j=1 Nj the number of neurons of the NN Φ, L(Φ) := L the number of layers or depth, Mj(Φ) := Aj 0 + bj 0 the number of nonzero weights in the j-th layer, and M(Φ) := L j=1 Mj(Φ) the number of nonzero weights of Φ, also referred to as its size. We refer to NL as the dimension of the output layer of Φ.
Concatenation, parallelization, emulation of identity
An essential component in the ensuing proofs is to construct NNs out of simpler building blocks. For instance, given two NNs, we would like to identify another NN so that the realization of it equals the sum or the composition of the first two NNs. To describe these operations precisely, we introduce a formalism of operations on NNs below. The first of these operations is the concatenation. Proposition 3.2 (NN concatenation, [38,Remark 2.6]). Let L1, L2 ∈ N, and let Φ 1 , Φ 2 be two NNs of respective depths L1 and L2 such that N 1 0 = N 2 L 2 =: d, i.e., the input layer of Φ 1 has the same dimension as the output layer of Φ 2 .
Then, there exists a NN Φ 1 Φ 2 , called the sparse concatenation of Φ 1 and Φ 2 , such that
Φ 1 Φ 2 has L1 + L2 layers, R(Φ 1 Φ 2 ) = R(Φ 1 ) • R(Φ 2 ) and M Φ 1 Φ 2 ≤ 2 M Φ 1 + 2 M Φ 2 .
The second fundamental operation on NNs is parallelization, achieved with the following construction. Proposition 3.3 (NN parallelization, [38,Definition 2.7]). Let L, d ∈ N and let Φ 1 , Φ 2 be two NNs with L layers and with d-dimensional input each. Then there exists a NN P(Φ 1 , Φ 2 ) with d-dimensional input and L layers, which we call the parallelization of Φ 1 and Φ 2 , such that
R P Φ 1 , Φ 2 (x) = R Φ 1 (x), R Φ 2 (x) , for all x ∈ R d and M(P(Φ 1 , Φ 2 )) = M(Φ 1 ) + M(Φ 2 ).
Proposition 3.3 requires two NNs to have the same depth. If two NNs have different depth, then we can artificially enlarge one of them by concatenating with a NN that implements the identity. One possible construction of such a NN is presented next. 
∈ N there exists a NN Φ Id d,L with L(Φ Id d,L ) = L and M(Φ Id d,L ) ≤ 2dL, such that R(Φ Id d,L ) = Id R d .
Finally, we sometimes require a parallelization of NNs that do not share inputs. Proposition 3.5 (Full parallelization of NNs with distinct inputs, [9,Setting 5.2]). Let L ∈ N and let
Φ 1 = A 1 1 , b 1 1 , . . . , A 1 L , b 1 L , Φ 2 = A 2 1 , b 2 1 , . . . , A 2 L , b 2 L
be two NNs with L layers each and with input dimensions N 1 0 = d1 and N 2 0 = d2, respectively. Then there exists a NN, denoted by FP(Φ 1 , Φ 2 ), with d-dimensional input where d = (d1 + d2) and L layers, which we call the full parallelization of Φ 1 and Φ 2 , such that for all
x = (x1, x2) ∈ R d with xi ∈ R d i , i = 1, 2 R FP Φ 1 , Φ 2 (x1, x2) = R Φ 1 (x1), R Φ 2 (x2) and M(FP(Φ 1 , Φ 2 )) = M(Φ 1 ) + M(Φ 2 ). Proof. Set FP Φ 1 , Φ 2 := A 3 1 , b 3 1 , . . . , A 3 L , b 3 L
where, for j = 1, . . . , L, we define
A 3 j := A 1 j 0 0 A 2 j and b 3 j := b 1 j b 2 j .
All properties of FP Φ 1 , Φ 2 claimed in the statement of the proposition follow immediately from the construction.
Emulation of multiplication and piecewise polynomials
In addition to the basic operations above, we use two types of functions that we can approximate especially efficiently with NNs. These are high dimensional multiplication functions and univariate piecewise polynomials. We first give the result of an emulation of a multiplication in arbitrary dimension. 
0  0 depending on Cp > 0, Cv > 0, bv > 0 and (b − a) only. In addition, R (Φ v i ε 1 ) (xj) = vi(xj) for all i ∈ {1, . . . , N 1d } and j ∈ {0, . . . , Nint}, where {xj} N int j=0 are the nodes of G 1d .
Proof. Let i = 1, . . . , N 1d . For vi as in the assumption of the corollary, we have that either supp(vi) = J for a unique J ∈ G 1d or supp(vi) = J ∪ J for two neighboring intervals J, J ∈ G 1d . Hence, there exists a partition Ti of I of at most four subintervals so that vi ∈ Sp(I, Ti), where p = (p) i∈{1,...,4} . Because of this, an application of Proposition 3.7 with q = 2 and Remark 3.8 yields that for every
0  0 such that (4.2) holds. Furthermore,
M (Φ v i ε 1 ) ≤ 4C(1 + log(p)) (p + |log (ε1)|) + C 4 i=1 p(p + |log (ε1)|) p 2 + |log (ε1)| p + (1 + log(p)) (p + |log (ε1)|) .
We use p 1 + |log(ε1)| and obtain that there exists C5 > 0 such that (4.3) holds.
Emulation of functions with singularities in cubic domains by NNs
Below we state a result describing the efficiency of re-approximating continuous, piecewise tensor product polynomial functions in a cubic domain, as introduced in Theorem 2.1, by realizations of NNs with the ReLU activation function. 
H 1 (Q) ≤ ε. (4.4) Furthermore, there holds R (Φε,c) L ∞ (Q) ≤ (2 d + 1)Cc(1 + |log ε| 2d ), M(Φε,c) ≤ C(1 + |log ε| 2d+1 ), L(Φε,c) ≤ C(1 + |log ε| log(|log ε|)),
where C > 0 depends on Cp, CN 1d , Cv, Cc, bv, d, and (b − a) only.
Proof. Assume I = ∅ as otherwise there is nothing to show. Let CI ≥ 1 be such that C −1
I ≤ (b−a) ≤ CI . Let cv,max := max{ vi H 1 (I) : i ∈ {1, . . . , N 1d }} ≤ Cvε −bv , let ε1 := min{ε/(2 · d · (cv,max + 1) d · c 1), 1/2, C −1/2 I C −1 v ε bv }, and let ε2 := min{ε/(2 · ( √ d + 1) · (cv,max + 1) · c 1), 1/2}.
Construction of the neural network. Invoking Corollary 4.1 we choose, for i = 1, . . . , N 1d , NNs
Φ v i ε 1 so that R(Φ v i ε 1 ) − vi H 1 (I) ≤ Cvε1ε −bv ≤ 1. It follows that for all i ∈ {1, . . . , N 1d } R (Φ v i ε 1 ) H 1 (I) ≤ R (Φ v i ε 1 ) − vi H 1 (I) + vi H 1 (I) ≤ 1 + cv,max (4.5)
and that, by Sobolev imbedding, Then, let Φ basis be the NN defined as
R (Φ v i ε 1 ) ∞ ≤ R (Φ v i ε 1 ) − vi ∞ + vi ∞ ≤ C 1/2 I R (Φ v i ε 1 ) − vi H 1 (I) + 1 ≤ C 1/2 I Cvε1ε −bv + 1 ≤ 2.Φ basis := FP P(Φ v 1 ε 1 , . . . , Φ v N 1d ε 1 ), . . . , P(Φ v 1 ε 1 , . . . , Φ v N 1d ε 1 ) , (4.7)
where the full parallelization is of d copies of
P(Φ v 1 ε 1 , . . . , Φ v N 1d ε 1
). Note that Φ basis is a NN with d-dimensional input and dN 1d -dimensional output. Subsequently, we introduce the N d 1d matrices E (i 1 ,...,i d ) ∈ {0, 1} d×dN 1d such that, for all (i1, . . . , i d ) ∈ {1, . . . , N 1d } d ,
E (i 1 ,...,i d ) a = {a (j−1)N 1d +i j : j = 1, . . . , d}
for all a = (a1, . . . , a dN 1d ) ∈ R dN 1d .
Note that, for all (i1, . . . , i d ) ∈ {1, . . . , N 1d } d ,
R((E (i 1 ,...,i d ) , 0) Φ basis ) : (x1, . . . , x d ) → R(Φ v i j ε 1 )(xj) : j = 1, . . . , d .
Then, we set
Φε := P Π d ε 2 ,2 (E (i 1 ,...,i d ) , 0) : (i1, . . . , i d ) ∈ {1, . . . , N 1d } d Φ basis ,(4.8)
where Π d ε 2 ,2 is according to Proposition 3.6. Note that, by (4.6), the inputs of Π d ε 2 ,2 are bounded in absolute value by 2. Finally, we define  
Φε,c := ((vec(c) , 0)) Φε, where vec(c) ∈ R N d 1d is the reshaping such that, for all (i1, . . . , i d ) ∈ {1, . . . , N 1d } d (vec(c))i = ci 1 ...i d , with i = 1 + d j=1 (ij − 1)N j−1 1d .Φ i,k (x 1 , x 2 ) = R Π 2 ε2,2 ((E (i,k) , 0)) Φ basis (x 1 , x 2 ).
Approximation accuracy. Let us now analyze if Φε,c has the asserted approximation accuracy. We estimate by the triangle inequality that
Define, for all (i1, . . . , i d ) ∈ {1, . . . , N 1d } d φi 1 ...i d = d j=1 vi j ,N 1d i 1 ,...,i d =1 ci 1 ...i d φi 1 ...i d − R(Φε,c) H 1 (Q) = N 1d i 1 ,...,i d =1 ci 1 ...i d φi 1 ...i d − N 1d i 1 ,...,i d =1 ci 1 ...i d R(Φi 1 ...i d ) H 1 (Q) ≤ N 1d i 1 ,...,i d =1 |ci 1 ...i d | φi 1 ...i d − R(Φi 1 ...i d ) H 1 (Q) .
(4.10) We have that
φi 1 ...i d − R(Φi 1 ...i d ) H 1 (Q) = d j=1 vi j − R Π d ε 2 ,2 • R Φ v i 1 ε 1 , . . . , R Φ v i d ε 1 H 1 (Q)
and, by another application of the triangle inequality, we have that
φi 1 ...i d − R (Φi 1 ...i d ) H 1 (Q) ≤ d j=1 vi j − d j=1 R Φ v i j ε 1 H 1 (Q) + d j=1 R Φ v i j ε 1 − R Π d ε 2 ,2 • R Φ v i 1 ε 1 , . . . , R Φ v i d ε 1 H 1 (Q) ≤ d j=1 vi j − d j=1 R Φ v i j ε 1 H 1 (Q) + ( √ d + 1)ε2(cv,max + 1),(4.11)
where the last estimate follows from Proposition 3.6 and the chain rule:
d j=1 R Φ v i j ε 1 − R Π d ε 2 ,2 • R Φ v i 1 ε 1 , . . . , R Φ v i d ε 1 L 2 (Q) ≤ ε2 and d j=1 R Φ v i j ε 1 − R Π d ε 2 ,2 • R Φ v i 1 ε 1 , . . . , R Φ v i d ε 1 2 H 1 (Q) = d k=1 ∂ ∂x k d j=1 R Φ v i j ε 1 − ∂ ∂x k R Π d ε 2 ,2 • R Φ v i 1 ε 1 , . . . , R Φ v i d ε 1 2 L 2 (Q) = d k=1     d j=1 j =k R Φ v i j ε 1 − ∂ ∂x k R Π d ε 2 ,2 • R Φ v i 1 ε 1 , . . . , R Φ v i d ε 1     ∂ ∂x R Φ v i k ε 1 2 L 2 (Q) ≤ d k=1 ε 2 2 ∂ ∂x R Φ v i k ε 1 2 L 2 (I) ≤ dε 2 2 (cv,max + 1) 2 ,
where we used (4.5). We now use (4.6) to bound the first term in (4.11): for d = 3, we have that, for all
(i1, . . . , i d ) ∈ {1, . . . , N 1d } d , d j=1 vi j − d j=1 R Φ v i j ε 1 H 1 (Q) ≤ (vi 1 − R(Φ v i 1 ε 1 )) ⊗ d j=2 vi j H 1 (Q) + R Φ v i j ε 1 ⊗ vi 2 − R Φ v i 2 ε 1 ⊗ vi d H 1 (Q) + d−1 j=1 R(Φ v i j ε 1 ) ⊗ (vi d − R(Φ v i d ε 1 )) H 1 (Q) =: (I).
For d = 2, we end up with a similar estimate with only two terms. By the tensor product structure, it is clear that (I) ≤ dε1(cv,max + 1) d . We have from (4.10) and the considerations above that
N 1d i 1 ,...,i d =1 ci 1 ...i d φi 1 ...i d − R (Φε,c) H 1 (Q) ≤ c 1 dε1(cv,max + 1) d + ( √ d + 1)ε2(cv,max + 1) ≤ ε.
This yields (4.12).
Bound on the L ∞ norm of the neural network. As we have already shown, R (Φ v i ε 1 ) ∞ ≤ 2. Therefore, by Proposition 3.6, R (Φε) ∞ ≤ 2 d + ε2. It follows that
R (Φε,c) ∞ ≤ c 1 2 d + ε2 ≤ (2 d + 1)Cc(1 + |log ε| 2d ).
Size of the neural network. Bounds on the size and depth of Φε,c follow from Proposition 3.6 and Corollary 4.1. Specifically, we start by remarking that there exists a constant C1 > 0 depending on Cv, bv, CI and d only, such that |log(ε1)| ≤ C1(1 + |log ε|). Then, by Corollary 4.1, there exist constants C4, C5 > 0 depending on Cp, Cv, bv, (b − a), and d only such that for all i = 1, . . . , N 1d ,  Using also the fact that N 1d ≤ C(1 + |log ε| 2 ) for C > 0 independent of ε and since d ≥ 2,
L (Φ v i ε 1 ) ≤ C4(1 + |log ε|)(1 + log(1 + |log ε|)) and M (Φ v i ε 1 ) ≤ C5(1 + |log ε| 2 ).M(Φε,c) ≤ C11(1 + |log ε|) 2d+1 ,
for a constant C11 > 0 depending on Cp, Cv, bv, (b − a), d and Cc only.
Next, we state our main approximation result, which describes the approximation of singular functions in (0, 1) d by realizations of NNs. Assume furthermore that C f , A f > 0, and
γ = {γc : c ∈ C}, with γc > 1, for all c ∈ C if d = 2, γ = {γc, γe : c ∈ C, e ∈ E}, with γc > 3/2 and γe > 1, for all c ∈ C and e ∈ E if d = 3.
Then, for every
f ∈ J γ (Q; C, E; C f , A f ) and every 0  0, and γ = {γc : c ∈ C}, with γc > 1, for all c ∈ C if d = 2, γ = {γc, γe : c ∈ C, e ∈ E}, with γc > 3/2 and γe > 1, for all c ∈ C and e ∈ E if d = 3.
Then, for all f = (f1, . . . ,
fN f ) ∈ J γ (Q; C, E; C f , A f ) N f and every 0  0 independent of N f and ε.
Exponential expression rates for solution classes of PDEs
In this section, we develop Theorem 4.3 into several exponentially decreasing upper bounds for the rates of approximation, by realizations of NNs with ReLU activation, for solution classes to elliptic PDEs with singular data (such as singular coefficients or domains with nonsmooth boundary). In particular, we consider elliptic PDEs in two-dimensional general polygonal domains, in three-dimensional domains that are a union of cubes, and elliptic eigenvalue problems with isolated point singularities in the potential which arise in models of electron structure in quantum mechanics.
In each class of examples, the solution sets belong to the class of weighted analytic functions introduced in Subsection 2.2. However, the approximation rates established in Section 4 only hold on tensor product domains with singularities on the boundary. Therefore, we will first extend the exponential NN approximation rates to functions which exhibit singularities on a set of isolated points internal to the domain, arising from singular potentials of nonlinear Schrödinger operators. In Section 5.2, we demonstrate, using an argument based on a partition of unity, that the approximation problem on general polygonal domains can be reduced to that on tensor product domains and Fichera-type domains, and establish exponential NN expression rates for linear elliptic source and eigenvalue problems. In Section 5.3, we show exponential NN expression rates for classes of weighted analytic functions on twoand three-dimensional Fichera-type domains.
Nonlinear eigenvalue problems with isolated point singularities
Point singularities emerge in the solutions of elliptic eigenvalue problems, as arise, for example, for electrostatic interactions between charged particles that are modelled mathematically as point sources in R 3 . Other problems that exhibit point singularities appear in general relativity, and for electron structure models in quantum mechanics. We concentrate here on the expression rate of "ab initio" NN approximation of the electron density near isolated singularities of the nuclear potential. Via a ReLU-based partition of unity argument, an exponential approximation rate bound for a single, isolated point singularity in Theorem 5.1 is extended in Corollary 5.4 to electron densities corresponding to potentials with multiple point singularities at a priori known locations, modeling (static) molecules.
The numerical approximation in ab initio electron structure computations with NNs has been recently reported to be competitive with other established, methodologies (e.g. [39,20] and the references there). The exponential ReLU expression rate bounds obtained here can, in part, underpin competitive performances of NNs in (static) electron structure computations.
We recall that all NNs are realized with the ReLU activation function, see (3.1).
Nonlinear Schrödinger equations
Let
Ω = R d /(2Z) d , where d ∈ {2, 3}
, be a flat torus and let V : Ω → R be a potential such that V (x) ≥ V0 > 0 for all x ∈ Ω and there exists δ > 0 and AV > 0 such that
r 2+|α|−δ ∂ α V L ∞ (Ω) ≤ A |α|+1 V |α|! ∀α ∈ N d 0 ,(5.1)
where r(x) = dist(x, (0, . . . , 0)). For k ∈ {0, 1, 2}, we introduce the Schrödinger eigenproblem that consists in finding the smallest eigenvalue λ ∈ R and an associated eigenfunction u ∈ H 1 (Ω) such that
(−∆ + V + |u| k )u = λu in Ω, u L 2 (Ω) = 1. (5.2)
There holds the following approximation result.
Theorem 5.1. Let k ∈ {0, 1, 2} and (λ, u) ∈ R × H 1 (Ω)\{0} be a solution of the eigenvalue problem (5.2) with minimal λ, where V satisfies (5.1).
Then, for every 0  N − 1,
there exists a ground state ϕ1, . . . , ϕN of (5.4), solution to (5.5).
The following statement gives exponential expression rate bounds of the NN-based approximation of electronic wave functions in the vicinity of one singularity (corresponding to the location of a nucleus) of the potential. Theorem 5.3. Assume that (5.5) has N real eigenvalues λ1, . . . , λN with associated eigenfunctions ϕ1, . . . , ϕN , such that
R 3 ϕiϕj = δij. Fix k ∈ {1, . . . , M }, let R k be one of the singularities of V and let a > 0 such that |Rj − R k | > 2a for all j ∈ {1, . . . , M } \ {k}. Let Ω k be the cube Ω k = x ∈ R 3 : x − R k ∞ ≤ a .
Then there exists a NN Φε,ϕ such that R(Φε,ϕ) :
R 3 → R N , satisfies ϕi − R(Φε,ϕ)i H 1 (Ω k ) ≤ ε, ∀i ∈ {1, . . . , N }. (5.6)
In addition, as
ε → 0, R(Φε,ϕ)i L ∞ (Ω k ) = O(|log ε| 6 ) for every i = {1, . . . , N }, M(Φε,ϕ) = O(|log(ε)| 7 + N |log(ε)| 6 ), L(Φε,ϕ) = O(|log(ε)| log(|log(ε)|)).
Proof. Let C = {(0, 0, 0)} and E = ∅ and fix k ∈ {1, . . . , M }. From the regularity result in [31,Corollary 1], see also [10,11], there exist Cϕ, Aϕ, and γc > 3/2 such that (ϕ1, . . . , ϕN ) ∈ J γc (Ω k ; C, E; Cϕ, Aϕ) N .
Then, (5.6), the L ∞ bound and the depth and size bounds on the NN Φε,ϕ follow from the hp approximation result in Proposition A.25 (centered in R k by translation), from Theorem 4.2, as in Corollary 4.4.
The arguments in the preceding subsections applied to wave functions for a single, isolated nucleus modelled by the singular potential V as in (5.1) can then be extended to give upper bounds on the approximation rates achieved by realizations of NNs of the wave functions in a bounded, sufficiently large domain containing all singularities of the nuclear potential in (5.4).
Corollary 5.4. Assume that (5.5) has N real eigenvalues λ1, . . . , λN with associated eigenfunctions ϕ1, . . . , ϕN , such that
R 3 ϕiϕj = δij. Let ai, bi ∈ R, i = 1, 2, 3, and Ω = × d i=1 (ai, bi) such that {Rj} M j=1 ⊂ Ω. Then, for every 0  1. Then, for all u ∈ J γ (Ω; C, ∅) and for every 0  1 and since the maps ψi are affine, we observe that for every i ∈ {1, . . . , Np}, there exists γ such that min γ > 1 and
ui ∈ J γ ( Ωi, {(0, 0)}, ∅), because of [P2] and [P3]. Let ε1 := ε 2Np max i∈{1,...,Np} φi W 1,∞ (Ω) det J ψ i L ∞ ((0,1) 2 ) 1 + J ψ −1 i 2 2 L ∞ (Ω i ) 1/2 .
By Theorem 4.3 and by Lemma 5.21 and Theorem 5.16 in the forthcoming Subsection 5.3, there exist Np
NNs Φû i ε 1 , i ∈ {1, . . . , Np}, such that ûi − R(Φû i ε 1 ) H 1 ( Ω i ) ≤ ε1, ∀i ∈ {1, . . . , Np},(5.11)
and there exists C∞ > 0 independent of ε1 such that, for all i ∈ {1, . . . , Np} and allε ∈ (0, 1)
R(Φû î ε ) L ∞ ( Ω i ) ≤ C∞(1 + |logε| 4 ).
The NNs given by Theorem 4.3, Lemma 5.21 and Theorem 5.16, which we here denote by Φû i ε 1 for i = 1, . . . , Np, may not have equal depth. Therefore, for all i = 1, . . . , Np and suitable Li ∈ N we define Φû i
ε 1 := Φ Id 1,L 1 i Φû i ε 1 , so that the depth is the same for all i = 1, . . . , Np.
To estimate the size of the enlarged NNs, we use the fact that the size of a NN is not smaller than the depth unless the associated realization is constant. In the latter case, we could replace the NN by a NN with one nonzero weight without changing the realization. By this argument, we obtain for all i = 1, . . . , Np that
M(Φû i ε 1 ) ≤ 2 M(Φ Id 1,L i ) + 2 M( Φû i ε 1 ) ≤ C maxj=1,...,N p L( Φû j ε 1 ) + C M( Φû i ε 1 ) ≤ C maxj=1,...,N p M( Φû j ε 1 ). Furthermore, as shown in [19], there exist NNs Φ φ i , i ∈ {1, . . . , Np}, such that R(Φ φ i )(x) = φi(x), ∀x ∈ Ω, ∀i ∈ {1, . . . , Np}.
Here we use that T is a partition R 2 , so that φi is defined on all of R 2 and [19, Theorem 5.2] applies, which itself is based on [51]. Similarly to the previously handled case of Φû i ε 1 , we can assume that Φ φ i for i = 1, . . . , Np all have equal depth and that the size of Φ φ i is bounded independent of i.
Since by [P2] the mappings ψi are affine and invertible, it follows that ψ −1 i is affine for every
i ∈ {1, . . . , Np}. Thus, there exist NNs Φ ψ −1 i , i ∈ {1, . . . , Np}, of depth 1, such that R(Φ ψ −1 i )(x) = ψ −1 i (x), ∀x ∈ Ωi, ∀i ∈ {1, .
. . , Np}. 
), 0) P Π 2 ε × ,M × (ε 1 ) P(Φû i ε 1 Φ ψ −1 i , Φ Id 1,L Φ φ i ) Np i=1 , (5.13) where L ∈ N is such that L(Φû 1 ε 1 Φ ψ −1 1 ) = L(Φ Id 1,L Φ φ 1 ), which yields that M(Φ Id 1,L ) ≤ C L(Φû 1 ε 1 Φ ψ −1 1 ).
Approximation accuracy. By (5.13), we have for all x ∈ Ω,
R(Φε,u)(x) = Np i=1 R(Π 2 ε × ,M × (ε 1 ) ) R(Φû i ε 1 Φ ψ −1 i )(x), R(Φ φ i )(x) .
Therefore,
u − R(Φε,u) H 1 (Ω) ≤ u − Np i=1 φi R(Φû i ε 1 Φ ψ −1 i ) H 1 (Ω) + Np i=1 R(Π 2 ε × ,M × (ε 1 ) ) R(Φû i ε 1 Φ ψ −1 i ), φi − φi R(Φû i ε 1 Φ ψ −1 i ) H 1 (Ω) = (I) + (II). (5.14)
We start by considering term (I). For each i ∈ {1, . . . , Np}, thanks to (5.11), there holds, with J ψ −1 i 2 2
denoting the square of the matrix 2-norm of the Jacobian of ψ −1 i ,
u − R(Φû i ε 1 Φ ψ −1 i ) H 1 (Ω i ) = ûi • ψ −1 i − R(Φû i ε 1 ) • ψ −1 i H 1 (Ω i ) = Ω i |ûi| 2 + J ψ −1 i ∇ ûi − R(Φû i ε 1 ) 2 2 det J ψ i dx 1/2 ≤ ε1 det J ψ i L ∞ ( Ω i ) + det J ψ i L ∞ ( Ω i ) J ψ −1 i 2 2 L ∞ (Ω i ) 1/2 ≤ ε2 := ε1 max i det J ψ i L ∞ ( Ω i ) + det J ψ i L ∞ ( Ω i ) J ψ −1 i 2 2 L ∞ (Ω i ) 1/2 . (5.15) By [35, Theorem 2.1], (I) ≤ Npε2 max i∈{1,...,Np} φi W 1,∞ (Ω) ≤ ε 2 . (5.16)
We now consider term (II) in (5.14). There holds, by Theorem 4.3 and (5.12),
R(Φû i ε 1 Φ ψ −1 i ) L ∞ (Ω i ) = R(Φû i ε 1 ) L ∞ ( Ω i ) ≤ C∞(1 + |log ε1| 4 ) for all i ∈ {1, . . . , Np}. Furthermore, by [P1], φi(x) = 0 for all x ∈ Ω \ Ωi and, by Proposition 3.6, R(Π 2 ε × ,M × (ε 1 ) ) R(Φû i ε 1 Φ ψ −1 i )(x), φi(x) = 0, ∀x ∈ Ω \ Ωi. From (5.15), we also have | R(Φû i ε 1 Φ ψ −1 i )| H 1 (Ω i ) ≤ |u| H 1 (Ω i ) + u − R(Φû i ε 1 Φ ψ −1 i ) H 1 (Ω i ) ≤ 1 + |u| H 1 (Ω i ) .
Hence,
(II) = Np i=1 R(Π 2 ε × ,M × (ε 1 ) ) R(Φû i ε 1 Φ ψ −1 i ), φi − φi R(Φû i ε 1 Φ ψ −1 i ) H 1 (Ω i ) ≤ Np i=1 R(Π 2 ε × ,M × (ε 1 ) )(a, b) − ab W 1,∞ ([−M × (ε 1 ),M × (ε 1 )] 2 ) · |Ω| 1/2 + | R(Φû i ε 1 Φ ψ −1 i )| H 1 (Ω i ) + |φi| H 1 (Ω i ) ≤ Npε× |Ω| 1/2 + 1 + |u| H 1 (Ω i ) + |Ω| 1/2 max i=1,...,Np φi W 1,∞ (Ω) ≤ ε 2 .
(5.17)
The asserted approximation accuracy follows by combining (5.14), (5.16), and (5.17).
Size of the neural network.
To bound the size of the NN, we remark that Np and the sizes of Φ ψ −1 i and of Φ φ i only depend on the domain Ω. Furthermore, there exist constants CΩ,i, i = 1, 2, 3, that depend only on Ω and u such that
|log ε1| ≤ CΩ,1(1 + |log ε|), |log ε×| ≤ CΩ,2(1 + |log ε|), |log M×(ε1)| ≤ CΩ,3(1 + log(1 + |log ε|)). (5.18)
From Theorem 4.3 and Proposition 3.6, in addition, there exist constants C L u , C M u , C× > 0 such that, for all 0  1. For any connected component Γ of ∂Ω, let Γ > 0 be the length of Γ, such that there exists a continuous, piecewise affine parametrization θ : [0, Γ] → R 2 : t → θ(t) of Γ with finitely many affine linear pieces and d dt θ 2 = 1 for almost all t ∈ [0, Γ]. Then, for all u ∈ J γ (Ω; C, ∅) and for all 0  0 such that for all v ∈ H 1 (Ω) it holds T v H 1/2 (Γ) ≤ CΓ v H 1 (Ω) , and without loss of generality we may assume CΓ ≥ 1.
Next, for any ε ∈ (0, 1), let Φ ε/C Γ ,u be as given by Theorem 5.6. Define Φ ε,u,θ := Φ ε/C Γ ,u Φ θ . It follows that
T u − R(Φ ε,u,θ ) • θ −1 H 1/2 (Γ) = T u − R(Φ ε/C Γ ,u ) H 1/2 (Γ) ≤ CΓ u − R(Φ ε/C Γ ,u ) H 1 (Ω) ≤ ε.
The bounds on its depth and size follow directly from Proposition 3.2, Theorem 5.6, and the fact that the depth and size of Φ θ are independent of ε. This finishes the proof.
Remark 5.8. The exponent 5 in the bound on the NN size M(Φ ε,u,θ ) in Corollary 5.7 is likely not optimal, due to it being transferred from the NN rate in Ω.
The proof of Theorem 5.6 established exponential expressivity of realizations of NNs with ReLU activation for the analytic class J γ (Ω; C, ∅) in Ω. This implies that realizations of NNs can approximate, with exponential expressivity, solution classes of elliptic PDEs in polygonal domains Ω. We illustrate this by formulating concrete results for three problem classes: second order, linear, elliptic source and eigenvalue problems in Ω, and viscous, incompressible flow. To formulate the results, we specify the assumptions on Ω. Definition 5.9 (Linear, second order, elliptic divergence-form differential operator with analytic coefficients). Let d ∈ {2, 3} and let Ω ⊂ R d be a bounded domain. Let the coefficient functions aij, bi, c : Ω → R be real analytic in Ω, and such that the matrix function A = (aij) 1≤i,j≤d : Ω → R d×d is symmetric and uniformly positive definite in Ω. With these functions, we define the linear, second order, elliptic divergence-form differential operator L acting on w ∈ C ∞ 0 (Ω) via (summation over repeated indices i, j ∈ {1, . . . , d})  Proof. The proof is obtained by verifying weighted, analytic regularity of solutions. By [2,Theorem 3.1] there exists γ such that min γ > 1 and u ∈ J γ (Ω; C, ∅). Then, the application of Theorem 5.6 concludes the proof.
(Lw)(x) := −∂i(aij(x)∂jw(x)) + bj(x)∂jw(x) + c(x)w(x) , x ∈ Ω .
Next, we address NN expression rates for eigenfunctions of (L, B). The analytic regularity of solutions u in the proof of Theorem 5.6 also holds for certain nonlinear, elliptic PDEs. We illustrate it for the velocity field of viscous, incompressible flow in Ω. Proof. The velocity fields of Leray solutions of the Navier-Stokes equations in Ω satisfy the weighted, analytic regularity u ∈ J γ (Ω; C, ∅) 2 , with min γ > 1, see [33]. Then, the application of Theorem 5.6 concludes the proof.
Elliptic PDEs in Fichera-type polyhedral domains
Fichera-type polyhedral domains Ω ⊂ R 3 are, loosely speaking, closures of finite, disjoint unions of (possibly affinely mapped) axiparallel hexahedra with ∂Ω Lipschitz. In Fichera-type domains, analytic regularity of solutions of linear, elliptic boundary value problems from acoustics and linear elasticity in displacement formulation has been established in [7]. As an example of a boundary value problem covered by [7] and our theory, consider Ω := (−1, 1) d \ (−1, 0] d for d = 2, 3, displayed for d = 3 in Figure 2.
We recall that all NNs are realized with the ReLU activation function, see (3.1). We introduce the setting for elliptic problems with analytic coefficients in Ω. Note that the boundary of Ω is composed of 6 edges when d = 2 and of 9 faces when d = 3.  The following extension lemma will be useful for the approximation of the solution to (5.29) by NNs. We postpone its proof to Appendix B.2. Proof. By Lemma 5.15, we extend the function u to a functionũ such that
u ∈ W 1,1 mix ((−1, 1) d ) andũ|Ω = u.
Note that, by the stability of the extension, there exists a constant Cext > 0 independent of u such that By Theorem A.25 exist Cp > 0, C N 1d > 0, C N int > 0, Cṽ > 0, C c > 0, and bṽ > 0 such that, for all 0  0 is independent of (p1, p2, p3), (s1, s2, s3), and v.
Proof. Let (p1, p2, p3) ∈ N 3 , and (s1, s2, s3) ∈ N 3 , be as in the statement of the lemma. Also, let i ∈ {1, 2, 3} and {j, k} = {1, 2, 3} \ {i}. By Lemma A.4, there holds
∂x i (v − Πpv) 2 L 2 ( K) ≤ Cappx1 Ψp 1 ,s 1 α 1 ,α 2 ≤1 ∂ (s 1 +1,α 1 ,α 2 ) v 2 L 2 ( K) + Ψp 2 ,s 2 α 1 ,α 2 ≤1 ∂ (α 1 ,s 2 +1,α 2 ) v 2 L 2 ( K) + Ψp 3 ,s 3 α 1 ,α 2 ≤1 ∂ (α 1 ,α 2 ,s 3 +1) v 2 L 2 ( K) . (A.17)
With a Cappx1 > 0 independent of (p1, p2, p3), (s1, s2, s3), and v. Let now, vi :
I 2 → R such that vi(xj, x k ) = I v(x1, x2, x3)dxi.
We denoteṽ := v − vi and, remarking that ∂x i vi = ∂x i Πpvi = 0, we apply (A.17) toṽ, so that
∂x i (v − Πpv) 2 L 2 ( K) ≤ C Ψp 1 ,s 1 α 1 ,α 2 ≤1 ∂ (s 1 +1,α 1 ,α 2 )ṽ 2 L 2 ( K) + Ψp 2 ,s 2 α 1 ,α 2 ≤1 ∂ (α 1 ,s 2 +1,α 2 )ṽ 2 L 2 ( K) + Ψp 3 ,s 3 α 1 ,α 2 ≤1 ∂ (α 1 ,α 2 ,s 3 +1)ṽ 2 L 2 ( K) . (A.18)
By the Poincaré inequality, it holds for all α1 ∈ {0, 1} that
∂ s j +1 x j ∂ α 1 x kṽ 2 L 2 ( K) ≤ C ∂x i ∂ s j +1 x j ∂ α 1 x k v 2 L 2 ( K) and ∂ α 1 x j ∂ s k +1 x kṽ 2 L 2 ( K) ≤ C ∂x i ∂ α 1 x j ∂ s k +1 x k v 2 L 2 ( K) .
Using the fact that ∂x iṽ = ∂x i v in the remaining terms of (A.18) concludes the proof.
A.4.1 One dimensional estimate
The following result is a consequence of, e.g., [ 
v ∈ H s+1 (J k ) h −2 v − π k p v 2 L 2 (J k ) + ∇(v − π k p v) 2 L 2 (J k ) ≤ Cτ 2(s+1) σ Ψp,sh 2(min{γ−1,s}) |x| (s+1−γ) + v (s+1) 2 L 2 (J k ) (A.19) where h = |J k | σ −k .
Proof. From [43,Lemma 8.1], there exists C > 0 independent of p, k, s, and v such that
h −2 v − π k p v 2 L 2 (J k ) + ∇(v − π k p v) 2 L 2 (J k ) ≤ CΨp,sh 2s v (s+1) 2 L 2 (J k ) .
In addition, for all k = 1, . . . , , there holds x| J k ≥ σ 1−σ h. Hence, for all γ  0 such that for all v ∈ J 2 β (K)
α∈N 2 0 :|α|≤1 ∂ α (π 0 1 ⊗ π 0 1 )v L 2 (K) ≤ C1   v H 1 (K) + α∈N 2 0 :|α|=2 σ (β−1) r 2−β ∂ α v L 2 (K)   . (A.20)
and
α∈N 2 0 :|α|≤1 σ − (1−|α|) ∂ α (v − (π 0 1 ⊗ π 0 1 )v) L 2 (K) ≤ C2σ (β−1) α∈N 2 0 :|α|=2 r 2−β ∂ α v L 2 (K) . (A.21)
Proof. Denote by ci, i = 1, . . . , 4 the corners of K and by ψi, i = 1, . . . , 4 the bilinear functions such that ψi(cj) = δij. Then,
(π 0 1 ⊗ π 0 1 )v = 4 i=1 v(ci)ψi.
Therefore, writing h = σ , we have
(π 0 1 ⊗ π 0 1 )v L 2 (K) ≤ i=1,...,4 |v(ci)| ψi L 2 (K) ≤ 4 v L ∞ (K) |K| 1/2 ≤ 4h v L ∞ (K) . (A.22)
With the imbedding J 2 β ((0, 1) 2 ) → L ∞ ((0, 1) 2 ) which is valid for β > 1 (which follows e.g. from Lemma A.22 and W 1,1 mix ((0, 1) 2 ) → L ∞ ((0, 1) 2 )), a scaling argument gives
h 2 v 2 L ∞ (K) ≤ Ch 2   h −2 v 2 L 2 (K) + |v| 2 H 1 (K) + |α|=2 h 2β−2 r 2−β ∂ α v 2 L 2 (K)   ,
so that we obtain
(π 0 1 ⊗ π 0 1 )v 2 L 2 (K) ≤ C   v 2 L 2 (K) + h 2 |v| 2 H 1 (K) + |α|=2 h 2β r 2−β ∂ α v 2 L 2 (K)   . (A.23)
For any |α| = 1, denoting v0 = v(0, 0) and using the fact that (π 0 1 ⊗ π 0 1 )v0 = v0 hence ∂ α (π 0 1 ⊗ π 0 1 )v0 = 0,
∂ α (π 0 1 ⊗π 0 1 )v L 2 (K) = ∂ α (π 0 1 ⊗π 0 1 )(v−v0) L 2 (K) ≤ i=1,...,4 |(v−v0)(ci)| ∂ α ψi L 2 (K) ≤ C v−v0 L ∞ (K) . (A.24) With the imbedding J 2 β ((0, 1) 2 ) → L ∞ ((0, 1) 2 )
, Poincaré's inequality, and rescaling we obtain
∂ α (π 0 1 ⊗ π 0 1 )v 2 L 2 (K) ≤ C   |v| 2 H 1 (K) + |α|=2 h 2β−2 r 2−β ∂ α v 2 L 2 (K)   ,
which finishes the proof of (A.20). To prove (A.21), note that by the Sobolev imbedding of W 2,1 (K) into H 1 (K) and by scaling, we have
|α|≤1 h |α|−1 ∂ α (v − (π 0 1 ⊗ π 0 1 )v) L 2 (K) ≤ C |α|≤2 h |α|−2 ∂ α (v − (π 0 1 ⊗ π 0 1 )v) L 1 (K) .
By classical interpolation estimates [6,Theorem 4.4.4], we additionally conclude that
|α|≤1 h |α|−2 ∂ α (v − (π 0 1 ⊗ π 0 1 )v) L 1 (K) ≤ C|v| W 2,1 (K) .
Using the Cauchy-Schwarz inequality,
|α|≤1 h |α|−1 ∂ α (v − (π 0 1 ⊗ π 0 1 )v) L 2 (K) ≤ C |α|=2 ∂ α v L 1 (K) ≤ C |α|=2 r −2+β L 2 (K) r 2−β ∂ α v L 2 (K) ≤ C |α|=2 h β−1 r 2−β ∂ α v L 2 (K)
where we also have used, in the last step, the facts that r(x) ≤ √ 2h for all x ∈ K and that β > 1.
A.5 Interior estimates
The following lemmas give the estimate of the approximation error on the elements not belonging to edge or corner layers. For d = 3, all ∈ N, all k1, k2, k3 ∈ {0, . . . , } and all K = J k 1 × J k 2 × J k 3 , we denote, by h the length of K in the direction parallel to the closest singular edge, and by h ⊥,1 and h ⊥,2 the lengths of K in the other two directions. If an element has multiple closest singular edges, we choose one of those and consider it as "closest edge" for all points in that element. When considering functions from J d γ (Q), γe will refer to the weight of this closest edge. Similarly, we denote by ∂ (resp. ∂ ⊥,1 and ∂ ⊥,2 ) the derivatives in the direction parallel (resp. perpendicular) to the closest singular edge. Lemma A.9. Let d = 3, ∈ N and K = J k 1 × J k 2 × J k 3 for 0  0 dependent only on σ, Cappx2, Cv and A > 0 dependent only on σ, Av such that for all 1 ≤ s ≤ p
∂ (v − Π K p v) 2 L 2 (K) ≤ CΨp,sA 2s+6 (d K c ) 2 + (d K c ) 2(γc−1) ((s + 3)!) 2 , (A.25)
where ∂ is the derivative in the direction parallel to the closest singular edge.
Proof. We write da = d K a , a ∈ {c, e}. There holds
d 2 c = σ 1 − σ 2 (h 2 + h 2 ⊥,1 + h 2 ⊥,2 ), d 2 e = σ 1 − σ 2 (h 2 ⊥,1 + h 2 ⊥,2 ). Denotingv = v • Φ −1 K and Πpv = Π K p v • Φ −1 K = Πp(v • ΦK )
, using the result of Lemma A.6 and rescaling, we have
∂ (v − Πpv) 2 L 2 ( K) ≤ Cappx2Ψp,s h h ⊥,1 h ⊥,2   α 1 ,α 2 ≤1 h 2s h 2α 1 ⊥,1 h 2α 2 ⊥,2 ∂ s+1 ∂ α 1 ⊥,1 ∂ α 2 ⊥,2 v 2 L 2 (K) + α 1 ≤1 h 2s+2 ⊥,1 h 2α 1 ⊥,2 ∂ ∂ s+1 ⊥,1 ∂ α 1 ⊥,2 v 2 L 2 (K) + α 1 ≤1 h 2α 1 ⊥,1 h 2s+2 ⊥,2 ∂ ∂ α 1 ⊥,1 ∂ s+1 ⊥,2 v 2 L 2 (K)   = Cappx2Ψp,s h h ⊥,1 h ⊥,2 (I) + (II) + (III) . (A.26) Denote Kc = K ∩ Qc, Ke = K ∩ Qe, Kce = K ∩ Qce, and K0 = K ∩ Q0. Furthermore, we indicate (I)c = α 1 ,α 2 ≤1 h 2s h 2α 1 ⊥,1 h 2α 2 ⊥,2 ∂ s+1 ∂ α 1 ⊥,1 ∂ α 2 ⊥,2 v 2 L 2 (Kc) ,
and do similarly for the other terms of the sum (II) and (III) and the other subscripts e, ce, 0. Remark also that r i| K ≥ di, i ∈ {c, e}, and that for a, b ∈ R holds r a c r b e = r a+b c ρ b ce . We will also write γ = γc − γe. We start by considering the term (I)ce. Let α1 = α2 = 1; then,
h 2s h 2 ⊥,1 h 2 ⊥,2 ∂ s+1 ∂ ⊥,1 ∂ ⊥,2 v 2 L 2 (Kce) ≤ τ 2s+4 σ d 2s c d 4 e ∂ s+1 ∂ ⊥,1 ∂ ⊥,2 v 2 L 2 (Kce) ≤ τ 2s+4 σ d 2 γ−2 c d 2γe e r s+3−γc c ρ 2−γe ce ∂ s+1 ∂ ⊥,1 ∂ ⊥,2 v 2 L 2 (Kce) , where τσ is as in (A.6). Furthermore, if α1 + α2 ≤ 1 and s + 1 + α1 + α2 − γc ≥ 0, h 2s h 2α 1 ⊥,1 h 2α 2 ⊥,2 ∂ s+1 ∂ α 1 ⊥,1 ∂ α 2 ⊥,2 v 2 L 2 (Kce) ≤ τ 2s+2(α 1 +α 2 ) σ d 2s c d 2(α 1 +α 2 ) e ∂ s+1 ∂ α 1 ⊥,1 ∂ α 2 ⊥,2 v 2 L 2 (Kce) ≤ τ 2s+2(α 1 +α 2 ) σ d 2γc−2 c r s+1+α 1 +α 2 −γc c ∂ s+1 ∂ α 1 ⊥,1 ∂ α 2 ⊥,2 v 2 L 2 (Kce)
, where we have also used de ≤ dc. Therefore,
(I)ce ≤ τ 2s+4 σ d 2γc−2 c α 1 ,α 2 ≤1 r s+1+α 1 +α 2 −γc c ρ (α 1 +α 2 −γe) + ce ∂ s+1 ∂ α 1 ⊥,1 ∂ α 2 ⊥,2 v 2 L 2 (Kce) .
If s + 1 + α1 + α2 − γc  0 dependent only on Cv and σ and A > 0 dependent only on Av and σ such that
(I) ≤ CA 2s+6 ((s + 3)!) 2 d 2 c + d 2γc−2 c . (A.27)
We then apply the same argument to the terms (II) and (III). Indeed,
(II)ce = α 1 ≤1 h 2s+2 ⊥,1 h 2α 1 ⊥,2 ∂ ∂ s+1 ⊥,1 ∂ α 1 ⊥,2 v 2 L 2 (Kce) ≤ τ 2s+4 σ α 1 ≤1 d 2s+2+2α 1 e ∂ ∂ s+1 ⊥,1 ∂ α 1 ⊥,2 v 2 L 2 (Kce) ≤ τ 2s+4 σ α 1 ≤1 d 2 γ−2 c d 2γe e r s+2+α 1 −γc c ρ s+1+α 1 −γe ce ∂ ∂ s+1 ⊥,1 ∂ α 1 ⊥,2 v 2 L 2 (Kce)
and the estimate for (III)ce follows by exchanging h ⊥,1 and ∂ ⊥,1 with h ⊥,2 and ∂ ⊥,2 in the inequality above. The estimates for (II)c,e,0 and (III)c,e,0 can be obtained as for (I)c,e,0:
(II)c ≤ τ 2s+4 σ α 1 ≤1 d 2γc−2 c r s+2+α 1 −γc c ∂ ∂ s+1 ⊥,1 ∂ α 1 ⊥,2 v 2 L 2 (Kc) , (II)e ≤ τ 2s+4 σ α 1 ≤1 d 2γe e r s+1+α 1 −γe e ∂ ∂ s+1 ⊥,1 ∂ α 1 ⊥,2 v 2 L 2 (Ke) , (II)0 ≤ τ 2s+4 σ α 1 ≤1 d 2s+2 e ∂ ∂ s+1 ⊥,1 ∂ α 1 ⊥,2 v 2 L 2 (K 0 ) . Therefore, we have (II), (III) ≤ CA 2s+6 (d 2 c + d 2γc−2 c )((s + 3)!) 2 . (A.28)
We obtain, from (A.26), (A.27), and (A.28) that there exists C > 0 (dependent only on σ, Cappx2, Cv and A > 0 (dependent only on σ, Av) such that
∂ (v − Πpv) 2 L 2 ( K) ≤ C h h ⊥,1 h ⊥,2 Ψp,sA 2s+6 (d 2 c + d 2γc−2 c )((s + 3)!) 2 .
Considering that
∂ (v − Πpv) 2 L 2 (K) ≤ h ⊥,1 h ⊥,2 h ∂ (v − Πpv) 2 L 2 ( K)
completes the proof.
Lemma A.10. Let d = 3, ∈ N and K = J k 1 × J k 2 × J k 3 for 0  0 dependent only on σ, Cappx2, Cv and A > 0 dependent only on σ, Av such that for all p ∈ N and all 1 ≤ s ≤ p
∂ ⊥,1 (v − Π K p v) 2 L 2 (K) + ∂ ⊥,2 (v − Π K p v) 2 L 2 (K) ≤ CΨp,sA 2s+6 (d K c ) 2(γc−1) + (d K c ) 2(γe−1) ((s + 3)!) 2 , (A.29)
where ∂ ⊥,1 , ∂ ⊥,2 are the derivatives in the directions perpendicular to the closest singular edge.
Proof. The proof follows closely that of Lemma A.9 and we use the same notation. From Lemma A.6 and rescaling, we have As before, we will write γ = γc − γe. We start by considering the term (I)ce. When α1 = 1,
∂ ⊥,1 (v − Πpv) 2 L 2 ( K) ≤ Cappx2Ψp,s h ⊥,1 h h ⊥,2   α 1 ≤1 h 2s+2 h 2α 1 ⊥,2 ∂ s+1 ∂ ⊥,1 ∂ α 1 ⊥,2 v 2 L 2 (K) + α 1 ,α 2 ≤1 h 2α 1 h 2s ⊥,1 h 2α 2 ⊥,2 ∂ α 1 ∂ s+1 ⊥,1 ∂ α 2 ⊥,2 v 2 L 2 (K) + α 1 ≤1 h 2α 1 h 2s+2 ⊥,2 ∂ α 1 ∂ ⊥,1 ∂ s+1 ⊥,2 v 2 L 2 (K)   = Cappx2Ψp,s h ⊥,1 h h ⊥,2(Ih 2s+2 h 2 ⊥,2 ∂ s+1 ∂ ⊥,1 ∂ ⊥,2 v 2 L 2 (Kce) ≤ τ 2s+4 σ d 2s+2 c d 2 e ∂ s+1 ∂ ⊥,1 ∂ ⊥,2 v 2 L 2 (Kce) ≤ τ 2s+4 σ d 2 γ c d 2γe−2 e r s+3−γc c ρ 2−γe ce ∂ s+1 ∂ ⊥,1 ∂ ⊥,2 v 2 L 2 (Kce) , where d 2 γ c d 2γe−2 e ≤ d 2γc−2 c . Furthermore, if α1 = 0, h 2s+2 ∂ s+1 ∂ ⊥,1 v 2 L 2 (Kce) ≤ τ 2s+2 σ d 2s+2 c ∂ s+1 ∂ ⊥,1 v 2 L 2 (Kce) ≤ τ 2s+2 σ d 2γc−2 c r s+2−γc c ∂ s+1 ∂ ⊥,1 v 2 L 2 (Kce) .
Therefore,
(I)ce ≤ 1 − σ σ 2s+4 d 2γc−2 c α 1 ≤1 r s+2+α 1 −γc c ρ (1+α 1 −γe) + ce ∂ s+1 ∂ ⊥,1 ∂ α 1 ⊥,2 v 2 L 2 (Kce) .
The estimates for (I)c,e,0 follow from the same technique:
(I)e ≤ α 1 ≤1 τ 2s+4 σ d 2s+2 c r (1+α 1 −γe) + e ∂ s+1 ∂ ⊥,1 ∂ α 1 ⊥,2 v 2 L 2 (Ke) , (I)c ≤ α 1 ≤1 τ 2s+4 σ d 2γc−2 c r s+2+α 1 −γc c ∂ s+1 ∂ ⊥,1 ∂ α 1 ⊥,2 v 2 L 2 (Kc) , (I)0 ≤ α 1 ≤1 τ 2s+4 σ d 2s+2 c ∂ s+1 ∂ ⊥,1 ∂ α 1 ⊥,2 v 2 L 2 (K 0 ) .
Hence, from (2.1), there exists C > 0 dependent only on Cv and σ and A > 0 dependent only on Av and σ such that
(I) ≤ CA 2s+6 ((s + 3)!) 2 d 2γc−2 c . (A.31)
We then apply the same argument to the terms (II) and (III). Indeed, if s + 1 + α1 + α2 − γc ≥ 0
(II)ce = α 1 ,α 2 ≤1 h 2α 1 h 2s ⊥,1 h 2α 2 ⊥,2 ∂ α 1 ∂ s+1 ⊥,1 ∂ α 2 ⊥,2 v 2 L 2 (Kce) ≤ τ 2s+4 σ α 1 ,α 2 ≤1 d 2α 1 c d 2s+2α 2 e ∂ α 1 ∂ s+1 ⊥,1 ∂ α 2 ⊥,2 v 2 L 2 (Kce) ≤ τ 2s+4 σ α 1 ≤1 d 2 γ c d 2γe−2 e r s+1+α 1 +α 2 −γc c ρ s+1+α 2 −γe ce ∂ α 1 ∂ s+1 ⊥,1 ∂ α 2 ⊥,2 v 2 L 2 (Kce) ≤ τ 2s+4 σ α 1 ≤1 d 2γc−2 c r s+1+α 1 +α 2 −γc c ρ s+1+α 2 −γe ce ∂ α 1 ∂ s+1 ⊥,1 ∂ α 2 ⊥,2 v 2 L 2 (Kce) ,
where at the last step we have used that γe > 1 and de ≤ dc. If s + 1 + α1 + α2 − γc  0 dependent only on σ, Cappx2, Cv and A > 0 dependent only on σ, Av such that
∂ ⊥,1 (v − Πpv) 2 L 2 ( K) ≤ C h ⊥,1 h h ⊥,2 Ψp,sA 2s+6 d 2(γc−1) c + d 2(γe−1) c ((s + 3)!) 2 .
Considering that
∂ ⊥,1 (v − Πpv) 2 L 2 (K) ≤ h h ⊥,2 h ⊥,1 ∂ ⊥,1 (v − Πpv) 2 L 2 ( K)
and considering that the estimate for the other term at the left-hand side of (A.29) is obtained by exchanging {h, ∂} ⊥,1 with {h, ∂} ⊥,2 completes the proof.
Lemma A.11. Let d = 3, ∈ N and K = J k 1 × J k 2 × J k 3 for 0  0 dependent only on σ, Cappx1, Cv and A > 0 dependent only on σ, Av such that for all p ∈ N and all 1 ≤ s ≤ p
v − Π K p v 2 L 2 (K) ≤ CΨp,sA 2s+6 d 2(γc−1) c + d 2(γe−1) c ((s + 3)!) 2 . (A.33)
Proof. The proof follows closely that of Lemmas A.9 and A.10; we use the same notation. From Lemma A.4 and rescaling, we have
v − Πpv 2 L 2 ( K) ≤ Cappx1Ψp,s 1 h h ⊥,1 h ⊥,2   α 1 ,α 2 ≤1 h 2s+2 h 2α 1 ⊥,1 h 2α 2 ⊥,2 ∂ s+1 ∂ α 1 ⊥,1 ∂ α 2 ⊥,2 v 2 L 2 (K) + α 1 ,α 2 ≤1 h 2α 1 h 2s+2 ⊥,1 h 2α 2 ⊥,2 ∂ α 1 ∂ s+1 ⊥,1 ∂ α 2 ⊥,2 v 2 L 2 (K) + α 1 ,α 2 ≤1 h 2α 1 h 2α 2 ⊥,1 h 2s+2 ⊥,2 ∂ α 1 ∂ α 2 ⊥,1 ∂ s+1 ⊥,2 v 2 L 2 (K)   . (A.34)
Most terms at the right-hand side above have already been considered in the proofs of Lemmas A.9 and A.10, and the terms with α1 = α2 = 0 can be estimated similarly; the observation that
v − Πpv 2 L 2 (K) ≤ h h ⊥,1 h ⊥,2 v − Πpv 2 L 2 ( K)
concludes the proof.
We summarize Lemmas A.9 to A.11 in the following result.
Lemma A.12. Let d = 3, ∈ N and K = J k 1 × J k 2 × J k 3 such that 0  0 dependent only on σ, Cappx1, Cappx2, Cv and A > 0 dependent only on σ, Av such that for all p ∈ N and all 1 ≤ s ≤ p
v − Π K p v 2 H 1 (K) ≤ CΨp,sA 2s+6 d 2(γc−1) c + d 2(γe−1) c ((s + 3)!) 2 . (A.35)
We then consider elements on the faces (but not abutting edges) of Q.
Lemma A.13. Let d = 3, ∈ N and K = J k 1 × J k 2 × J k 3 such that kj = 0 for one j ∈ {1, 2, 3} and 0  0 dependent only on σ, Cappx1, Cappx2, Cv and A > 0 dependent only on σ, Av such that
v − Π K p 1 p 2 p 3 v 2 H 1 (K) C Ψp,sA 2s+6 (d K c ) 2(min(γc,γe)−1) ((s + 3)!) 2 + (d K e ) 2(min(γc,γe)−2) σ 2 A 8 . (A.36)
Proof. We write da = d K a , a ∈ {c, e}. Suppose, for ease of notation, that j = 3, i.e. k3 = 0. The projector is then given by Π K pp1 = π k 1 p ⊗ π k 2 p ⊗ π 0 1 . Also, we denote h ⊥,2 = σ and ∂ ⊥,2 = ∂x 3 . By (A. 16),
∂ (v − Π K pp1 v) 2 L 2 (K) ≤ Cappx2   Ψp,s α 1 ,α 2 ≤1 h 2s h 2α 1 ⊥,1 h 2α 2 ⊥,2 ∂ s+1 ∂ α 1 ⊥,1 ∂ α 2 ⊥,2 v 2 L 2 (K) + α 1 ≤1 h 2s+2 ⊥,1 h 2α 1 ⊥,2 ∂ ∂ s+1 ⊥,1 ∂ α 1 ⊥,2 v 2 L 2 (K) + α 1 ≤1 h 2α 1 ⊥,1 h 4 ⊥,2 ∂ ∂ α 1 ⊥,1 ∂ 2 ⊥,2 v 2 L 2 (K)   = Cappx2 (I) + (II) + (III) .
The bounds on the terms (I) and (II) can be derived as in Lemma A.9, and give
(I) + (II) ≤ CΨp,sA 2s+6 (d K c ) 2 + (d K c ) 2(γc−1) ((s + 3)!) 2 .
We consider then term (III): with the usual notation, writing γ = γc − γe,
(III)ce = α 1 ≤1 h 2α 1 ⊥,1 h 4 ⊥,2 ∂ ∂ α 1 ⊥,1 ∂ 2 ⊥,2 v 2 L 2 (Kce) ≤ α 1 ≤1 τ 4+2α 1 σ d 2 γ−2 c d 2γe−4 e σ 4 r 3+α 1 −γc c ρ 2+α 1 −γe ce ∂ ∂ α 1 ⊥,1 ∂ 2 ⊥,2 v 2 L 2 (Kce) ≤ Cτ 6 σ d 2 γ−2 c d 2γe−4 e σ 4 A 8 .
(A.37)
Note that dc ≥ de and
d γ c d γe e ≤ 1 γ d γe e if γ ≥ 0 d γ e d γe e if γ ≥ 0 ≤ d min(γc,γe) e , (A.38)
where we have also used that dc ≤ 1. Hence,
(III)ce ≤ Cτ 6 σ d 2 min(γe,γc)−6 e σ 4 A 8 ≤ Cτ 6 σ d 2 min(γe,γc)−4 e σ 2 A 8 . (A.39)
The bounds on the terms (III)c,e,0 follow by the same argument:
(III)e ≤ Cτ 6 σ d 2γe−4 e σ 4 A 8 , (III)c ≤ Cτ 6 σ d 2γc−6 c σ 4 A 8 ≤ Cτ 6 σ d 2γc−4 e σ 2 A 8 , (III)0 ≤ Cτ 6 σ σ 4 A 8 .
Then,
∂ ⊥,1 (v − Π K pp1 v) 2 L 2 (K) ≤ Cappx2   (p − s)! (p + s)! α 1 ≤1 h 2s+2 h 2α 1 ⊥,2 ∂ s+1 ∂ ⊥,1 ∂ α 1 ⊥,2 v 2 L 2 (K) + α 1 ,α 2 ≤1 h 2α 1 h 2s ⊥,1 h 2α 2 ⊥,2 ∂ α 1 ∂ s+1 ⊥,1 ∂ α 2 ⊥,2 v 2 L 2 (K) + α 1 ≤1 h 2α 1 h 4 ⊥,2 ∂ α 1 ∂ ⊥,1 ∂ 2 ⊥,2 v 2 L 2 (K)   ≤ Cappx2 (I) + (II) + (III) .
The bounds on the first two terms at the right-hand side above can be obtained as in Lemma A.10:
(I) + (II) ≤ CΨp,sA 2s+6 (d K c ) 2(γc−1) + (d K c ) 2(γe−1) ((s + 3)!) 2 ,
while the last term can be bounded as in (A.39),
(III)ce ≤ τ 6 σ d 2 γ c d 2γe−6 e σ 4 A 8 ≤ Cτ 6 σ d 2 min(γc,γe)−4 e σ 2 A 8 , (III)e ≤ τ 6 σ d 2γe−6 e σ 4 A 8 ≤ Cτ 6 σ d 2γe−4 e σ 2 A 8 , (III)c ≤ τ 6 σ d 2γc−6 c σ 4 A 8 ≤ Cτ 6 σ d 2γc−4 e σ 2 A 8 , (III)0 ≤ τ 6 σ σ 4 A 8 , so that α 1 ≤1 h 2α 1 h 4 ⊥,2 ∂ α 1 ∂ ⊥,1 ∂ 2 ⊥,2 v 2 L 2 (K) ≤ Cd 2 min(γc,γe)−4 e σ 2 A 8 .
The same holds true for the last term of the gradient of the approximation error, given by
∂ ⊥,2 (v − Π K pp1 v) 2 L 2 (K) ≤ Cappx2   Ψp,s α 1 ≤1 h 2s+2 h 2α 1 ⊥,1 ∂ s+1 ∂ α 1 ⊥,1 ∂ ⊥,2 v 2 L 2 (K) + α 1 ≤1 h 2α 1 h 2s+2 ⊥,1 ∂ α 1 ∂ s+1 ⊥,1 ∂ ⊥,2 v 2 L 2 (K) + α 1 ,α 2 ≤1 h 2α 1 h 2α 2 ⊥,1 h 2 ⊥,2 ∂ α 1 ∂ α 2 ⊥,1 ∂ 2 ⊥,2 v 2 L 2 (K)   ≤ Cappx2 (I) + (II) + (III) .
From Lemma A.10, we obtain
(I) + (II) ≤ CΨp,sA 2s+6 (d K c ) 2(γc−1) + (d K c ) 2(γe−1) ((s + 3)!) 2 ,
whereas for the third term, it holds that if α1 + α2
+ 2 − γc ≥ 0 (III)ce ≤ τ 6 σ d 2 γ c d 2γe−4 e σ 2 A 8 ≤ Cτ 6 σ d 2 min(γc,γe)−4 e σ 2 A 8 , (III)c ≤ τ 6 σ d 2γc−4 c σ 2 A 8 ,
and if α1 + α2 + 2 − γc  3/2, γe > 1. There exists a constant C0 > 0 such that if p ≥ C0 , there exist constants C, b > 0 such that for every ∈ N holds
K:d K e >0 v − Π ,p hp,d v 2 H 1 (K) ≤ Ce −b .
Proof. We suppose, without loss of generality, that γc ∈ (3/2, 5/2), and γe ∈ (1, 2). The general case follows from the inclusion J γ 1 (Q; C, E) ⊂ J γ 2 (Q; C, E), valid for γ1 ≥ γ2. Fix any C0 > 0 and choose p ≥ C0 . For all A > 0 there exist C1, b1 > 0 such that (see, e.g., [45,Lemma 5.9])
∀p ∈ N : min 1≤s≤p Ψp,sA 2s (s!) 2 ≤ C1e −b 1 p .
From (A.35) and (A.36), there holds
K:d K e >0 v − Π ,p hp,d v 2 H 1 (K) ≤ C2    K:d K e >0 e −b 1 (d K c ) 2(min(γc,γe)−1) + K:d K e >0,d K f =0 (d K e ) 2(min(γe,γc)−2) σ 2    = C2 (I) + (II) ,
where d K f indicates the distance of an element K to one of the faces of Q. There holds directly (I) ≤ C 2 e −b 1 . Furthermore, because (min(γc, γe) − 2)  1. There exists a constant C0 > 0 such that if p ≥ C0 , there exist constants C, b > 0 such that
K:d K c >0 v − Π ,p hp,d v 2 H 1 (K) ≤ Ce −b , ∀ ∈ N.
A.6 Estimates on elements along an edge in three dimensions
In the following lemma, we consider the elements K along one edge, but separated from the singular corner.
Lemma A.16. Let d = 3, e ∈ E and let K ∈ G 3 be such that d K c > 0 for all c ∈ C and d K e = 0. Let Cv, Av > 0. Then, if v ∈ J γ (Q; C, E; Cv, Av) with γc ∈ (3/2, 5/2), γe ∈ (1, 2), there exist C, A > 0 such that for all p ∈ N and all 1 ≤ s ≤ p, with (p1, p2, p3) ∈ N 3 such that p = p, p ⊥,
1 = 1 = p ⊥,2 , v − Π K p 1 p 2 p 3 v 2 H 1 (K) ≤ C σ 2{min{γc−1,s}( −k) Ψp,sA 2s ((s + 3)!) 2 + σ 2(min(γe,γc)−1) , (A.40) where k ∈ {1, . . . , } is such that d K c = σ −k+1 .
Proof. We suppose that K = J k × J 0 × J 0 for some k ∈ {1, . . . , }, the elements along other edges follow by symmetry. This implies that the singular edge is parallel to the first coordinate direction. Furthermore, we denote Π K p11 = π k p ⊗ (π 0 1 ⊗ π 0 1 ) = π ⊗ π ⊥ . For α = (α1, α2, α3) ∈ N 3 0 , we write α = (α1, 0, 0) and α ⊥ = (0, α2, α3). Also,
h = |J k | = σ −k (1 − σ) h ⊥ = σ . We have v − Π K p11 v = v − π ⊥ v + π ⊥ v − π v . (A.41)
We start by considering the first terms at the right-hand side of the above equation. We also compute the norms over Kce = K ∩ Qce; the estimate on the norms over Kc = K ∩ Qc and Ke = K ∩ Qe follow by similar or simpler arguments. By (A.21) from Lemma A.8, we have that if γc  0 depending only on σ and Av such that 
|α ⊥ |≤1 ∂ α ⊥ π ⊥ (v − π v) 2 L 2 (Kce) τ 2s+2 σ Ψp,sh 2 min{γc−1,s} h 2   |α ⊥ |≤1 h −2|α ⊥ | r (s+1+|α ⊥ |−γc) + c ∂ α v 2 L 2 (Kce) + |α ⊥ |=2 h 2(γe−1) ⊥ h −2γe r s+3−γc c ρ 2−γe ce ∂ α v 2 L 2 (Kce)   σ 2( −k) min{γc−1,s} Ψp,sA 2s+4 ((s + 3)!) 2 , (A.45a) and similarly |α ⊥ |≤1 ∂ α ⊥ π ⊥ (v − π v) 2 L 2 (Ke) σ 2( −k) min{γc,∂ α π ⊥ (v − π v) 2 L 2 (K) |α |≤1   |α ⊥ |≤1 h 2|α ⊥ | ⊥ ∂ α ⊥ ∂ α (v − π v) 2 L 2 (K) + |α ⊥ |=2 h 2γe ⊥ r 2−γe e ∂ α ⊥ ∂ α (v − π v) 2 L 2 (K)   τ 2s+2 σ Ψp,sh 2 min{γc−1,s}   |α |=s+1 |α ⊥ |≤1 h 2|α ⊥ | ⊥ r (s+1−γc) + c ∂ α ∂ α ⊥ v 2 L 2 (K) + |α |=s+1 |α ⊥ |=2 h 2γe ⊥ r 2−γe e r (s+1−γc) + c ∂ α ∂ α ⊥ v 2 L 2 (K)   .
As before, there exists A > 0 depending only on σ and Av such that Lemma A.17. Let d = 3 and v ∈ J γ (Q; C, E) with γc > 3/2, γe > 1. There exists a constant C0 > 0 such that if p ≥ C0 , there exist constants C, b > 0 such that
|α |≤1 ∂ α π ⊥ (v − π v) 2 L 2 (Kce) τ 2s+2 σ Ψp,sh 2 min{γc−1,s}   |α |=s+1 |α ⊥ |≤1 h 2|α ⊥ | ⊥ h −2|α ⊥ | r (s+1+|α ⊥ |−γc) + c ∂ α v 2 L 2 (Kce) + |α |=s+1 |α ⊥ |=2 h 2γe ⊥ h −2γe r s+3−γc c ρ 2−γe ce ∂ α v 2 L 2 (Kce)   σ 2( −k) min{γc−1,s} Ψp,sA 2s+4 ((s + 3)!) 2 , (A.46a) and |α |≤1 ∂ α π ⊥ (v − π v) 2 L 2 (Ke) τ 2s+2 σ Ψp,sh 2 min{γc−1,s}   |α |=s+1 |α ⊥ |≤1 h 2|α ⊥ | ⊥ r (s+1−γc) + c ∂ α v 2 L 2 (Ke) + |α |=s+1 |α ⊥ |=2 h 2γe ⊥ r (s+1−γc) + c r 2−γe e ∂ α v 2 L 2 (Ke)   σ 2( −k) min{γc−1,K:d K c >0, d K e =0 v − Π ,p hp,d v H 1 (K) ≤ Ce −b , ∀ ∈ N.
Proof. As in the proof of Lemma A.14, we may assume that γc ∈ (3/2, 5/2) and γe ∈ (1, 2). The proof of this statements follows by summing over the right-hand side of (A.40), i.e.,
K:d K c >0, d K e =0 v − Π ,p hp,d v 2 H 1 (K) ≤ C k=1 σ 2 min{γc−1,s}( −k) Ψp,sA 2s ((s + 3)!) 2 + σ 2(min(γc,γe)−1)
= C((I) + (II)).
We have (II) σ 2(min(γc,γe)−1) ; the observation that for all A > 0 there exist C1, b1 > 0 such that
min 1≤s≤p Ψp,s((s + 3)!) 2 A 2s ≤ C1e −b 1 p ,
(see, e.g., [45,Lemma 5.9]). Combining with p ≥ C0 concludes the proof.
A.7 Estimates at the corner
The lemma below follows from classic low-order finite element approximation results and from the embedding J 2 γ (Q; C, E) ⊂ H 1+θ (Q), valid for a θ > 0 if γc − d/2 > 0, for all c ∈ C, and, when d = 3, γe > 1 for all e ∈ E (see, e.g., [43,Remark 2.3]). 
Lemma A.18. Let d ∈ {2, 3}, K = × d i=1 J 0 . Then, if v ∈ J γ (Q; C, E) with γc > 1, for all c ∈ C, if d = 2,
A.8 Exponential convergence
The exponential convergence of the approximation in the full domain Q follows then from Lemmas A.14, A.15, A.17, and A.18.
Proposition A.19. Let d ∈ {2, 3}, v ∈ J γ (Q; C, E) with γc > 1, for all c ∈ C, if d = 2,
γc > 3/2 and γe > 1, for all c ∈ C and e ∈ E, if d = 3.
Then, there exist constants cp > 0 and C, b > 0 such that, for all ∈ N,
v − Π ,cp hp,d v H 1 (Q) ≤ Ce −b .
With respect to the dimension of the discrete space
N dof = dim(X ,cp hp,d ), the above bound reads v − Π ,cp hp,d v H 1 (Q) ≤ C exp(−bN 1/(2d) dof
).
A.9 Explicit representation of the approximant in terms of continuous basis functions
Let p ∈ N. Letζ1(x) = (1 + x)/2 andζ2 = (1 − x)/2. Let alsoζn(x) = 1 2
x −1 Ln−2(ξ)dξ, for n = 3, . . . , p + 1, where Ln−2 denotes the L ∞ ((−1, 1))-normalized Legendre polynomial of degree n − 2 introduced in Section A.2. Then, fix ∈ N and write ζ k n =ζn • φ k , n = 1, . . . , p + 1 and k = 0, . . . , , with the affine map φ k : J k → (−1, 1) introduced in Section A.2. We construct those functions explicitly: denoting J k = (x k , x k+1 ) and h k = |x k+1 − x k |, there holds, for x ∈ J k ,
ζ k 1 (x) = 1 h k (x − x k ), ζ k 2 (x) = 1 h k (x k+1 − x), (A.47) and ζ k n (x) = 1 h k x x k Ln−2(φ k (η))dη n = 3, . . . , p + 1. (A.48)
Then, for any element K ∈ G 3 , with
K = J k 1 × J k 2 × J k 3 , there exist coefficients c K i 1 ...i d such that Π ,p hp,d u | K (x1, x2, x3) = p+1 i 1 ,i 2 ,i 3 =1 c K i 1 ...i d ζ k 1 i 1 (x1)ζ k 2 i 2 (x2)ζ k 3 i 3 (x3), ∀(x1, x2, x3) ∈ K (A.49)
by construction. We remark that, whenever ij > 2 for all j = 1, 2, 3, the basis functions vanish on the boundary of the element: 
ζ k 1 i 1 ζ k 2 i 2 ζ k 3 i 3 | ∂K = 0 if ij ≥ 3, j = 1, 2, 3. Furthermore, write ψ K i 1 ...i d (x1, x2, x3) = ζ k 1 i 1 (x1)ζ k 2 i 2 (x2)ζ k 3 i 3 (
A.9.1 Explicit bounds on the coefficients
We derive here a bound on the coefficients of the local projectors with respect to the norms of the projected function. We will use that
Li • φ k L 2 (J k ) = h k 2 1/2
Li L 2 ((−1,1)) = h k 2i + 1 
|c K i 1 ...i d | ≤ C u W 1,1 mix (Q)              d j=1 ij if ti 1 ...i d = 0, ( + 1) d j 1 =1 d j 2 =j 1 +1 ij 1 ij 2 if ti 1 ...i d = 1, ( + 1) 2 d j=1 ij if ti 1 ...i d = 2, ( + 1) d if ti 1 ...i d = 3. (A.52) Proof. Let d = 3 and K = J k 1 × J k 2 × J k 3 ∈ G 3 .
Internal modes. We start by considering the case of the coefficients of internal modes, i.e., c K i 1 ,i 2 ,i 3 as defined in (A.49) for in ≥ 3, n = 1, 2, 3. Let then i1, i2, i3 ∈ {3, . . . , p + 1} and write L k n = Ln • φ k : there holds
c K i 1 ,i 2 ,i 3 = (2i1 − 3)(2i2 − 3)(2i3 − 3) K (∂x 1 ∂x 2 ∂x 3 u(x1, x2, x3)) L k 1 i 1 −2 (x1)L k 2 i 2 −2 (x2)L k 3 i 3 −2 (x3)dx1dx2dx3. (A.53)
If u ∈ W 1,1 mix (K), since Ln L ∞ (−1,1) = 1 for all n, we have
|c K i 1 ...i d | ≤ (2i1 − 3)(2i2 − 3)(2i3 − 3) ∂x 1 ∂x 2 ∂x 3 u L 1 (K) in ≥ 3, n = 1, 2, 3, (A.54) hence, K∈G 3 |c K i 1 ...i d | ≤ (2i1 − 3)(2i2 − 3)(2i3 − 3) ∂x 1 ∂x 2 ∂x 3 u L 1 (Q) in ≥ 3, n = 1, 2, 3. (A.55)
Face modes. We continue with face modes and fix, for ease of notation, i1 = 1. We also denote F = J k 2 × J k 3 . The estimates will then also hold for i1 = 2 and for any permutation of the indices by symmetry. We introduce the trace inequality constant C T,1 , independent of K, such that, for all v ∈ W 1,1 (Q) andx ∈ (0, 1),
v(x, ·, ·) L 1 (F ) ≤ v(x, ·, ·) L 1 ((0,1) 2 ) ≤ C T,1 v L 1 (Q) + ∂x 1 v L 1 (Q) . (A.56)
This follows from the trace estimate in [44,Lemma 4.2] and from the fact that
v(x, ·, ·) L 1 ((0,1) 2 ) ≤ C min 1 |1 −x| v L 1 ((x,1)×(0,1) 2 ) + ∂x 1 v L 1 ((x,1)×(0,1) 2 ) , 1 |x| v L 1 ((0,x)×(0,1) 2 ) + ∂x 1 v L 1 ((0,x)×(0,1) 2 ) .
There holds, for i2, i3 ∈ {3, . . . , p + 1},
c K 1,i 2 ,i 3 = (2i2 − 3)(2i3 − 3) F ∂x 2 ∂x 3 u(x k 1 , x2, x3) L k 2 i 2 −2 (x2)L k 3 i 3 −2 (x3)dx2dx3. (A.57)
Since the Legendre polynomials are L ∞ normalized and using the trace inequality (A.56),
|c K 1,i 2 ,i 3 | ≤ (2i2 − 3)(2i3 − 3) (∂x 2 ∂x 3 u)(x k 1 , ·, ·) L 1 (F ) ≤ C T,1 (2i2 − 3)(2i3 − 3) u W 1,1 mix (Q)
. (A.58) Summing over all internal faces, furthermore,
K∈G 3 |c K 1,i 2 ,i 3 | ≤ (2i2 − 3)(2i3 − 3) k 1 =0 (∂x 2 ∂x 3 u)(x k 1 , ·, ·) L 1 ((0,1) 2 ) ≤ C T,1 ( + 1)(2i2 − 3)(2i3 − 3) u W 1,1 mix (Q) .
(A.59)
Edge modes. We now consider edge modes. Fix for ease of notation i1 = i2 = 1; as before, the estimates will hold for (i1, i2) ∈ {1, 2} 2 and for any permutation of the indices. By the same arguments as for (A.56), there exists a trace constant C T,2 such that, denoting e = J k 3 , for all v ∈ W 1,1 ((0, 1) 2 ) and for allx ∈ (0, 1),
v(x, ·) L 1 (e) ≤ v(x, ·) L 1 ((0,1)) ≤ C T,2 u L 1 ((0,1) 2 ) + ∂x 2 u L 1 ((0,1) 2 ) . (A.60) By definition,
c K 1,1,i 3 = (2i3 − 3) e ∂x 3 u(x k 1 , x k 2 , x3) L k 3 i 3 −2 (x3)dx3. (A.61)
Using (A.56) and (A.60)
|c K 1,1,i 3 | ≤ (2i3 − 3) (∂x 3 u)(x k 1 , x k 2 , ·) L 1 (e) ≤ C T,1 C T,2 (2i3 − 3) u W 1,1 mix (Q) . (A.62)
Summing over edges, in addition,
K∈G 3 |c K 1,1,i 3 | ≤ (2i3 − 3) k 1 =0 k 2 =0 (∂x 3 u)(x k 1 , x k 2 , ·) L 1 ((0,1)) ≤ C T,1 C T,2 ( + 1) 2 (2i3 − 3) u W 1,1 mix (Q) .
(A.63)
Node modes. Finally, we consider the coefficients of nodal modes, i.e., c K i 1 ,i 2 ,i 3 for i1, i2, i3 ∈ {1, 2}, which by construction equal function values of u, e.g.
c111 = u(x k 1 , x k 2 , x k 3 ).
(A.64)
The Sobolev imbedding W 1,1 mix (Q) → L ∞ (Q) and scaling implies the existence of a uniform constant C imb such that, for any v ∈ W 1,1
mix (Q) v L ∞ (K) ≤ v L ∞ (Q) ≤ C imb v W 1,1 mix (Q)
. Then, by construction,
|c K i 1 ,i 2 ,i 3 | ≤ u L ∞ (K) ≤ C imb u W 1,1 mix (Q) ∀i1, i2, i3 ∈ {1, 2}. (A.65)
Summing over nodes, it follows directly that There exists a constant C > 0 such that, for all u ∈ J d γ (Q; C, E),
K∈G 3 |c K i 1 ,i 2 ,i 3 | ≤ K∈G 3 u L ∞ (K) ≤ C imb ( + 1) 3 u W 1,1 mix (Q) ∀i1, i2, i3 ∈ {1,u W 1,1 mix (Q) ≤ C u J d γ (Q) .
Proof. We recall the decomposition of Q as
Q = Q0 ∪ QC ∪ QE ∪ QCE , where QE = QCE = ∅ if d = 2. There holds u W 1,1 mix (Q 0 ) ≤ C|Q0| 1/2 u H d (Q 0 ) ≤ C|Q0| 1/2 u J d γ (Q) . (A.67)
We now consider the subdomain Qc, for any c ∈ C. There holds, with constant C that depends only on γc and on |Qc|,
u W 1,1 mix (Qc) = u W 1,1 (Qc) + 2≤|α|≤d |α|∞≤1 ∂ α u L 1 (Qc) ≤ C|Qc| 1/2 u H 1 (Qc) + C 2≤|α|≤d |α|∞≤1 r −(|α|−γc) + c L 2 (Qc) r (|α|−γc) + c ∂ α u L 2 (Qc) ≤ C u J d γ (Q) , (A.68)
where the last inequality follows from the fact that γc > d/2, hence the norm r −(|α|−γc) + c L 2 (Qc) is bounded for all |α| ≤ d. Consider then d = 3 and any e ∈ E. Suppose also, without loss of generality, that γc − γe > 1/2 and γe  1, then r −|α ⊥ |+γe e L 2 (Qe) is bounded by a constant depending only on γe and |Qe| as long as α is such that |α ⊥ | ≤ 2. Hence, denoting by ∂ the derivative in the direction parallel to e,
u W 1,1 mix (Qe) = u W 1,1 (Qe) + |α ⊥ |=1 ∂ ∂ α ⊥ u L 1 (Qe) + α 1 =0,1 ∂ α 1 ∂ ⊥,1 ∂ ⊥,2 u L 1 (Qe) ≤ C|Qe| 1/2   u H 1 (Qe) + |α ⊥ |=1 ∂ ∂ α ⊥ v L 2 (Qe)   + C α 1 =0,1 r −2+γe e L 2 (Qe) r 2−γe e ∂ α 1 ∂ ⊥,1 ∂ ⊥,2 u L 2 (Qe) ≤ C u J 3 γ (Q) . (A.69)
Since x ≤ rc(x) ≤ ε for all x ∈ Qce, and because Qce ⊂ x ∈ (0, ε), (x ⊥,1 , x ⊥,2 ) ∈ (0, ε 2 ) 2 , there holds r
−(γe+1−γc) + c r −2+γe e L 2 (Qce) ≤ x −(γe+1−γc) + L 2 ((0, ε)) r −2+γe e L 2 ((0, ε 2 ) 2 ) ≤ C,
for a constant C that depends only on ε, γc, and γe. Hence,
u W 1,1 mix (Qce) = u W 1,1 (Qce) + |α ⊥ |=1 ∂ ∂ α ⊥ u L 1 (Qce) + α 1 =0,1 ∂ α 1 ∂ ⊥,1 ∂ ⊥,2 u L 1 (Qce) ≤ C|Qce| 1/2 u H 1 (Qe) + C |α ⊥ |=1 r −(2−γc) + c L 2 (Qce) r (2−γc) + c ∂ ∂ α ⊥ u L 2 (Qce) + C α 1 =0,1 r −(α 1 +γe−γc) + c r −2+γe e L 2 (Qce) r (α 1 +2−γc) + c ρ 2−γe ce ∂ α 1 ∂ ⊥,1 ∂ ⊥,2 u L 2 (Qce) ≤ C u J 3 γ (Q) , (A.70) with C independent of u.
Combining inequalities (A.67) to (A.70) concludes the proof.
The following statement is a direct consequence of the two lemmas above and the fact that 
ψ K i 1 ...i d L ∞ (K) ≤ 1hp,d ) v − v hp H 1 (Q) ≤ C hp e −b hp .
For ε > 0, we choose Then, by (A.73),
L := 1 b hp |log(ε/C hp )| , (A.73) so that v − v L hp H 1 (Q) ≤ ε. Furthermore, v L hp = N 1d i 1 ,σ −L ≤ σ − 1 b hp log(C hp ) ε − 1 b hp log(1/σ) .
This concludes the proof of Items 1 and 2. Finally, Item 3 follows from Proposition A.24 and the fact that p ≤ Cp (1 + |log(ε)|) for a constant Cp > 0 independent of ε.
A.10 Combination of multiple patches
The approximation results in the domain Q = (0, 1) d can be generalized to include the combination of multiple patches. We give here an example, relevant for the PDEs considered in Section 5. For the sake of conciseness, we show a single construction that takes into account all singularities of the problems in Section 5. We will then use this construction to prove expression rate bounds for realizations of NNs. 
ψ1,±(x) = ± a 2 x, ψ2,±(x) = ± a − a 2 x .
For all ∈ N, define then
G 1 = i∈{1,2}, ∈{+.−} ψi, (G 1 ). Consequently, for d = 2, 3, denote G d = { × d i=1
Ki : K1, . . . , K d ∈ G 1 }, see Figure 3. The hp space in Ω = (−a, a) d is then given by
X ,p hp,d = {v ∈ H 1 (Ω) : v | K ∈ Qp(K), for all K ∈ G d }.
Finally, recall the definition of π ,p hp from (A.12) and construct π ,p hp : W 1,1 ((−a, a)) → X ,p hp,1 such that, for all v ∈ W 1,1 ((−a, a)),  where for all j = 1, . . . , d and ij = 1, . . . , N 1d ,ṽi j ∈ X ,p hp,1 with support in at most two, neighboring elements of G 1 . Finally, there exist constants C1, C2 > 0 independent of such that Proof. The statement is a direct consequence of Propositions A.19 and A.24. We start the proof by showing that for any function v ∈ W 1,1 mix (Ω), the approximation Π ,p hp,d v is continuous; the rest of the theorem will then follow from the results in each sub-patch. Let now w ∈ W 1,1 ((−a, a)). Then, it holds that π ,p hp w |I ∈ C(I), for all I ∈ {(0, a/2), (a/2, a), (−a/2, 0), (−a, −a/2)}, by definition (A.76). Furthermore, by the nodal exactness of the local projectors, forx ∈ {−a/2, 0, a/2}, there holds
π ,p hp v | (0, a 2 ) = π ,p hp (v| (0, a 2 ) • ψ1,+) • ψ −1 1,+ , π ,p hp v | ( a 2 ,a) = π ,p hp (v| ( a 2 ,a) • ψ2,+) • ψ −1 2,+ , π ,p hp v | (− a 2 ,0) = π ,p hp (v| (− a 2 ,0) • ψ1,− • ψ −1 1,− , π ,p hp v | (−a,− a 2 ) = π ,p hp (v| (−a,− a 2 ) • ψ2,−) • ψ −1 2,− . (A.Π ,p hp,d = d i=1 π ,p hp .lim x→x − ( π ,p hp w)(x) = w(x) = lim x→x + ( π ,p hp w)(x),
implying then that π ,p hp w is continuous. Since Π ,p hp,d = d j=1 π ,p hp , this implies that Π ,p hp,d v is continuous for all v ∈ W 1,1 mix (Ω). Fix k ∈ {1, . . . , 4 d } such that v ∈ J γ k (Ω k ; C k , E k ). There exist then, by Proposition A.19, constants C, b, cp > 0 such that for all ∈ N v − Π Proof of Lemma 5.5. For any two sets X, Y ⊂ Ω, we denote by distΩ(X, Y ) the infimum of Euclidean lengths of paths in Ω connecting an element of X with one of Y . We introduce several domain-dependent quantities to be used in the construction of the triangulation T with the properties stated in the lemma.
Let E denote the set of edges of the polygon Ω. For each corner c ∈ C at which the interior angle of Ω is smaller than π (below called convex corner), we fix a parallelogram Gc ⊂ Ω and a bijective, affine transformation Fc : (0, 1) 2 → Gc such that • Fc((0, 0)) = c,
• two edges of Gc coincide partially with the edges of Ω abutting at the corner c.
If at c the interior angle of Ω is greater than or equal to π (both are referred to by slight abuse of terminology as nonconvex corner), the same properties hold, with Fc : (−1, 1) × (0, 1) → Gc if the interior angle equals π, and Fc : (−1, 1) 2 \ (−1, 0] 2 → Gc else, and with Gc having the corresponding shape. Let now dc,1 := sup{r > 0 : Br(c) ∩ Ω ⊂ Gc}, dC,1 := min c∈C dc,1.
Then, for each c ∈ C, let e1 and e2 be the edges abutting c, and define dc,2 := distΩ e1 ∩ B √ Then, in case Ω is a triangle, let d0 be half of the radius of the inscribed circle, else let d0 := 1 3 dE  0.
For any shape regular triangulation T of R 2 , such that for all K ∈ T , K ∩ ∂Ω = ∅, denote TΩ = {K ∈ T : K ⊂ Ω} and h(TΩ) = maxK∈T Ω h(K), where h(K) denotes the diameter of K. Denote by NΩ the set of nodes of T that are in Ω.  Let T be a triangulation of R 2 such that
h(TΩ) ≤ min d0 √ 2 , dC,1 √ 2 + 1 , dC,2 2 √ 2 , dE 2 √ 2 , (B.1)
and such that for all K ∈ T it holds K ∩ ∂Ω = ∅.
The hat-function basis {φn}n∈N Ω is a basis for P1(TΩ) such that supp(φn) ⊂ patch(n) for all n ∈ NΩ, and it is a partition of unity.
We will show that, for each n ∈ NΩ, there exists a subdomain Ωn with the desired properties, such that patch(n) ∩ Ω ⊂ Ωn. We point to Figure 4 for an illustration of the patches Ωn that will be introduced in the proof, for different sets of nodes.
For each c ∈ C, let Nc = {n ∈ NΩ : patch(n) ∩ Ω ⊂ Gc}. There holds For all n ∈ N0, there exists a square Qn such that patch(n) ⊂ B h(T Ω ) (n) ⊂ Qn ⊂ B √ 2h(T Ω ) (n) ⊂ Ω, see Figure 4c. Hence, for all n ∈ N0, we take Ωn := Qn. Define NE := N \ (N0 ∪ NC) = n ∈ N : distΩ(n, c) > dC,1 − h(TΩ), ∀c ∈ C, and distΩ(n, ∂Ω)  0 independent of u such that ue W 1,1 mix (R) ≤ C u W 1,1 mix (Ω) , u f W 1,1 mix (R) ≤ C u W 1,1 mix (Ω) , for all e ∈ E, f ∈ F . Then, by (B.2) and (B.3), v W 1,1 mix ((−1,1) d ) ≤ C u W 1,1 mix (Ω) , for an updated constant C independent of u. This concludes the proof when d = 3. The case d = 2 can be treated by the same argument.
(Ω) := v ∈ L p (Ω) : v W 1,p mix (Ω)  0 such that, for all p = (pi) i∈{1,...,N int } ⊂ N, for all partitions T of I = (0, 1) into Nint open, disjoint, connected subintervals Ii, i = 1, . . . , Nint, for all v ∈ Sp(I, T ) := {v ∈ H 1 (I) : v|I i ∈ Pp i (Ii), i = 1, . . . , Nint}, and for every 0  0 depending on (b − a). Indeed, for any v ∈ H 1 ((a, b)) the concatenation of v with the invertible, affine map T : x → (x − a)/(b − a) is in H 1 ((0,1)). Applying Proposition 3.7 yields NNs {Φ v,T ,p ε } ε∈(0,1) approximating v • T to an appropriate accuracy. Concatenating these networks with the 1-layer NN (A1, b1), where A1x + b1 = T −1 x yields the result. The explicit dependence of C > 0 on (b − a) can be deduced from the error bounds in (0, 1) by affine transformation.
Corollary 4. 1 .
1Let I = (a, b) ⊂ R be a bounded interval. Fix Cp > 0, Cv > 0, and bv > 0. Let 0  0 depending on Cp, Cv, bv, (b − a), and d only such thatL(Φ basis ) ≤ C6(1 + |log ε|)(1 + log(1 + |log ε|)) and M(Φ basis ) ≤ C7dN 1d (1 + |log ε| 2 ).Then, remarking that for all (i1, . . . , i d ) ∈ {1, . . . , N 1d } d there holds E (i 1 ,...,i d ) 0 = d and, by Propositions 3.2, 3.6, and 3.3, we have L(Φε) ≤ C8(1 + |log ε|)(1 + log(1 + |log ε|)), M(Φε) ≤ C9 N d 1d (1 + |log ε|) + M(Φ basis ) .
For
C8, C9 > 0 depending on Cp, Cv, bv, (b − a), d and Cc only. Finally, we conclude that there exists a constant C10 > 0 depending on Cp, Cv, bv, (b − a), d and Cc only such that L(Φε,c) ≤ C10(1 + |log ε|)(1 + log(1 + |log ε|)).
Theorem 4. 3 .
3Let d ∈ {2, 3} and Q := (0, 1) d . Let C = {c} where c is one of the corners of Q and let E = Ec contain the edges adjacent to c when d = 3, E = ∅ when d = 2.
Corollary 4. 4 .
4Let d ∈ {2, 3} and Q := (0, 1) d . Let C = {c} where c is one of the corners of Q and let E = Ec contain the edges adjacent to c when d = 3;
.
Let C = {(0, . . . , 0)} and E = ∅. The regularity of u is a consequence of [29, Theorem 2] (see also [30, Corollary 3.2] for the linear case k = 0): there exists γc > d/2 and Cu, Au > 0 such that u ∈ J γc (Ω; C, E; Cu, Au). Here, γc and the constants Cu and Au depend only on, V0, AV and δ in (5.1), and on k in (5.2).Then, for all 0  0 dependent only on V0, AV , δ, and k, such that
Lemma 5. 5 .
5Let Ω ⊂ R 2 be a polygon with Lipschitz boundary, consisting of straight sides, and with a finite set C of corners. Then, there exists Np ∈ N, a regular triangulation T of R 2 , such that for all K ∈ T either K ⊂ Ω or K ⊂ Ω c . Moreover, there exists a partition of unity {φi}Np i=1 ⊂ [S1(Ω, T )] Np such that [P1] supp(φi) ∩ Ω ⊂ Ωi for all i = 1, . . . , Np, [P2] for each i ∈ {1, . . . , Np}, there exists an affine map ψi : R 2 → R 2 such that ψ −1 i (Ωi) = Ωi for Ωi ∈ {(0, 1) 2 , ΩDN , Ω}, with ΩDN := (−1, 1) × (0, 1), Ω := (−1, 1) 2 \ (−1, 0] 2 ;
|u| H 1 (Ω) + maxi=1,...,N p φi W 1,∞ (Ω) |Ω| 1/2 ) and M×(ε1) := C∞(1 + |log ε1| 4 ). Finally, we set Φε,u := ((1, . . . , 1 Np times
Setting 5. 10 .
10We assume that Ω ⊂ R 2 is an open, bounded polygon with boundary ∂Ω that is Lipschitz and connected. In addition, ∂Ω is the closure of a finite number J ≥ 3 of straight, open sides Γj, i.e., Γi ∩ Γj = ∅ for i = j and ∂Ω = 1≤j≤J Γj. We assume the sides are enumerated cyclically, according to arc length, i.e. ΓJ+1 = Γ1. By nj, we denote the exterior unit normal vector to Ω on Γj and by cj := Γj−1 ∩ Γj the corner j of Ω.With L as in Definition 5.9, we associate on boundary segment Γj a boundary operator Bj ∈ {γ j 0 , γ j 1 }, i.e. either the Dirichlet trace γ0 or the distributional (co-)normal derivative operator γ1, acting on w ∈ C 1 (Ω) via γ j 0 w := w|Γ j , γ j 1 w := (A∇w) · nj|Γ j , j = 1, ..., J .(5.22)We collect the boundary operators Bj in B := {Bj} J j=1 . The first corollary addresses exponential ReLU expressibility of solutions of the source problem corresponding to (L, B).
Corollary 5. 11 .
11Let Ω, L, and B be as in Setting 5.10 with d = 2. For f analytic in Ω, let u denote a solution to the boundary value problem Lu = f in Ω, Bu = 0 on ∂Ω . (5.23)Then, for every 0  0 and let u ∈ H 1 0 (Ω) 2 be the velocity field of the Leray solutions of the viscous, incompressible Navier-Stokes equations in Ω, with homogeneous Dirichlet ("no slip") boundary conditions− ν∆u + (u · ∇)u + ∇p = f in Ω, ∇ · u = 0 in Ω, u = 0 on ∂Ω,(5.27)where the components of f are analytic in Ω and such that f H −1 (Ω) /ν 2 is small enough so that u is unique. Then, for every 0  1, for all c ∈ C if d = 2, γ = {γc, γe : c ∈ C, e ∈ E}, with γc > 3/2 and γe > 1, for all c ∈ C and e ∈ E if d = 3. Then, for any 0  1 for all c ∈ CDN .Then, for any 0  0 independent of u such thatũ W 1,1 mix ((−1,1) d ) ≤ Cext u W 1,1 mix (Ω DN ) . (5.35)Because u ∈ J γ (ΩDN ; CDN , ∅), it holds with CS = S ∩ CDN that u ∈ J γ (S; CS, ∅) for allThe remaining steps are the same as those in the proof of Theorem 5.16.
For positive integers p and s such that 1 ≤ s ≤ p,
Lemma A. 4 ([ 45 ,
445Propositions 5.2 and 5.3]). Let d = 3, (p1, p2, p3) ∈ N 3 , and (s1, s2, s3) ∈ N 3 with 1 ≤ si ≤ pi. Then the projector Πp 1 p 2 p 3 :
) + (II) + (III) .(A.30)
s} Ψp,sA 2s+4 ((s + 3)!) 2 , (A.46b) and the estimate on Kc is the same as that on Kce. The assertion now follows from (A.42), (A.43), (A.45), and (A.46), upon possibly adjusting the value of the constant A.
γc > 3/2 and γe > 1, for all c ∈ C and e ∈ E, if d = 3, there exists a constant C0 > 0 independent of such that if p ≥ C0 , there exist constants C, b > 0 such that v − Π ,p hp,d v H 1 (K) ≤ Ce −b .
(A.51) from (A.54), (A.58), (A.62), and (A.65). Furthermore, (A.52) follows from (A.55), (A.59), (A.63), and (A.66). The estimates for the case d = 2 follow from the same argument. The following lemma shows the continuous imbedding of J d γ (Q; C, E) into W 1,1 mix (Q), given sufficiently large weights γ. Lemma A.22. Let d ∈ {2, 3}. Let γ be such that γc > d/2, for all c ∈ C and (if d = 3) γe > 1 for all e ∈ E.
...,i d ci 1 ...i d φi 1 ...i d and, for all (i1, . . . , i d ) ∈ {1, . . . , N 1d } d , there exists vi j , j = 1, . . . , d such that φi 1 ...i d = d j=1 vi j , see Section A.9.2 and Proposition A.24. By construction of vi in (A.71), and by using (A.47) and (A.48), we observe that vi L ∞ (I) ≤ 1 for all i = 1, . . . , N 1d . In addition, (A.50), demonstrates that vi H 1 (I) ≤ 2 |supp(vi)| 1/2 deg(vi) 1/2 ≤ 2σ −L/2 ∀i ∈ {1, . . . , N 1d }.
Let a > 0 and Ω = (−a, a) d . Denote the set of corners set of edges EΩ = ∅ if d = 2, and, if d the affine transformations ψ1,+ : (0, 1) → (0, a/2), ψ2,+ : (0, 1) → (a/2, a), ψ1,− : (0, 1) → (−a/2, 0), ψ2,− : (0, 1) → (−a, −a/2) such that
Figure 3 :
3Multipatch geometric tensor product meshes G d , for d = 2 (left) and d = 3 (right).
Theorem A. 25 .
25For a > 0, let Ω = (−a, a) d , d = 2, 3. Denote by Ω k , k = 1, . . . , 4 d the patches composingΩ, i.e., the sets Ω k = × d j=1 (a k j , a k j + a/2) with a k j ∈ {−a, −a/2, 0, a/2}. Denote also C k = CΩ ∩ Ω k and E k = {e ∈ EΩ : e ⊂ Ω k }, which contain one singular corner, and three singular edges abutting that corner, as in(A.1) and (A.2). Let I ⊂ {1, . . . , 4 d } and let v ∈ W 1,1 mix (Ω) be such that, for all k ∈ I, there holds v |Ω k ∈ J γ k (Ω k ; C k , E k ) with γ k c > 1, for all c ∈ C k , if d = 2,γ k c > 3/2 and γ k e > 1, for all c ∈ C k and e ∈ E k , if d = 3.Then, there exist constants cp > 0 and C, b > 0 such that, for all ∈ N,with p = cp , v − Π ,p hp,d v H 1 (Ω k ) ≤ Ce −b ≤ C exp(−b 2d √ N dof ). (A.77)Here, N dof = O( 2d ) denotes the overall number of degrees of freedom in the piecewise polynomial approximation. Furthermore, writing N 1d = 4( + 1)p + 1, there exists an array of coefficientsc = ci 1 ...i d : (i1, . . . , i d ) ∈ {1, . . . , N 1d } d such that Π ,p hp,d v (x1, . . . , x d ) = N 1d i 1 ,...,i d =1 ci 1 ...i d d j=1ṽ i j (xj) ∀(x1, . . . , x d ) ∈ Ω,
vi H 1 ((−a,a)) ≤ C1σ − /2 ∀i = 1, . . . , N 1d , (A.78) and N 1d i 1 ,...,i d =1 | ci 1 ...i d | ≤ C2 d j=0 ( + 1) j (p + 1) 2(d−j) v W 1,1 mix (Ω) . (A.79)
d H 1 (Ω k ) ≤ Ce −b . Equation (A.77) follows. The bounds (A.78) and (A.79) follow from the construction of the basis functions (A.47)-(A.48) and from the application of Lemma A.21 in each patch, respectively. B Proofs of Section 5 B.1 Proof of Lemma 5.5
Furthermore, for each e ∈ E, denote de := ∞ if Ω is a triangle, otherwise de := min {distΩ(e, e1) : e1 ∈ E and e ∩ e1 = ∅} , dE = min e∈E de. Finally, for all x ∈ Ω, let ne(x) := #{e1, e2, . . . ∈ E : distΩ(x, ∂Ω) = distΩ(x, e1) = distΩ(x, e2) = . . . }.
Figure 4 :
4Patches Ω n for nodes near a convex (a), nonconvex corner (b), for nodes in the interior of Ω (c), and near an edge (d).
Nc := {n ∈ NΩ : distΩ(n, c) ≤ dC,1 − h(TΩ)} ⊂ Nc.Therefore, all the nodes n ∈ Nc are such that patch(n) ∩ Ω ⊂ Gc =: Ωn. Denote thenNC = c∈C Nc. Note that, due to (B.1), there holds √ 2h(TΩ) ≤ √ 2 √ 2+1 dC,1 ≤ dC,1 − h(TΩ). We consider the nodes in N \ NC. First, consider the nodes in N0 := {n ∈ N \ NC : distΩ(n, ∂Ω) ≥ √ 2h(TΩ)}.
To show that n ∈ NC ∪ N0, it suffices to show that dist(x0, ei) ≥ √ 2h(TΩ) for i = 1, 2. Because
Figure 5 :
5dC,1 ≤ dC,1 − h(TΩ), it a fortiori holds for all points y in Ω on the bisector intersected withB d C,1 −h(T Ω ) (c) c , that dist(y, ei) ≥ √ 2h(TΩ), which shows that if distΩ(n, c) ≥ dC,1 − h(TΩ), then n ∈ N0. If c is a nonconvex corner, then dist(x0, ei) ≥ √ 2h(TΩ) for i = 1, 2 follows immediately from dist(x0, ei) = dist(x0, c) =√ 2 √ 2+1dC,1 and (B.1). To show that dist(x0, ei) ≥ √ 2h(TΩ), i = 1, 2 in case c is a convex corner, we make the following definitions (seeFigure 5): Situation near a convex corner c.• For i = 1, 2, let xi be the intersection of ei and B
with N 1d = (L + 1)p + 1, and, for all f ∈ J γ (Q; C, E; C f , A f ) there exists a d-dimensional array of coefficientsc = ci 1 ...i d : (i1, . . . , i d ) ∈ {1, . . . , N 1d } d suchthat 1. For every i = 1, . . . N 1d , supp(vi) intersects either a single interval or two neighboring subintervals of G L 1 . Furthermore, there exist constants Cv, bv depending only on C f , A f , σ such that vi ∞ ≤ 1, vi H 1 (I) ≤ Cvε −bv , ∀i = 1, . . . , N 1d .so that there exist piecewise polynomials 
vi ∈ Qp(G L 
1 ) ∩ H 1 (I), 
i = 1, . . . , N 1d , 
(2.3) 
35, Theorem 2.1], (5.8), and Proposition 3.6, for all i ∈ {1, . . . , N },
Proof. By [7, Corollary 7.1, Theorems 7.3 and 7.4] if d = 3 and [2, Theorem 3.1] if d = 2, there exists γ such that γc − d/2 > 0 for all c ∈ C and γe > 1 for all e ∈ E such that u ∈ J γ (Ω; C, E). An application of Theorem 5.16 concludes the proof. PINNs' methodology. There, the loss function was based on minimization of the residual of the NN approximation in the strong form of the PDE. Evidently, a different (smoother) activation than the ReLU activations considered here had to be used. Starting from the approximation of products by NNs with smoother activation functions introduced in [46, Sec.3.3] and following the same line of reasoning as in the present paper, the results we obtain for ReLU-based realizations of NNs can be extended to large classes of NNs with smoother activations and similar architecture.Remark 5.19. By [7, Corollary 7.1 and Theorem 7.4], Corollary 5.18 holds verbatim also under the hypothesis 
that the right-hand side f is weighted analytic, with singularities at the corners/edges of the domain; specifically, 
(5.33) and the size bounds on the NN Φε,u hold under the assumption that there exists γ such that γc − d/2 > 0 
for all c ∈ C and γe > 1 for all e ∈ E such that 
f ∈ J γ−2 (Ω; C, E). 
Remark 5.20. The numerical approximation of solutions for (5.29) with a NN in two dimensions has been 
investigated e.g. in [28] using the so-called 'Furthermore, in [8, Section 3.1
d ). The continuity can also be observed from the expansion in terms of continuous, globally defined basis functions given in Proposition A.24.    Remark A.3. The projector Π ,p 
hp,d is defined on a larger space than H 1 
mix (Q) as specified below (e.g. Remark 
A.20). 
x3) and consider ti 1 ...i d = #{ij ≤ 2, j = 1, 2, 3}. We have • if ti 1 ...i d = 1, then ψ K i 1 ...i d is not zero only on one face of the boundary of K, • if ti 1 ...i d = 2, then ψ K i 1 ...i d is not zero only on one edge and neighboring faces of the boundary of K, • if ti 1 ...i d = 3, then ψ K i 1 ...i d is not zero only on one corner and neighboring edges and faces of the boundary of K. Similar arguments hold when d = 2.
1 / 2
/∀i ∈ N0, ∀k ∈ {0, . . . , }. (A.50) Remark A.20. As mentioned in Remark A.3, the hp-projector Π ,p hp,d can be defined for more general functions than u ∈ H 1 mix (Q). As follows from Equations (A.53), (A.57), (A.61) and (A.64) below, the projector is also defined for u ∈ W 1,1 mix (Q). Lemma A.21. There exist constants C1, C2 such that, for all u ∈ W 1,1 mix (Q), all ∈ N, all p ∈ N |c K i 1 ...i d | ≤ C ∀K ∈ G d , ∀(i1, . . . , i d ) ∈ {1, . . . , p + 1} d (A.51) and for all (i1, . . . , i d ) ∈ {1, . . . , p + 1} d K∈G 3d 
j=1 
ij 
u W 1,1 
mix (Q) 
for all K ∈ G 3 and all i1, . . . , i d ∈ {1, . . . , p + 1}. Corollary A.23. Let γ be such that γc − d/2 > 0, for all c ∈ C and, if d = 3, γe > 1 for all e ∈ E. There exists a constant C > 0 such that for all , p ∈ N and for all u ∈ J d γ (Q; C, E),Π ,p hp,d u L ∞ (Q) ≤ Cp 2d u J d γ (Q) .Proof of Theorem 2.1. Fix A f , C f , and γ as in the hypotheses. Then, by Proposition A.19, there exists cp, C hp , b hp > 0 such that for every ∈ N and for all v ∈ J γ (Q; C, E; C f , A f ), there exists v hp ∈ X ,cp hp,d such that (see Section A.1 for the definition of the space XA.9.3 Proof of Theorem 2.1 
,cp 
76 )
76Then, the global hp projection operator Π ,p hp,d : W 1,1 mix (Ω) → X ,p hp,d is defined as
We assume isotropic tensorization, i.e. the same σ and the same number of geometric mesh layers in each coordinate direction; all approximation results remain valid (with possibly better numerical values for the constants in the error bounds) for anisotropic, co-ordinate dependent choices of and of σ.
A.9.2 Basis of continuous functions with compact supportIt is possible to construct a basis for Π ,p hp,d in Q such that all basis functions are continuous and have compact support. For all ∈ N and all p ∈ N, extend to zero outside of their domain of definition the functions ζ k n defined in (A.47) and (A.48), for k = 0, . . . , and n = 1, . . . , p + 1. We introduce the univariate functions with compact support vj : (0, 1) → R, for j = 1, . . . ,and v +2+k(p−1)+n = ζ k n+2 , for all k = 0, . . . , and n = 1, . . . , p − 1.There exists an array of coefficientsFurthermore, there exist constants C1, C2 > 0 independent of , p, and u, such that
The h-p version of the finite element method for domains with curved boundaries. I Babuška, B Guo, SIAM J. Numer. Anal. 254I. Babuška and B. Guo. The h-p version of the finite element method for domains with curved boundaries. SIAM J. Numer. Anal., 25(4):837-861, 1988.
Regularity of the solution of elliptic problems with piecewise analytic data. I. Boundary value problems for linear elliptic equation of second order. I Babuška, B Guo, SIAM J. Math. Anal. 191I. Babuška and B. Guo. Regularity of the solution of elliptic problems with piecewise analytic data. I. Boundary value problems for linear elliptic equation of second order. SIAM J. Math. Anal., 19(1):172-203, 1988.
Regularity and numerical solution of eigenvalue problems with piecewise analytic data. I Babuška, B Q Guo, J E Osborn, SIAM J. Numer. Anal. 266I. Babuška, B. Q. Guo, and J. E. Osborn. Regularity and numerical solution of eigenvalue problems with piecewise analytic data. SIAM J. Numer. Anal., 26(6):1534-1560, 1989.
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. J Berner, P Grohs, A Jentzen, SIAM J. Math. Data Sci. 23J. 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. SIAM J. Math. Data Sci., 2(3):631- 657, 2020.
Régularité Gevrey pour le problème de Dirichlet dans des domainesà singularités coniques. P Bolley, M Dauge, J Camus, Comm. Partial Differential Equations. 104P. Bolley, M. Dauge, and J. Camus. Régularité Gevrey pour le problème de Dirichlet dans des domainesà singularités coniques. Comm. Partial Differential Equations, 10(4):391-431, 1985.
The mathematical theory of finite element methods. S C Brenner, L R Scott, Texts in Applied Mathematics. 15Springerthird editionS. C. Brenner and L. R. Scott. The mathematical theory of finite element methods, volume 15 of Texts in Applied Mathematics. Springer, New York, third edition, 2008.
Analytic regularity for linear elliptic systems in polygons and polyhedra. M Costabel, M Dauge, S Nicaise, Math. Models Methods Appl. Sci. 22863M. Costabel, M. Dauge, and S. Nicaise. Analytic regularity for linear elliptic systems in polygons and polyhedra. Math. Models Methods Appl. Sci., 22(8):1250015, 63, 2012.
The deep Ritz method: a deep learning-based numerical algorithm for solving variational problems. W E , B Yu, Commun. Math. Stat. 61W. E and B. Yu. The deep Ritz method: a deep learning-based numerical algorithm for solving variational problems. Commun. Math. Stat., 6(1):1-12, 2018.
DNN expression rate analysis of highdimensional PDEs: Application to option pricing. D Elbrächter, P Grohs, A Jentzen, C Schwab, 1809.07669arXiv preprintD. Elbrächter, P. Grohs, A. Jentzen, and C. Schwab. DNN expression rate analysis of high- dimensional PDEs: Application to option pricing. arXiv preprint 1809.07669, 2018.
Asymptotic regularity of solutions to Hartree-Fock equations with Coulomb potential. H.-J Flad, R Schneider, B.-W Schulze, Math. Methods Appl. Sci. 3118H.-J. Flad, R. Schneider, and B.-W. Schulze. Asymptotic regularity of solutions to Hartree-Fock equations with Coulomb potential. Math. Methods Appl. Sci., 31(18):2172-2201, 2008.
Analytic structure of solutions to multiconfiguration equations. S Fournais, M Hoffmann-Ostenhof, T Hoffmann-Ostenhof, T Østergaard Sørensen, J. Phys. A. 423111S. Fournais, M. Hoffmann-Ostenhof, T. Hoffmann-Ostenhof, and T. Østergaard Sørensen. Analytic structure of solutions to multiconfiguration equations. J. Phys. A, 42(31):315208, 11, 2009.
Caratterizzazioni delle tracce sulla frontiera relative ad alcune classi di funzioni in n variabili. E Gagliardo, Rend. Sem. Mat. Univ. Padova. 27E. Gagliardo. Caratterizzazioni delle tracce sulla frontiera relative ad alcune classi di funzioni in n variabili. Rend. Sem. Mat. Univ. Padova, 27:284-305, 1957.
Error bounds for approximations with deep ReLU neural networks in W s,p norms. I Gühring, G Kutyniok, P Petersen, Anal. Appl. (Singap.). 1805I. Gühring, G. Kutyniok, and P. Petersen. Error bounds for approximations with deep ReLU neural networks in W s,p norms. Anal. Appl. (Singap.), 18(05):803-859, 2020.
On the regularity of elasticity problems with piecewise analytic data. B Guo, I Babuška, Adv. in Appl. Math. 143B. Guo and I. Babuška. On the regularity of elasticity problems with piecewise analytic data. Adv. in Appl. Math., 14(3):307-347, 1993.
Regularity of the solutions for elliptic problems on nonsmooth domains in R 3 . I. Countably normed spaces on polyhedral domains. B Guo, I Babuška, Proc. Roy. Soc. Edinburgh Sect. A. 1271B. Guo and I. Babuška. Regularity of the solutions for elliptic problems on nonsmooth domains in R 3 . I. Countably normed spaces on polyhedral domains. Proc. Roy. Soc. Edinburgh Sect. A, 127(1):77-126, 1997.
Regularity of the solutions for elliptic problems on nonsmooth domains in R 3 . II. Regularity in neighbourhoods of edges. B Guo, I Babuška, Proc. Roy. Soc. Edinburgh Sect. A. 1273B. Guo and I. Babuška. Regularity of the solutions for elliptic problems on nonsmooth domains in R 3 . II. Regularity in neighbourhoods of edges. Proc. Roy. Soc. Edinburgh Sect. A, 127(3):517-545, 1997.
Analytic regularity of Stokes flow on polygonal domains in countably weighted Sobolev spaces. B Guo, C Schwab, J. Comput. Appl. Math. 1901-2B. Guo and C. Schwab. Analytic regularity of Stokes flow on polygonal domains in countably weighted Sobolev spaces. J. Comput. Appl. Math., 190(1-2):487-519, 2006.
Solving many-electron Schrödinger equation using deep neural networks. J Han, L Zhang, W E , J. Comput. Phys. 399108929J. Han, L. Zhang, and W. E. Solving many-electron Schrödinger equation using deep neural networks. J. Comput. Phys., 399:108929, 8, 2019.
ReLU deep neural networks and linear finite elements. J He, L Li, J Xu, C Zheng, J. Comp. Math. 382020J. He, L. Li, J. Xu, and C. Zheng. ReLU deep neural networks and linear finite elements. J. Comp. Math., 38, 2020.
Deep neural network solution of the electronic Schrödinger equation. J Hermann, Z Schätzle, F Noé, arXiv:1909.08423arXiv preprintJ. Hermann, Z. Schätzle, and F. Noé. Deep neural network solution of the electronic Schrödinger equation. arXiv preprint arXiv:1909.08423, 2019.
A proof that deep artificial neural networks overcome the curse of dimensionality in the numerical approximation of Kolmogorov partial differential equations with constant diffusion and nonlinear drift coefficients. A Jentzen, D Salimova, T Welti, arXiv:1809.07321arXiv preprintA. Jentzen, D. Salimova, and T. Welti. A proof that deep artificial neural networks overcome the curse of dimensionality in the numerical approximation of Kolmogorov partial differential equations with constant diffusion and nonlinear drift coefficients. arXiv preprint arXiv:1809.07321, 2018.
QTT-Finite-Element approximation for multiscale problems I: model problems in one dimension. V Kazeev, I Oseledets, M Rakhuba, C Schwab, Adv. Comput. Math. 432V. Kazeev, I. Oseledets, M. Rakhuba, and C. Schwab. QTT-Finite-Element approximation for multiscale problems I: model problems in one dimension. Adv. Comput. Math., 43(2):411-442, 2017.
Quantized tensor-structured finite elements for second-order elliptic PDEs in two dimensions. V Kazeev, C Schwab, Numer. Math. 1381V. Kazeev and C. Schwab. Quantized tensor-structured finite elements for second-order elliptic PDEs in two dimensions. Numer. Math., 138(1):133-190, 2018.
Efficient approximation of solutions of parametric linear transport equations by ReLU DNNs. F Laakmann, P Petersen, arXiv:2001.11441arXiv preprintF. Laakmann and P. Petersen. Efficient approximation of solutions of parametric linear transport equations by ReLU DNNs. arXiv preprint arXiv:2001.11441, 2020.
The Hartree-Fock theory for Coulomb systems. E H Lieb, B Simon, Comm. Math. Phys. 533E. H. Lieb and B. Simon. The Hartree-Fock theory for Coulomb systems. Comm. Math. Phys., 53(3):185-194, 1977.
Solutions of Hartree-Fock equations for Coulomb systems. P.-L Lions, Comm. Math. Phys. 1091P.-L. Lions. Solutions of Hartree-Fock equations for Coulomb systems. Comm. Math. Phys., 109(1):33-97, 1987.
Deep network approximation for smooth functions. J Lu, Z Shen, H Yang, S Zhang, arXiv:2001.03040arXiv preprintJ. Lu, Z. Shen, H. Yang, and S. Zhang. Deep network approximation for smooth functions. arXiv preprint arXiv:2001.03040, 2020.
DeepXDE: A deep learning library for solving differential equations. L Lu, X Meng, Z Mao, G E Karniadakis, arXiv:1907.04502arXiv e-printsL. Lu, X. Meng, Z. Mao, and G. E. Karniadakis. DeepXDE: A deep learning library for solving differential equations. arXiv e-prints, page arXiv:1907.04502, July 2019.
Y Maday, C Marcati, arXiv:1912.07483Analyticity and hp discontinuous Galerkin approximation of nonlinear Schrödinger eigenproblems. arXiv e-prints. Y. Maday and C. Marcati. Analyticity and hp discontinuous Galerkin approximation of nonlinear Schrödinger eigenproblems. arXiv e-prints, page arXiv:1912.07483, Dec. 2019.
Regularity and hp discontinuous Galerkin finite element approximation of linear elliptic eigenvalue problems with singular potentials. Y Maday, C Marcati, Math. Models Methods Appl. Sci. 298Y. Maday and C. Marcati. Regularity and hp discontinuous Galerkin finite element approximation of linear elliptic eigenvalue problems with singular potentials. Math. Models Methods Appl. Sci., 29(8):1585-1617, 2019.
Weighted analyticity of Hartree-Fock eigenfunctions. Y Maday, C Marcati, 2020-59Seminar for Applied Mathematics. Technical ReportY. Maday and C. Marcati. Weighted analyticity of Hartree-Fock eigenfunctions. Technical Report 2020-59, Seminar for Applied Mathematics, ETH Zürich, Switzerland, 2020.
Tensor Rank bounds for Point Singularities in R 3. C Marcati, M Rakhuba, C Schwab, 2019-68Seminar for Applied Mathematics. Technical ReportC. Marcati, M. Rakhuba, and C. Schwab. Tensor Rank bounds for Point Singularities in R 3 . Technical Report 2019-68, Seminar for Applied Mathematics, ETH Zürich, Switzerland, 2019.
Analytic regularity for the incompressible Navier-Stokes equations in polygons. C Marcati, C Schwab, SIAM J. Math. Anal. 523C. Marcati and C. Schwab. Analytic regularity for the incompressible Navier-Stokes equations in polygons. SIAM J. Math. Anal., 52(3):2945-2968, 2020.
Elliptic Equations in Polyhedral Domains. V Maz&apos;ya, J Rossmann, Mathematical Surveys and Monographs. 162American Mathematical SocietyV. Maz'ya and J. Rossmann. Elliptic Equations in Polyhedral Domains, volume 162 of Mathematical Surveys and Monographs. American Mathematical Society, Providence, Rhode Island, April 2010.
The partition of unity finite element method: basic theory and applications. J M Melenk, I Babuška, Comput. Methods Appl. Mech. Engrg. 1391-4J. M. Melenk and I. Babuška. The partition of unity finite element method: basic theory and applications. Comput. Methods Appl. Mech. Engrg., 139(1-4):289-314, 1996.
Deep ReLU networks and high-order finite element methods. J A A Opschoor, P C Petersen, C Schwab, Analysis and Applications. 1805J. A. A. Opschoor, P. C. Petersen, and C. Schwab. Deep ReLU networks and high-order finite element methods. Analysis and Applications, 18(05):715-770, 2020.
Exponential ReLU DNN expression of holomorphic maps in high dimension. J A A Opschoor, C Schwab, J Zech, 2019-35 rev1Seminar for Applied Mathematics. Technical ReportAccepted for publication in Constructive ApproximationJ. A. A. Opschoor, C. Schwab, and J. Zech. Exponential ReLU DNN expression of holomorphic maps in high dimension. Technical Report 2019-35 rev1, Seminar for Applied Mathematics, ETH Zürich, Switzerland, 2020. Accepted for publication in Constructive Approximation.
Optimal approximation of piecewise smooth functions using deep ReLU neural networks. P Petersen, F Voigtlaender, Neural Netw. 108P. Petersen and F. Voigtlaender. Optimal approximation of piecewise smooth functions using deep ReLU neural networks. Neural Netw., 108:296 -330, 2018.
Ab-initio solution of the many-electron Schrödinger equation with deep neural networks. D Pfau, J S Spencer, A G , G Matthews, W M C Foulkes, arXiv:1909.02487D. Pfau, J. S. Spencer, A. G. de G. Matthews, and W. M. C. Foulkes. Ab-initio solution of the many-electron Schrödinger equation with deep neural networks, 2019. arXiv:1909.02487.
Hidden physics models: machine learning of nonlinear partial differential equations. M Raissi, G E Karniadakis, J. Comput. Phys. 357M. Raissi and G. E. Karniadakis. Hidden physics models: machine learning of nonlinear partial differential equations. J. Comput. Phys., 357:125-141, 2018.
Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. M Raissi, P Perdikaris, G E Karniadakis, J. Comput. Phys. 378M. Raissi, P. Perdikaris, and G. E. Karniadakis. Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. J. Comput. Phys., 378:686-707, 2019.
Exponential convergence for hp-version and spectral finite element methods for elliptic problems in polyhedra. D Schötzau, C Schwab, Math. Models Methods Appl. Sci. 259D. Schötzau and C. Schwab. Exponential convergence for hp-version and spectral finite element methods for elliptic problems in polyhedra. Math. Models Methods Appl. Sci., 25(9):1617-1661, 2015.
Exponential convergence of hp-FEM for elliptic problems in polyhedra: mixed boundary conditions and anisotropic polynomial degrees. D Schötzau, C Schwab, Found. Comput. Math. 183D. Schötzau and C. Schwab. Exponential convergence of hp-FEM for elliptic problems in polyhedra: mixed boundary conditions and anisotropic polynomial degrees. Found. Comput. Math., 18(3):595- 660, 2018.
hp-dGFEM for second-order elliptic problems in polyhedra I: Stability on geometric meshes. D Schötzau, C Schwab, T P Wihler, SIAM J. Numer. Anal. 513D. Schötzau, C. Schwab, and T. P. Wihler. hp-dGFEM for second-order elliptic problems in polyhedra I: Stability on geometric meshes. SIAM J. Numer. Anal., 51(3):1610-1633, 2013.
hp-DGFEM for second order elliptic problems in polyhedra II: Exponential convergence. D Schötzau, C Schwab, T P Wihler, SIAM J. Numer. Anal. 514D. Schötzau, C. Schwab, and T. P. Wihler. hp-DGFEM for second order elliptic problems in polyhedra II: Exponential convergence. SIAM J. Numer. Anal., 51(4):2005-2035, 2013.
Deep learning in high dimension: Neural network expression rates for generalized polynomial chaos expansions in UQ. C Schwab, J Zech, Anal. Appl. (Singap.). 171C. Schwab and J. Zech. Deep learning in high dimension: Neural network expression rates for generalized polynomial chaos expansions in UQ. Anal. Appl. (Singap.), 17(1):19-55, 2019.
PFNN: A Penalty-Free Neural Network Method for Solving a Class of Second-Order Boundary-Value Problems on Complex Geometries. H Sheng, C Yang, arXiv:2004.06490arXiv e-printsH. Sheng and C. Yang. PFNN: A Penalty-Free Neural Network Method for Solving a Class of Second- Order Boundary-Value Problems on Complex Geometries. arXiv e-prints, page arXiv:2004.06490, Apr. 2020.
Y Shin, Z Zhang, G E Karniadakis, arXiv:2010.08019Error estimates of residual minimization using neural networks for linear PDEs. arXiv e-prints. Y. Shin, Z. Zhang, and G. E. Karniadakis. Error estimates of residual minimization using neural networks for linear PDEs. arXiv e-prints, page arXiv:2010.08019, Oct. 2020.
DGM: a deep learning algorithm for solving partial differential equations. J Sirignano, K Spiliopoulos, J. Comput. Phys. 375J. Sirignano and K. Spiliopoulos. DGM: a deep learning algorithm for solving partial differential equations. J. Comput. Phys., 375:1339-1364, 2018.
Adaptivity of deep ReLU network for learning in Besov and mixed smooth Besov spaces: optimal rate and curse of dimensionality. T Suzuki, 7th International Conference on Learning Representations, ICLR 2019. New Orleans, LA, USAT. Suzuki. Adaptivity of deep ReLU network for learning in Besov and mixed smooth Besov spaces: optimal rate and curse of dimensionality. In 7th International Conference on Learning Representations, ICLR 2019, New Orleans, LA, USA, May 6-9, 2019, 2019.
Region configurations for realizability of lattice piecewise-linear models. J M Tarela, M V Martínez, Math. Comput. Modelling. 30J. M. Tarela and M. V. Martínez. Region configurations for realizability of lattice piecewise-linear models. Math. Comput. Modelling, 30(11-12):17-27, 1999.
Error bounds for approximations with deep ReLU networks. D Yarotsky, Neural Netw. 94D. Yarotsky. Error bounds for approximations with deep ReLU networks. Neural Netw., 94:103 - 114, 2017.
Optimal approximation of continuous functions by very deep ReLU networks. D Yarotsky, PMLRProceedings of Machine Learning Research. Machine Learning ResearchD. Yarotsky. Optimal approximation of continuous functions by very deep ReLU networks. vol- ume 75 of Proceedings of Machine Learning Research, pages 639-649. PMLR, 06-09 Jul 2018.
Exponential ReLU Neural Network Approximation Rates for Point and Edge Singularities
https://arxiv.org/abs/2010.12217
https://stats.stackexchange.com/questions/532309/why-pure-exponent-is-not-used-as-activation-function-for-neural-networks
Why pure exponent is not used as activation function for neural networks?
2021-06-26
MefAldemisov MefAldemisov
<div><div>
 
<div>
 
 
<div>
 <p>$\begingroup$</p><div>
 
<p>The ReLU function is commonly used as an activation function in machine learning, as well, as its modifications (ELU, leaky ReLU).</p>
<p>The overall idea of these functions is the same: before <code>x = 0</code> the value of the function is small (its limit to infinity is zero or <code>-1</code>), after <code>x = 0</code> the function grows proportionally to x.</p>
<p>The exponent function (<code>e^x</code> or <code>e^x-1</code>) has similar behavior, and its derivative in <code>x = 0</code> is greater than for sigmoid.</p>
<p>The visualization below illustrates the exponent in comparison with ReLU and sigmoid activation functions.</p>
<p><a href="https://i.stack.imgur.com/tT6no.png"></a></p>
<p>So, why the simple function <code>y=e^x</code> is not used as an activation function in neural networks?</p>
 </div>
 
 <div>
 <p>
 asked Jun 26, 2021 at 16:18
 </p>
 <div>
 <a href="/users/326632/mefaldemisov"><p></p></a>
 </div>
 
</div>
 <p>$\endgroup$
</p></div>
 <p>3</p>
 </div>
 <div>
 
 
<div>
 
 
<div>
 <p>$\begingroup$</p><div>
<p>I think the most prominent reason is stability. Think about having consequent layers with exponential activation, and what happens to the output when you input a small number to the NN (e.g. $x=1$), the forward calculation will look like:
$$o=\exp(\exp(\exp(\exp(1))))\approx e^{3814279}$$</p>
<p>It can go crazy very quickly and I don't think you can train deep networks with this activation function unless you add other mechanisms like clipping.</p>
 </div>
 <div>
 <p>
 answered Jun 26, 2021 at 16:47
 </p>
 <div>
 <a href="/users/204068/gunes"><p></p></a>
 </div>
 <div>
 <p><a href="/users/204068/gunes">gunes</a>gunes</p><p>54.3k4 gold badges45 silver badges80 bronze badges
 </p>
 </div>
</div>
 <p>$\endgroup$
</p></div>
 <p>4</p>
 </div>
 
 <h2>
 
 </h2>
 </div>
 </div>||||I|||| Stack Exchange Network

 Stack Exchange network consists of 181 Q&A communities including Stack Overflow, the largest, most trusted online community for developers to learn, share their knowledge, and build their careers.

 Visit Stack Exchange
 Loading…
 1. 
 + Tour Start here for a quick overview of the site
 + Help Center Detailed answers to any questions you might have
 + Meta Discuss the workings and policies of this site
 + About Us Learn more about Stack Overflow the company
 + Business Learn more about our products
 2. 

 3. current community

 + Cross Validated
 help chat
 + Cross Validated Meta

 your communities

 Sign up or log in to customize your list.

 more stack exchange communities

 company blog
 4. 
 5. Log in
 6. Sign up

 Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It only takes a minute to sign up.

 Sign up to join this community
 Anybody can ask a question
 Anybody can answer
 The best answers are voted up and rise to the top
 1. Home
 2. 
 1. Public
 2. Questions
 3. Tags
 4. Users
 5. Unanswered
 3. 
 1. Teams
 Stack Overflow for Teams – Start collaborating and sharing organizational knowledge. Create a free Team Why Teams?
 2. Teams
 3. Create free Team

 Teams

 Q&A for work

 Connect and share knowledge within a single location that is structured and easy to search.

 Learn more about Teams

 Why pure exponent is not used as activation function for neural networks?

 Ask Question
 Asked 1 year, 7 months ago
 Modified 1 year, 7 months ago
 Viewed 4k times
 17
 $\begingroup$

 The ReLU function is commonly used as an activation function in machine learning, as well, as its modifications (ELU, leaky ReLU).

 The overall idea of these functions is the same: before x = 0 the value of the function is small (its limit to infinity is zero or -1), after x = 0 the function grows proportionally to x.

 The exponent function (e^x or e^x-1) has similar behavior, and its derivative in x = 0 is greater than for sigmoid.

 The visualization below illustrates the exponent in comparison with ReLU and sigmoid activation functions.

 So, why the simple function y=e^x is not used as an activation function in neural networks?

 * machine-learning
 * neural-networks
 * perceptron
 Share
 Cite
 Improve this question
 Follow
 asked Jun 26, 2021 at 16:18
 MefAldemisov MefAldemisov
 273 2 2 silver badges 6 6 bronze badges
 $\endgroup$
 3
 * 3
 $\begingroup$ You can . There is nothing stopping you. If you are interested in a smooth variant of the ReLU function, why not use the Softplus function? $\endgroup$
 – mhdadk
 Jun 26, 2021 at 16:33
 * 3
 $\begingroup$ @mhdadk, I am not asking, if I can do this or not. I'm asking about the drawback(s) of the relatively naive solution to the problem of activation functions. $\endgroup$
 – MefAldemisov
 Jun 26, 2021 at 16:38
 * 3
 $\begingroup$ If you view Neural Networks in the abstract as a sort of black-box tool for function approximation, a good toy problem in this domain is to consider the convergence of continued exponentials. Bender and Vinson wrote about this in the context of pertubation theory, where you use only a small number of terms of the Taylor expansion. The upshot was that the convergence, while quite fascinating, has an intricate structure and doesn't work everywhere, which explains a bit why this kind of thing might not be suitable as an activation function. $\endgroup$
 – Cade Reinberger
 Jun 28, 2021 at 3:29
 Add a comment |

 1 Answer

 Sorted by: Reset to default
 Highest score (default) Date modified (newest first) Date created (oldest first)
 31
 $\begingroup$

 I think the most prominent reason is stability. Think about having consequent layers with exponential activation, and what happens to the output when you input a small number to the NN (e.g. $x=1$ ), the forward calculation will look like: $$o=\exp(\exp(\exp(\exp(1))))\approx e^{3814279}$$

 It can go crazy very quickly and I don't think you can train deep networks with this activation function unless you add other mechanisms like clipping.

 Share
 Cite
 Improve this answer
 Follow
 edited Jun 26, 2021 at 19:45
 answered Jun 26, 2021 at 16:47
 gunes gunes
 54.3k 4 4 gold badges 45 45 silver badges 80 80 bronze badges
 $\endgroup$
 4
 * 5
 $\begingroup$ Well, when using floating point this quickly becomes an “indicator function” $f(x) = \begin{cases}0 & \text{if $x\leq 0$} \\ \infty &\text{otherwise.}\end{cases}$ That kind of “function” can actually be used with nonsmooth convex-optimisation algorithms like Douglas-Rachford (which don't really operate on the function itself or its gradient, but instead on proximals and convex conjugates ). For GD-based nonlinear optimisation as used in deep learning, this is not an option. $\endgroup$
 – leftaroundabout
 Jun 27, 2021 at 13:34
 * 4
 $\begingroup$ Even functions that saturate at 1 (sigmoids, arctan) have similar issues; maybe they are simply less apparent because they don't result in an inf , NaN , or error message. $\endgroup$
 – Federico Poloni
 Jun 27, 2021 at 14:17
 * $\begingroup$ that's a noteworthy point, but if this were the only problem, then defining the activation as g(x)= e^(x-1) would solve it since g(g(g(1))) is just 1. $\endgroup$
 – KishKash
 Jul 4, 2021 at 15:45
 * 1
 $\begingroup$ The input $1$ is just an example, it could be $2$ as well. The point is with a very small input, the network can output very large numbers in a handful of layers. $\endgroup$
 – gunes
 Jul 4, 2021 at 15:47
 Add a comment |

 Your Answer

 Thanks for contributing an answer to Cross Validated!

 * Please be sure to answer the question. Provide details and share your research!

 But avoid …

 * Asking for help, clarification, or responding to other answers.
 * Making statements based on opinion; back them up with references or personal experience.

 Use MathJax to format equations. MathJax reference.

 To learn more, see our tips on writing great answers.

 Sign up or log in

 Sign up using Google
 Sign up using Facebook
 Sign up using Email and Password
 Submit

 Post as a guest

 Name
 Email

 Required, but never shown

 Post as a guest

 Name
 Email

 Required, but never shown

 Post Your Answer Discard

 By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

 Not the answer you're looking for? Browse other questions tagged
 * machine-learning
 * neural-networks
 * perceptron
 or ask your own question.

 * Featured on Meta
 * Accessibility Update: Colors

 Related

 51
 How does rectilinear activation function solve the vanishing gradient problem in neural networks?
 5
 What kind of deep neural networks are (not) data-intensive?
 21
 Why is step function not used in activation functions in machine learning?
 9
 Does the universal approximation theorem for neural networks hold for any activation function?
 4
 Why are nonlinear activations used in convolutional neural networks?
 1
 Can a neural network have an activation function that is a transformation of the parent function?
 2
 Can I use *any* activation function in my Neural network as long as it is non-linear?
 6
 Why aren't neural networks used with RBF activation functions (or other non-monotonic ones)?
 3
 Neural networks - what is the point of having sigmoid activation function $\sigma(.)$ AND sigmoid g(t)?

 Hot Network Questions

 * Sum numbers which are strings
 * Agile Principles: What does it mean sustainable development?
 * Why is mind interacting with matter any more problematic than matter interacting with matter?
 * Can one connect two technic pins end-to-end?
 * A Murder Mystery
 * Does 'peu delicate' mean 'snobbish'
 * Recombination of protons in solar wind
 * What does "invalidate" mean in this context?
 * DisplayPort downstream limited to 1920×1080
 * Hangman game, with secret word provided by user
 * The conjunction "autem" expressing indignation or surprise in questions and exclamations
 * Computers lose power, lights dim. I suspect the mains power supply; is there anything I can do?
 * How does putting an LED in parallel with an 8-bit computer load affect the voltage/current of a circuit?
 * A SF humor short story about the Navy finding a plug at the bottom of the sea, then a giant bar of soap
 * How to play a measure with no chord notation above it?
 * Black-and-white SF movie with giant octopus- or starfish-like creature attacking a city
 * Are there other countries besides Türkiye that insist their name be written in English using diacritics (or even other letters) not common in English?
 * Addressing people in the Japanese Christian/Catholic Church
 * How to get the cd shell-builtin to stop guessing?
 * When did Kievan Rus' rulers become culturally Slavic?
 * We don't own the land so is our house nearly worthless?
 * Searching text about school girl growing wings and trying to hide them at first
 * orbital speed to maintain relative position to the Earth
 * Which country has the highest number of diplomatic missions in a single other country?
 more hot questions
 Question feed

 Subscribe to RSS

 Question feed

 To subscribe to this RSS feed, copy and paste this URL into your RSS reader.

 Cross Validated

 * Tour
 * Help
 * Chat
 * Contact
 * Feedback

 Company

 * Stack Overflow
 * Teams
 * Advertising
 * Collectives
 * Talent
 * About
 * Press
 * Legal
 * Privacy Policy
 * Terms of Service
 * Cookie Settings
 * Cookie Policy

 Stack Exchange Network

 * Technology
 * Culture & recreation
 * Life & arts
 * Science
 * Professional
 * Business
 * API
 * Data
 * Blog
 * Facebook
 * Twitter
 * LinkedIn
 * Instagram

 Site design / logo © 2023 Stack Exchange Inc; user contributions licensed under CC BY-SA. rev 2023.2.7.43225

 Your privacy

 By clicking “Accept all cookies”, you agree Stack Exchange can store cookies on your device and disclose information in accordance with our Cookie Policy.

 Accept all cookies Customize settings
$\begingroup$ 
 
The ReLU function is commonly used as an activation function in machine learning, as well, as its modifications (ELU, leaky ReLU).
The overall idea of these functions is the same: before x = 0 the value of the function is small (its limit to infinity is zero or -1), after x = 0 the function grows proportionally to x.
The exponent function (e^x or e^x-1) has similar behavior, and its derivative in x = 0 is greater than for sigmoid.
The visualization below illustrates the exponent in comparison with ReLU and sigmoid activation functions.
   
So, why the simple function y=e^x is not used as an activation function in neural networks?
 
 
 
 
 asked Jun 26, 2021 at 16:18
 
 
    
 
 
 $\endgroup$
 
 3
 
 
 
 
 
 
 $\begingroup$ 
I think the most prominent reason is stability. Think about having consequent layers with exponential activation, and what happens to the output when you input a small number to the NN (e.g. $x=1$), the forward calculation will look like:
$$o=\exp(\exp(\exp(\exp(1))))\approx e^{3814279}$$
It can go crazy very quickly and I don't think you can train deep networks with this activation function unless you add other mechanisms like clipping.
 
 
 
 answered Jun 26, 2021 at 16:47
 
 
    
 
 
  gunesgunes 54.3k4 gold badges45 silver badges80 bronze badges
 
 
 $\endgroup$
 
 4
Why pure exponent is not used as activation function for neural networks?
https://stats.stackexchange.com/questions/532309/why-pure-exponent-is-not-used-as-activation-function-for-neural-networks