{'fraction_non_alphanumeric': 0.04898767689465364, 'fraction_numerical': 0.018913204959716587, 'mean_word_length': 4.094013777044727, 'pattern_counts': {'":': 0, '<': 20, '<?xml version=': 0, '>': 20, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Deep Neural Networks Are Effective At Learning High-Dimensional Hilbert-Valued Functions From Limited Data

The accurate approximation of scalar-valued functions from sample points is a key task in mathematical modeling and computational science. Recently, machine learning techniques based on Deep Neural Networks (DNNs) have begun to emerge as promising tools for function approximation in scientific computing problems, with impressive results achieved on problems where the dimension of the underlying data or problem domain is large. In this work, we broaden this perspective by focusing on approximation of functions that are Hilbert-valued, i.e. they take values in a separable, but typically infinite-dimensional, Hilbert space. This problem arises in many science and engineering problems, in particular those involving the solution of parametric Partial Differential Equations (PDEs). Such problems are challenging for three reasons. First, pointwise samples are expensive to acquire. Second, the domain of the function is usually high dimensional, and third, the range lies in a Hilbert space. Our contributions are twofold. First, we present a novel result on DNN training for holomorphic functions with so-called hidden anisotropy. This result introduces a DNN training procedure and a full theoretical analysis with explicit guarantees on the error and sample complexity. This error bound is explicit in the three key errors occurred in the approximation procedure: best approximation error, measurement error and physical discretization error. Our result shows that there is a procedure for learning Hilbert-valued functions via DNNs that performs as well as current best-in-class schemes. Second, we provide preliminary numerical results illustrating the practical performance of DNNs on Hilbert-valued functions arising as solutions to parametric PDEs. We consider different parameters, modify the DNN architecture to achieve better and competitive results and compare these to current best-in-class schemes.


Introduction
Driven by their success in many historically-challenging machine learning problems, Deep Neural Networks (DNNs) and Deep Learning (DL) are beginning to be applied successfully to challenging tasks in computational science and engineering. Such tasks are often characterized by the high dimensionality of their data or problem. Examples include inverse problems in imaging Adcock and Hansen (2021); Ongie et al. (2020), molecular dynamics simulations Faber et al. (2017), protein structure prediction Jumper et al. (2020), discovery of unknown dynamical systems Lagergren et al. (2020), Partial Differential Equations (PDEs) Berg and Nyström (2018) and, as discussed below, parameterized PDEs for Uncertainty Quantification (UQ).
The above list represents only a selection of the many recent results in this area. However, it is notable that these results generally fall into the category of existence theory: namely, they assert the existence of a DNN with desirable approximation properties, but not a constructive means to compute such a network (we note in passing several exceptions Dereventsov et al. (2019); Fokina and Oseledets (2019); Daws and Webster (2019a), although these generally lack theoretical guarantees on trainability). As has discussed in Geist et al. (2020) and shown in Adcock and Dexter (2020), there can often be a substantial gap between theoretical existence results and practical performance of DNNs when trained using standard tools. Motivated by such performance gap, this paper seeks to address the following three key issues: (1) The vast majority of previous work considers only scalar-valued function approximation.
(2) Existence theory says little about whether such a DNN can be obtained by training and the number of samples of the function needed to do so.
(3) Many applications in computational science are relatively data starved. Hence it is critical to understand the sample complexity DNN approximation: namely, how much data is needed to train an accurate DNN for a given function.
Specifically, the focus on this paper is on learning high-dimensional, Hilbert-valued functions from limited datasets using DNNs. We next describe the motivations for considering this problem.

Motivations
An important task in UQ involves constructing a surrogate model of a physical system that depends on a set of parameters y ∈ R d . The physical system is typically modelled via a PDE (or system of PDEs) in terms of the relevant spatial and temporal variables. In other words, its solution is a function u = u(x, y), where x ∈ R n , n = 1, 2, 3, 4, denotes the physical variables (space and time) and y denotes the parameters. The function u is the solution to a PDE system where D x (·, y) is a differential operator in the physical variable x that depends on the parametric variable y. The objective in surrogate model construction is to understand the parametric dependence of the solution u. Since the PDE problem can often be posed (in weak form) in a separable Hilbert space V, this is equivalent to approximating the Hilbert-valued function (we suppress the x dependence for ease of notation). This problem have several key features that motivate this work (see, e.g. Cohen and DeVore (2015); Gunzburger et al. (2014) for more detail): (i) The input dimension d is typically high. Indeed, the more parameters, the better the model for the physical system. Typically, d may range from ten to over a hundred. Moreover, in some situations, one may also consider functions with infinitely-many parameters, i.e. d = ∞.
(ii) The output u(y) takes values in a Hilbert space V. In some applications one may only wish to approximate some scalar Quantity-of-Interest (QoI) of u (e.g. the spatial mean f (y) = Ω u(x, y) dx over the physical domain Ω ⊆ R n ). However, other applications call for approximating the whole solution u(y) (from which one can obviously approximate any number of QoI's).
(iii) Computing samples is expensive. For each value of y, evaluating u(y) requires solving the PDE (1) via a computationally-intensive numerical simulation. In practice, generating the samples may take a time ranging from minutes to days or even weeks. Furthermore, since the PDE is never solved exactly, this process always commits an error -the measurement error as we term it.
To summarize, surrogate model construction involves approximating a high-dimensional, Hilbertvalued function from limited data. Such a task is impossible without further assumptions. Fortunately, u is often smooth. To be precise, under certain conditions on the problem (1) one can show that u is a holomorphic function of the parameters y Cohen et al. (2010Cohen et al. ( , 2011Chkifa et al. (2015) (see also Cohen and DeVore (2015) for an in-depth review). This opens the door for approximating it efficiently from limited samples, even when d is large, and thereby lessening the curse of dimensionality. We shall exploit this assumption throughout.
Even with this assumption, though, there is another hurdle to overcome; namely the infinite dimensionality of V. The usual way to address this is to introduce a finite-dimensional discretization V h , where h is a discretization parameter, and compute an approximation taking values in V h . Typically, since u is the solution of a PDE, one takes V h to be a finite element discretization. However, regardless of how V h is chosen, it is important that the discretization error (the effect of replacing V by V h ) be quantified explicitly in the overall error bound. We also address this matter.

Contributions
In this work we study the approximation of a holomorphic, Hilbert-valued function f : U → V (we now switch notation from u to f since the problem we consider does not necessarily need to arise as the solution of a parametric PDE) from m noisy sample values We assume throughout that U = [−1, 1] d is the unit hypercube and the y i are drawn identically and independently from the uniform measure on U. Note that this choice of sampling is not only typical in practice, it is critical in allowing for theoretical sample complexity estimates that scale efficiently with the dimension d. The choice of the hypercube implies that the parameters y 1 , . . . , y d (where ) are independent and vary in finite intervals (which, up to rescaling, can be taken to be equal to [−1, 1]). Both assumptions are also standard in practice.
The values n i are the measurement errors. They constitute the errors involved in computing f , whether it be via a numerical PDE solve or some other unspecified process. Throughout, we assume that V is a separable Hilbert space and that it is discretized via a finite-dimensional subspace V h ⊆ V. We make the additional assumption that the samples d i are elements of V h , i.e. d i ∈ V h , ∀i = 1, . . . , m, and we seek to compute an approximationf : Our first main contribution is a novel theoretical result, Theorem 5. It shows that there exists a DNN architecture (of a given size and depth depending on m and d) and a training procedure (i.e. a loss function) such that any minimizer of the corresponding loss function minimization problem approximates f up to an explicit error bound. This error bound splits into three key terms that fully describe the effects of all main errors involved in the approximation process: (a) an approximation error that is exponentially-small in m, up to log factors; (b) a measurement error that is proportional to the 2 -norm of the measurement noise (n i ) m i=1 ; (c) a physical discretization error that is proportional to the best approximation error of f in V h . We term this result a practical existence theorem for Hilbert-valued function approximation. Going beyond standard existence theory, as discussed above, it shows not only the existence of a DNN with desirable approximation properties, but both a means of obtaining it via training and an estimate on the sample complexity via the exponentially-decaying approximation error (a).
A key facet of this work is the assumption on f . While we assume f is holomorphic, we do not assume any further information on it. In particular, f may have anisotropic dependence on the parameters y 1 , . . . , y d -i.e. it may vary more rapidly in some directions than in others -and such anisotropy may be hidden -i.e. it is not used to construct the approximation. This hidden anisotropy assumption is highly relevant in practice, yet substantially harder to tackle than scenarios where such behaviour is known a priori. We formalize exactly what we mean by hidden anisotropy in §3.
As we note below, current best-in-class methods for holomorphic, Hilbert-valued function approximation are based on multivariate orthogonal polynomials and also achieve exponential rates of convergence in m. Through (a) we show that the DNN training procedure achieves the same rate of convergence, up to a constant. Hence this work shows that DNNs can already perform as well as current best-in-class methods -which, being based on polynomials, are strongly tied to the underlying smoothness -thus opening the door for better performance across a variety of different function classes, including nonsmooth functions.
Further, we also categorize all the errors in the DNN training procedure, including the measurement error (b) and physical discretization error (c). The latter is ubiquitous in simulations (due to the need to work with V h instead of V), but often not included in theoretical analyses. Our work thus sheds light on understanding how to tune the method parameters (m, h, and so forth) to balance the various errors, thus leading to optimal practical performance. It also raises the potential for the use of multilevel or multifidelity schemes, where m and h are simultaneously refined.
In our second contribution, we present initial numerical experiments showing the effectiveness of trained DNNs for parametric PDE problems. We also compare with best-in-class schemes based on polynomial approximation, with a focus on the sample complexity of both approaches. Theorem 5 shows that theres exists a DNN training procedure that performs at least as well such schemes in this regard. However, as we discuss, this procedure deviates from standard DNN architectures and training in several senses, and, perhaps more importantly, it is not generally expected to yield any better performance. Following Adcock and Dexter (2020), in our experiments we use standard DNN architectures coupled with standard loss functions and training algorithms. We show that with proper architecture and hyperparameter selection we are able to achieve competitive results. In particular, we can use smaller architectures than those suggested by Theorem 5 and the standard 2 -loss functions. Further, we show that we can actually outperform state-of-the-art polynomial methods of the problem considered by switching from the ReLU (as used in Theorem 5) to smoother activations functions. All in all, these preliminary results demonstrate the promise of the DNN approach for parametric PDEs in terms of sample complexity. This complements recent results shown in Geist et al. (2020), which showed favourable performance of DNNs with respect to dimension d.

Related work
There are various ways to approximate the solution map of a parametric PDE, including reduced basis methods Hesthaven et al. (2015) and polynomial chaos expansions Xiu and Karniadakis (2002). In this paper, we focus on comparing DNN performance against the latter. Polynomial expansions are well-suited to smooth function approximation, with the so-called best s-term polynomial approximation offering exponential rates of convergence for holomorphic functions (see §3). For a function with isotropic or known anisotropic behaviour in its variables, such an approximation can be computed in a number of ways, including interpolation or least squares Cohen and Migliorati (2018). The situation becomes more challenging when the anisotropy is hidden. Adaptive interpolation or least-squares schemes Chkifa et al. (2013Chkifa et al. ( , 2014; Cohen and Migliorati (2018); Gittelson (2013); Migliorati (2015) may work in practice, though they generally lack theoretical guarantees.
In the last five years, techniques based on compressed sensing have emerged as viable tools for this problem (see, e.g., Adcock et al. (2017Adcock et al. ( , 2019; Chkifa et al. (2018); Doostan and Owhadi (2011); Hampton and Doostan (2015) and references therein). Theoretical guarantees show that such techniques provide quasi-best s-term polynomial approximations, with favourable sample complexities Adcock (2018); Adcock et al. (2017); Chkifa et al. (2018). However, standard compressed sensing only allows for recovery of real or complex sparse vectors, e.g., scalar QoI's of the solution map (2), and therefore does not allow for recovery guarantees for Hilbert-valued functions. The Simultaneous Compressed Sensing (SCS) method Dexter et al. (2019) enables fully discrete approximation of the solution map by combining spatial discretization, e.g., finite elements, with jointsparse vector recovery techniques modified for recovery in V h . The SCS framework also extends theoretical recovery guarantees from compressed sensing to the Hilbert-valued setting, thereby inheriting the same quasi-best s-term approximation rates and sample complexity estimates. Because of its favourable behaviour and theoretical guarantees, we refer to SCS as the current best-in-class method, and seek to match (or beat) this performance with a DNN procedure.
Recently  Laakmann and Petersen (2020) and references therein. These works generally lack theoretical analysis. On the theoretical side, Opschoor et al. (2019) provides an existence theorem for parametric PDEs with holomorphic solution maps. The work  considers the reduced basis approach to parametric PDEs, showing the existence of a DNN whose size depends on the intrinsic low-dimensionality of the solution manifold. Neither result addresses whether such a DNN can be trained, nor the sample complexity in doing so. Our work in particular complements Opschoor et al. (2019) by showing that DNNs admitting the same approximation error can be efficiently trained. Furthermore, our result also makes all errors committed by the process explicit, including the aforementioned measurement and physical discretization errors, which are typically not addressed in previous works.

Learning Hilbert-valued functions via DNNs
We first require some notation. Throughout d ∈ N denotes the dimension of the input space. We write N d 0 := {ν = (ν k ) d k=1 : ν k ∈ N 0 } for the set of nonnegative integer multi-indices and 0 and 1 for the multi-indices consisting of all zeros and all ones respectively. The inequality µ ≤ ν between two multi-indices is understood componentwise.
For 1 ≤ p ≤ ∞, we write · p for the usual vector p -norm. Moreover, for 1 ≤ p, q < ∞ we define the matrix p,q -norm as G q p,q := n j=1 ( m i=1 |G ij | p ) q/p . We also use the notation A B to mean that there exists a numerical constant c > 0 independent of A and B such that A ≤ cB, and likewise for A B.

Setup
We now describe the setup in detail. We write y = (y 1 , . . . , y d ) for the variables and let U = [−1, 1] d . We consider the uniform probability measure on U, i.e.
and write L p (U) for the corresponding weighted Lebesgue spaces of scalar-valued functions over U and · L p (U ) for their norms. Throughout, we let V be a separable Hilbert space over the field R (we could also consider the field C with few additional difficulties), with inner product ·, · V and corresponding norm We let V N be the vector space of Hilbert-valued vectors of length When p = 2 this is a Hilbert space with inner product Next, we define the weighted (Lebesgue-)Bochner space L p (U; V) as the space consisting of (equivalence classes of) strongly -measurable functions f : In general, we cannot work directly in the space V, since it is usually infinite dimensional. Hence, we consider a finite-dimensional discretization Here h > 0 denotes a discretization parameter, e.g. the mesh size in the case of a finite element discretization. In the context of finite elements, assuming (5) corresponds to considering so-called conforming discretizations. We let

DNNs
Let f : U → V and write its projection in terms of the basis Notice that the coefficients c k are scalar-valued functions of y, i.e. c k : U → R.
Our objective is to approximate the coefficients with a DNN. We consider standard feedforward DNN architectures of the form Φ : Here A l : R N l → R N l+1 , l = 0, . . . , L + 1 are affine maps and σ is the activation function, which we assume acts componentwise, i.e. σ(y) : The values {N l } L l=1 are the widths of the hidden layers. By definition N 0 = d and N L+2 = K. Given (6), we write for the resulting approximation to f . We write N for a class of DNNs of the above form. We term L the depth, and denote this as depth (N ). We write parameters (N ) for the number of trainable parameters in N (i.e. the number of weights and biases that parameterize N ). We also write size (N ) for the size of the DNNs in N . This is equal to the total number of nonzero weights and biases.

Problem statement
Let f : U → V be a Hilbert-valued function and consider m sample points y 1 , . . . , y m drawn independently from the uniform measure (3). We consider noisy evaluations of f of the form where n i ∈ V is the ith noise term. What is important is that the samples d i are elements of the finite-dimensional space V h . Hence, the term n i encompasses the error involved in approximating f (y i ) ∈ V by an element of V h . This takes into account that fact that these measurements are computed by some unspecified routine (e.g. a PDE solver in the case of parametric PDEs) which returns a value in the finite-dimensional Hilbert space V h (e.g. a finite element space). We do not specify precisely how this approximation is performed, nor how large this error is, but merely aim to show an error bound that scales linearly with respect to the n i . A particular case is when , but in what follows it is not necessary to assume this, as the numerical computation that leads to d i may not involve computing the projection P h . We remark in passing that one can easily extend this analysis to the case where the space V h used for constructing the approximation f Φ to f differs from that used to compute the evaluation of f . Since any algorithm for learning a DNN needs to take finite inputs, we assume the evaluations (8) are provided to us via the basis {ϕ k }. To be precise, we assume we have access to the data With this in hand, the problem is as follows: Problem 1 Use the data (9) to learn a DNNΦ : R d → R K from the class N and therefore an approximation fΦ to f of the form (7).
3. Holomorphy, best s-term polynomial approximation and hidden anisotropy

Holomorphy
We start by recalling the definition of holomorphy and of holomorphic extension for Hilbert-valued functions. The definition of holomorphy employed here is based on the notion of Gateaux partial derivative, although other equivalent definitions are possible (see (Hervé, 2011, Chapter 2)). The holomorphic extension assumption has played a crucial role in the context of parametric PDEs since the seminal work Cohen et al. (2011); see also Cohen and DeVore (2015) and references therein. : Definition 2 (Holomorphic extension) Let V be a Hilbert space. The function f : We are interested in approximating Hilbert-valued functions f : U → V that admit a holomorphic extension to suitable open neighborhoods of U. Specifically, regions defined by Bernstein (poly)ellipses. When d = 1 the Bernstein ellipse of parameter ρ > 1 is defined by E ρ = 1 2 (z + z −1 ) : z ∈ C, |z| ≤ ρ ⊂ C. This is an ellipse with ±1 as its foci and major and minor semi-axis lengths given by 1 2 (ρ ± ρ −1 ). For d ≥ 1, given ρ = (ρ j ) d j=1 ∈ R d with ρ > 1, we define the Bernstein polyellipse as Given this, we denote the class of unit-norm, Hilbert-valued functions that are holomorphic in E ρ as Note that the parameter ρ dictates the smoothness of f in the different coordinate directions. The larger ρ j , the smoother f is in the jth variable y j . We say that f has anisotropic dependence on the variables y if the ρ j are not all equal and that the anisotropy is hidden if the ρ j 's are unknown.

Polynomial approximation of holomorphic functions
The approximation theory of functions that are holomorphic in Bernstein ellipses is a classical topic, especially in d = 1 dimensions. However, polynomial approximations have proved to be an extremely effective tool in d 1 dimensions as well. Let {Ψ ν } ν∈N 0 be the univariate orthonormal Legendre polynomial basis of L 2 ([−1, 1]). Let f ∈ B(ρ) for some ρ > 1 and consider its expansion where the coefficients {c ν } ν∈N 0 are elements of V. It is well-known that the truncated expansion One aims to construct a polynomial approximation by selecting s terms from this expansion, i.e.
where S ⊂ N d 0 is a multi-index set with |S| = s. Unfortunately, standard, isotropic choices for S such as the (isotropic) total degree set S = {ν = (ν k ) d k=1 : ν 1 + · · · + ν d ≤ n} generally lead to less favourable convergence rates in terms of s. In particular, when the dependence of f on its variables y is anisotropic, the index set S may include many indices that correspond to smallnorm coefficients c ν . This motivates an alternative approach based on best s-term approximation, in which a polynomial approximation is formed by selecting the indices in (10) corresponding to the largest s of the coefficient norms c ν V ; see DeVore (1998); Cohen et al. (2010).
Best s-term approximation is particularly well suited to approximating holomorphic functions with anisotropic dependence. If f ∈ B(ρ) then it is known that, for any > 0, where f s is the best s-term polynomial approximation to f (see Theorem 10). This shows that the best s-term approximation error converges exponentially-fast in s without any prior knowledge of the ρ j 's. In particular, as the right-hand side of (11) depends on the product of the logarithms of the ρ j 's, it is completely independent of their ordering. We note that, although a closed formula for s(d, , ρ) in (11) is not available, an inspection of the proof of Theorem 10 (see Appendix B.3) and of (Opschoor et al., 2019, Theorem 3.5) reveals a constructive characterization of this quantity.
Note that estimates such as (11) says nothing about how to actually construct the approximation f s from sample values. Nevertheless, we consider the rate (11) as the benchmark against which we compare the developed DNN procedure.
Remark 3 There are various results on the convergence rate of best s-term polynomial approximation of a holomorphic function. Algebraic rates of convergence can be found in, for instance, Cohen et al. (2011); Cohen and DeVore (2015). These are attractive since they also hold when d = ∞, thus theoretically permitting the approximation of functions of infinitely-many variables. However, in finite dimensions the constants in these error bounds may be large Tran et al. (2017). The bound (11) is based on (Cohen and DeVore, 2015, Sec. 3.9) and Opschoor et al. (2019). However, the scaling 1/(d + 1) is not sharp. For quasi-optimal error bounds, see Beck et al. (2014Beck et al. ( , 2012; Tran et al. (2017). The results of Tran et al. (2017) are asymptotically sharp as s → ∞. The challenge, however, for our purposes is that they do not generally lead to an approximation in a so-called lower set, unlike (11). This is an important component of subsequent analysis (See §B.3).

Hidden anisotropy
Motivated by this discussion, we now define the following class of functions: Definition 4 Let d ≥ 1, γ > 0 and > 0. Then HA = HA(γ, , d) is the set of Hilbert-valued functions f : U → V that have a holomorphic extension to a Bernstein polyellipse E ρ and satisfy f L ∞ (Eρ;V) ≤ 1, and where the parameters ρ = (ρ j ) d j=1 satisfy Observe that, thanks to (11), for this class of functions we have It is this rate we seek to obtain with the DNN approximation, with some polynomial scaling between m and s. Specifically, we consider the following problem: Problem 2 Devise a DNN architecture and training procedure that solves Problem 1 and for which the error decays exponentially fast in m for all f ∈ HA(γ, , d).

Main result
In this section we present our main result that resolves Problem 2. Its proof can be found in Appendix B.4. Given m ≥ 1, 0 < ε < 1 and d ≥ 1, we first define L = L(m, d, ε) as where c 0 > 0 is a universal constant, and m = m(m, d, ε) := m/L.
Theorem 5 There are universal constants c 0 , c 1 , c 2 , c 3 > 0 such that the following holds. Let d ≥ 1, 0 < , ε < 1, γ > 0, m ≥ 1 and m ≥ 2 d be as in (15). Let y 1 , . . . , y m be drawn independently from the uniform measure (3) on U. Then there is (a) a class of neural networks N with N 0 = d, N L+2 = K and were ∆ N is the number of trainable parameters in N , which satisfies (with the convention that 0! = 1), (b) a regularization functional J : N → [0, ∞) equal to a certain norm of the parameters, and (c) a choice of regularization parameter λ involving only m and d, such that the following holds with probability at least 1 − ε. For all f ∈ HA(γ, , d) with noisy evaluations d i = f (y i ) + n i ∈ V h as in (8), every minimizerΦ of the training problem satisfies for m ≥ 2 d−2s2 , wheres =s(d, , ρ) is as in (13), f Φ is as in (7) and, for e = 1 √ m (n i ) m i=1 , This result asserts that there is a DNN architecture and training procedure with an explicit error bound comprising three terms. First, the approximation error E 1 . This quantifies both how well f is approximated by a DNN in terms of the number of samples m. It is exponentially small in m, up to log terms. Furthermore, the log term L = L(m, d, ) effectively behaves roughly like log 3 (m) log(d), i.e. polylogarithmic in m but only logarithmic in d. The second term is the measurement error E 2 = e V,2 . This is, as discussed, the error in the pointwise evaluations of f at the points y i . Note that we do not assume e stems from any random distribution -the term E 2 implies that the noise can be adversarial, which is generally suitable in problems such as parametric PDEs, where the noise stems from the (deterministic) error incurred in the numerical PDE solve. The third term is the physical discretization error E 3 . Recall that that we cannot work directly in the infinite-dimensional space V. This term accounts for the error induced by working in V h instead. Since the projection P h (f )(y) is the best approximation to f (y) ∈ V from V h in the V-norm, the term E 3 is optimal, up to a minor switch from the L 2 (U; V)-norm on the left-hand side of (17) to the L ∞ (U; V)-norm on its right-hand side.
It is notable that the approximation error E 1 achieves a similar exponential rate of decay as the best s-term polynomial approximation error bound (13), except with s replaced by √ m and an additional factor of 1/ √ 2. Note that the former translates into m ∝ s 2 , which is precisely the quadratic sample complexity that arises when computing s-sparse polynomial approximations from uniformly-distributed samples via compressed sensing Chkifa et al. (2018); Adcock et al. (2017). Hence, the presence of the term E 1 implies that the DNN procedure performs as well as current best-in-class polynomial approximation schemes based on compressed sensing.
Note that the training problem may at first appear somewhat unusual, since it involves fitting in the continuous · V -norm. However, it can be reformulated as a discrete norm.
Gc, c , which involves a standard (weighted) 2 -norm fit over R K . However, we also note that the loss function in (16) is is rather different to the standard 2 -loss used commonly in DL, since the data fidelity term is not squared. The reason for this, as discussed in §B.2, is to ensure property (c). If the loss term were squared, the optimal choice regularization parameter λ would depend on the unknown error terms E 1 , E 2 and E 3 . It is rather remarkable that forgoing the squaring leads to such a desirable property for λ. We are not aware of any existing results in DNN training where the regularization parameter can be set optimally in a function-independent way.

Numerical exploration
We now present preliminary experiments on the efficiency of fully-connected DNNs for approximating solutions to parametric PDEs from limited data. A more comprehensive study of the efficacy of DL techniques for problems in computational UQ will presented in an upcoming work. Differing from previous works -including Geist et al. (2020), which focused on the scaling of accuracy with the dimension d in the presence of large amounts of training data -in this study we investigate the trade-off between accuracy and the number of samples m of the solution of the PDE. We also include a comparison to the best-in-class SCS scheme Dexter et al. (2019) (see §1.3 and §A.4).

Setup
We approximate a parametric elliptic PDE defined over the spatial domain Ω = (0, 1) 2 and parametric domain U = [−1, 1] d with the uniform probability measure (3). Specifically, given g ∈ L 2 (Ω) we seek a function u : Ω × U → R satisfying − ∇ · (a(x, y)∇u(x, y)) = g(x), ∀x ∈ Ω, u(x, y) = 0, ∀x ∈ ∂Ω We choose the function g to be independent of y for simplicity. We use a finite element method for spatial discretization, based on a triangulation with K = 1089 nodes and meshshze h = √ 2/32. In this study we consider only the error in approximating the parametric PDE, therefore we use the same finite element discretization in generating all of the sample training and testing data, and also in discretization and solution of the SCS problem. See §A.1 for further details.
The training data for the DNNs and for the SCS method is generated by solving (19) at a set of uniform random points {y i } mmax i=1 ⊂ U, yielding a set of solutions {u h (·, y i )} mmax i=1 , where u h (·, y) ∈ V h is the computed solution of (19). Here we examine both the efficiency of DL techniques on fixed subsets of training data of size m ≤ m max and the effect of increasing m up to m max on the accuracy of the approximations. We also examine the average performance of the DNNs and SCS method by running multiple trials on different sets of uniform random points generated with the trial number as the seed. The testing data is generated in the same way by evaluating (19) at a set of points {y i } mtest i=1 . However, instead of using random points to test, we use a high-order stochastic collocation reference solution (see Nobile et al. (2008)) based on a deterministically-generated sparse grid quadrature rule. The testing error is recorded in the Bochner norms L 2 (U; L 2 (Ω)) and L 2 (U; H 1 0 (Ω)) (see (4)), approximated by this rule. See §A.2 for further details. The DNN strategy is based on §2.2. For y ∈ U, we define the Finite Element (FE) and DNN approximations by respectively. We consider several types of fully-connected DNN architectures with input dimension N 0 = d, output dimension N L+2 = K, and constant hidden-layer widths, i.e. N 1 = N 2 = . . . = N L = M for some M ∈ N. We follow the convention from Adcock and Dexter (2020) of referring to a DNN with activation function σ, L hidden layers and M nodes per hidden layer as a σ L × M DNN. In this study we focus on either the tanh, ρ(x) = tanh(x), or Leaky-ReLU, ρ(x) = max{x, 0.2x}, activation function. 2 For details on the training approach, see §A.3.

Effective architectures and loss functions and efficiency of training
We first study the performance of the DL approach based on two different loss functions. The first loss function is the well-known mean square error loss, given in terms of the coefficients (20) as On the other hand, motivated by the recent work (Geist et al., 2020) we also consider the squared V-norm loss function, which incorporates information about the complexity of the domain in the training procedure. This is defined by In Figure 1, we study the efficacy of minimizing either (21) or (22) given solution data from problem (19) with a simple affine coefficient given by a(x, y) = 3 + x 1 y 1 + x 2 y 2 .
We make several observations. First, DNNs trained with loss (22) are unable to achieve training errors of the same order as networks trained with (21) given the same budget of epochs. DNNs trained with (22) (22) underperform identical architectures trained with the MSE loss in both metrics, with errors stagnating at order 10 −1 after 600 epochs of training. However, it is interesting that the difference is largest in the case of the L 2 (U; L 2 (Ω)) testing error, where best-performing tanh and Leaky-ReLU 5 × 50 DNNs trained with the MSE loss achieved order 10 −3 error. Figure 2 displays a visualization comparing solutions output by the DNN at a fixed y both early in the training and with a much improved result after achieving order 10 −6 loss error.

Sample complexity of training
In this section we run a large-scale study of the average testing errors and training times over a range of 10 trials, with increasing number of training samples m up to m max = 675, and testing with a sparse grid quadrature rule with m test = 1861 testing points. For these experiments, we consider a modification of the example from Nobile et al. (2008) of a diffusion coefficient with one-dimensional (layered) spatial dependence given by .
Here we let β c = 1/8, and β p = max{1, 2β c }, β = β c /β p . We consider this problem with parameter dimension d = 30. Figure 2 displays the result solving problem (19) with coefficient (24) with SCS and a variety of DNN architectures. For the SCS method, we use the Legendre basis and hyperbolic cross of order p = 6 with cardinality N = 1486. The DNN architectures have depth parameter L and number of hidden layer nodes M chosen so that the ratio L/M = 0.1. In testing we observe competitive performance using DNNs with Leaky-ReLU activation function. We also observe superior performance over SCS and Leaky-ReLU DNNs using DNNs with tanh activation function, with such networks achieving testing errors on average approximately 3.2 times lower in the L 2 (U; L 2 (Ω)) error and 2.15 times lower in the L 2 (U; H 1 0 (Ω)) error than the SCS method given 675 training samples. We also include a comparison of the average training time of the SCS and DNN approaches. While the average times are overall quite similar between the SCS method and training the DNNs, we note that the DNNs are trained on a GPU with accelerated matrix-vector product operations. Therefore this comparison does not provide a good estimate of computational complexity. We leave a study of the computational efficiency of DL techniques for such problems to a future work.
These preliminary results indicate the promise of the DNN approach. The architecture and training differ from that described in Theorem 5; in particular, the networks are smaller and loss function simpler. Theorem 5 showed that learning with DNNs can perform as well as best-in-class schemes. Yet, through a suitable choice of architecture and optimizer, we are not only able to match such performance using a simpler setup, but also outperform it in this example. These preliminary results indicate significant practical promise of the DNN approach, as well as improvements to the theory. Results of a more comprehensive study will be reported in a future work.
M. Stoyanov, D. Lebrun-Grandie, J. Burkardt, and D. Munster. Tasmanian, 9 2013 For the spatial discretization we rely on the finite element method as implemented by Dolfin (Logg and Wells, 2010), and accessed through the python FEniCS project (Alnaes et al., 2015). We generate a regular triangulation T h of Ω composed of triangles T of equal diameter h T = h. We consider a conforming discretization, meaning an arbitrary finite dimensional subspace V h ⊂ V := H 1 0 (Ω), where V h is chosen as the usual Lagrange Finite Elements (FE) of order k = 1. We rely on the Dolfin UnitSquareMesh method to generate a mesh with 33 nodes per side, corresponding to a finite element triangulation with K = 1089 nodes, 2048 elements and meshsize h = √ 2/32. See (Dexter et al., 2019) for further implementation details.

A.2. Testing error
As discussed, rather than using a random points to evaluate the test error, we use a high-order stochastic collocation reference solution (see Nobile et al. (2008)) based on a deterministicallygenerated quadrature rule. Specifically, we apply a Smolyak sparse-grid quadrature rule based on Clenshaw-Curtis points to generate the testing data with the TASMANIAN software package Stoyanov et al. (2013);Stoyanov (2015). In testing, we choose the level of the sparse grid rule such that m test m max . The testing error is recorded in the Bochner norms L 2 (U; L 2 (Ω)) and L 2 (U; H 1 0 (Ω)) (see (4)) and approximated as the square root of the result of the quadrature formulas respectively, where w i are the quadrature weights associated with the sparse grid rule andũ h are the approximations obtained with either the DNNs or SCS. When testing multiple trials, we report the average of the errors over all of the trials.

A.3. DNN training
We follow a similar training methodology to that of Adcock and Dexter (2020). All weights and biases of the DNN are initialized as normal random variables with mean 0 and variance 0.01, with the same seed 0 for each network. We perform calculations in single precision, using the Adam optimizer Kingma and Ba (2014) with an exponentially-decaying learning rate and training for 50,000 epochs or until a stopping tolerance of 5 × 10 −7 is met. Due to the large size of the training data, we use a batch size of min{m, 256} for each training set of size m. As discussed, the loss function is chosen either as the MSE loss (21) or the MVNSE loss (22).

A.4. The SCS method
The SCS method uses multivariate orthonormal Legendre polynomial approximation as described in §3.2 and Appendix B.1, with the hyperbolic cross index set Λ = Λ HC s−1 as in (26). Given the measurement matrix A ∈ R m×N from (27) and measurement vector Following §3.2, we first consider approximating f using the tensor Legendre polynomial basis and let c = (c ν ) ν∈N d 0 ∈ 2 (N d 0 ; V) be the infinite vector of coefficients of f . Fix s ∈ N and let be the hyperbolic cross index set of index s − 1. Let N = |Λ| and ν 1 , . . . , ν N be an indexing of the multi-indices in Λ and define define the normalized measurement matrix.
We also define the the normalized measurement and error vectors as the truncated expansion of f based on the index set Λ and as the finite vector of coefficients of f with indices in Λ. Then we have and therefore where Hence, the recovery of the coefficients c Λ of f is equivalent to solving the noisy linear system of equations (29).

B.2. Hilbert-valued compressed sensing
To solve this, we consider techniques from Hilbert-valued compressed sensing Dexter et al. (2019).
In classical (scalar-valued) compressed sensing, one aims to solve (29) by finding a solution with minimal 1 -norm. In Hilbert-valued compressed sensing, we use the 1 (Λ; V) instead. Both classical and Hilbert-valued compressed sensing are usually formulated as either a quadratically-constrained basis pursuit or LASSO problem. However, we instead consider the Square Root LASSO (SR-LASSO) problem min where λ > 0 is a parameter. This is based on ideas of Adcock et al. (2019). While the other approaches are more common, they have the undesirable feature that the optimal choices for their parameters (in the sense that they give the best error bounds) depend on the magnitude of the terms e and e . These terms are unknown (in particular, they typically depend on f ). Hence, they would not give rise to a result such as Theorem 5, where the parameter choice is independent of f . Fortunately, as shown in Adcock et al. (2019), the SR-LASSO (30) has such a desirable property.
The theory of Hilbert-valued compressed sensing has been developed in Dexter et al. (2019). We recall several key definitions. The support of a Hilbert-valued vector We write Σ s ⊆ V N for the set of s-sparse Hilbert-valued vectors. For nonsparse vectors, we also define the p -norm best s-term approximation error as We now focus on properties of the matrix A that ensure recovery of approximately sparse vectors via (30). To this end, we now define some additional notation. Given a set S ⊆ {1, . . . , N } we write P S : V N → V N for the projection that, given a vector x ∈ V N , produces a vector P S (x) whose ith entry is equal to x i if i ∈ S and zero otherwise.
Definition 6 The matrix A ∈ R m×N satisfies the robust Null Space Property (rNSP) of order 1 ≤ s ≤ N over V N with constants 0 < ρ < 1 and τ > 0 if for any S ⊆ [N ] with |S| ≤ s.
See (Dexter et al., 2019, Defn. 4.1). In (Dexter et al., 2019, Prop. 4.2) the authors show that the rNSP is sufficient to provide an error bound for minimizers of the quadratically-constrained basis pursuit problem. We need an analogous result for the SR-LASSO problem (30). The following is a straightforward extension of (Adcock and Hansen, 2021, Thm. 6.4) from the scalar to the Hilbertvalued case: Lemma 7 Suppose that A ∈ R m×N has the rNSP of order 1 ≤ s ≤ N with constants 0 < ρ < 1 and τ > 0. Let and C 2 = (3ρ+5)τ 2(1−ρ) . Then every minimizerx ∈ V N of the Hilbert-valued SR-LASSO problem min Now we recall that a matrix A ∈ R m×N satisfies the Restricted Isometry Property (RIP) of order 1 ≤ s ≤ N with constant 0 < δ 2s < 1 if δ 2s is the smallest constant δ for which The following result connects the RIP and rNSP. It is based on a combination of (Dexter et al., 2019, Prop 4.3) (which applies to the Hilbert-valued case) and (Adcock and Hansen, 2021, Lem. 5.17) (which gives explicit values for ρ and τ ): Lemma 8 Suppose that A ∈ R m×N satisfies the RIP of order 2s with constant δ 2s < √ 2 − 1. Then A satisfies the rNSP of order s over V N with constants ρ and τ given by With this in mind, we end this section with a result asserting that conditions under which the matrix defined in (27) satisfies the RIP: be the orthonormal tensor Legendre polynomial basis of L 2 (U) and y 1 , . . . , y m be drawn independently from the uniform measure on U. Let 0 < δ, ε < 1 and suppose that where c > 0 is a universal constant. Then, with probability at least 1 − ε, the matrix A defined in satisfies the RIP of order s with constant δ s ≤ 1/4. Note that the choice of 1/4 here is arbitrary. Any value less than √ 2 − 1 ≈ 0.41 will suffice. Proof Theorem 2.2 of Chkifa et al. (2018) implies that any matrix of this form satisfies the RIP of order s with constant δ s ≤ 1/2, provided (after simplifying) m Θ 2 · s · log(2Θ 2 s) · max log(2Θ 2 s log(2Θ 2 s)) · log(2N ), log(2ε −1 log(2Θ 2 s)) , where N = |Λ| and Θ is defined by The Legendre polynomials attain their maxima at y = 1 and, due to the normalization with respect to the uniform measure, satisfy |Ψ ν (1)| = d k=1 √ 2ν k + 1. Hence, since Λ is the hyperbolic cross of index s − 1 (recall (26)), we have Furthermore, it was also show in (Chkifa et al., 2018, Lem. 3.5) gives that Θ 2 ≤ s log(3)/ log(2) . Substituting this into the above expression and using the fact that log(2Θ 2 s) ≤ log(2) + (log(3)/ log(2) + 1) log(s) log(2s) ≤ 2s. now shows that (B.2) is implied by m 2 d · s 2 · log(2s) · log(2s) · log(2N ) + log(ε −1 ) .
Substituting this into (31) now gives the result.

B.3. Polynomial approximation error bounds for holomorphic functions
In order to establish exponential rates of convergence, we need the following result for estimates the error of the best s-term polynomial approximation of a holomorphic function f . For this next result, we recall that a set S ⊆ N d 0 is lower if whenever ν ∈ S then µ ∈ S for all multi-indices µ satisfying µ ≤ ν.
Theorem 10 Let d ∈ N, f : U → V be holomorphic in a Bernstein polyellipse E ρ for some ρ > 1. Then for every > 0 there existss =s(d, , ρ) such that, for every s ≥s, there is a lower set Furthermore, the same bound also applies to the coefficient error c is the sequence of coefficients of f as in (10) and c S is the infinite sequence with νth entry equal to c ν if ν ∈ S and zero otherwise.
Note that this result immediately implies (11). In fact, it is stronger, since it asserts an L ∞ -norm bound, and shows that this can be achieved by a set that is also lower. Proof The proof is mainly based on (Opschoor et al., 2019, Theorem 3.5). We start by proving the theorem in the scalar-valued case, i.e. for V = R. Note that in this case the coefficients c ν ∈ R.
Observing that (35) defines a monotone decreasing relation between τ (s) and s, we can find the minimum values =s(α, d, ρ) ∈ N withs ≥ 2, such that for every s ≥s, we have τ (s) ≤τ . In this way, we have , ∀s ≥s.
Now, we observe that the cardinality of S τ can be explicitly bounded as in (Opschoor et al., 2019, Lemma 3.3) using a volumetric argument (note that S τ corresponds to Λ ε in Opschoor et al. (2019)). This cardinality bound can be written as , ∀τ ∈ (0, 1).
Combining the above inequalities and recalling that |S τ (s) | ≥ s/(d + 1) d , we obtain that, for any s ≥s, Finally, we let α = 1/(1 + ε) 1/d and observe that This concludes the proof in the case V = R. This proof can be generalized to the Hilbert-valued case by replacing coefficients' magnitude |c ν | with coefficients' norm c ν V . With this modification, the analogous of (33) holds. Indeed, the proof of (33) given in (Opschoor et al., 2019, Theorem 3.5) relies on coefficient bounds for Legendre polynomials that hold in the Hilbert-valued case as well (see, e.g., (Cohen and DeVore, 2015, Lemma 3.15). The rest of the argument is identical to the scalar-valued case. We also observe that the exponential bound holds for c − c S 1,V because Ψ ν L ∞ (U ) ≥ 1 for every ν ∈ N d 0 . This completes the proof.

B.4. Proof of Theorem 5
We are now ready to prove the main result. Our strategy is based on reformulating the above compressed sensing formulation as a DNN training problem. We first require the following result (Opschoor et al., 2019, Prop. 2

.10):
Proposition 11 There exists a universal constant c > 0 such that the following holds. For every finite subset Λ ⊂ N d 0 and every 0 < δ < 1 there exists a ReLU neural network Φ Λ,δ : The depth and size of this network satisfy depth(Φ) ≤ c · (1 + d log(d)) · (1 + log(m(Λ))) · (m(Λ) + log(δ −1 )), The general idea of the proof is to use Proposition 11 to approximate the matrix-vector multiplication Az -which is a polynomial evaluated at the sample points y i -as a neural network Φ evaluated at the sample points, or equivalently, to approximate the matrix A, which is built from polynomials, as a matrix A built from DNNs. We then use compressed sensing results applied to A to establish an error bound. Since this process commits an error, we first require the following result, which shows that the rNSP is robust to small matrix perturbations: Lemma 12 Suppose that A has the rNSP of order 1 ≤ s ≤ N over V N with constants 0 < ρ < 1 and τ > 0. Let A ∈ R m×N be such that A − A 2 ≤ δ, where Then A has the rNSP of order s over V N with constants 0 < ρ < 1 and τ > 0 given by This is a straightforward extension of (Adcock and Hansen, 2021, Lemma 8.5) from the scalarvalued case to the Hilbert-valued case. Proof (of Theorem 5) We divide the proof into several steps: Step 1: We first define the class of DNNs N . Let m be as in (15) and define Now let Λ = Λ HC s−1 and Φ Λ,δ be as in Proposition 11. We will choose δ in Step 5. Letting N = |Λ|, we define the class N as Note that Z is the matrix of trainable parameters (it is the weight matrix in the output layer). We will establish the size, depth and number of trainable parameters of this class (part (a) of Theorem 5) in Step 6.
Step 2: We next show that (16) can be reinterpreted as a SR-LASSO minimization problem. Analogously to (27), we define the matrix A as where Φ ν,δ is the νth component of Φ Λ,δ . Let Φ = Z Φ Λ,δ ∈ N and f Φ be as in (7). Then We now define the regularization functional J as We verify that J is a norm over the trainable parameters Z (part (b) of Theorem 5). Notice that where e j ∈ R N is the jth coordinate vector, G 1/2 is the unique positive definite square root of G (defined in §4, along with the norm · G,2 ) and · 2,1 is the matrix 2,1 -norm, defined in §2. Hence, J is a norm over the trainable parameters. Using this and (37), we now see that (16) can be expressed as the Hilbert-valued SR-LASSO problem min We will choose the value of λ in Step 4(ii). In particular, any minimizerĉ = (ĉ ν ) ν∈Λ of (38) yields a minimizerΦ of (16), given bŷ is such that and vice versa. This completes Step 2.
Lemma 9 now implies that the matrix A satisfies the RIP of order 2s with constant δ 2s ≤ 1/4 and Lemma 8 implies that A satisfies the rNSP of order s over V N h with constants ρ = √ 2/3, τ = 2 √ 5/3, with probability at least 1 − ε. To show that A satisfies the rNSP we use Lemma 12. We first bound A − A 2 . Let z ∈ R N . Then we observe that (the choice of δ that will be made in Step 5 will ensure this condition). Using Lemma 12 and the values (40) for ρ and τ we now see that A has the rNSP of order s over V N h with constants with probability at least 1 − ε.
Step 4: We now derive an error bound for f − fΦ, whereΦ is a minimizer of (16) and fΦ is as in (7). The problem (38), or equivalently (16), involves three discretizations. Truncation of the infinite expansion via the index set Λ, discretization of the space V via V h and replacement of the polynomial functions by DNNs. First, givenΦ = C Φ Λ,δ we letĉ be the corresponding minimizer of (38) and set fΨ = ν∈Λĉ ν Ψ ν .
Here P h is as defined in §2.1 and f Λ is as in (28). We now bound A 1 , A 2 and A 3 .
Step 4(i): We commence with A 1 . This is elementary. Since P h is an orthogonal projection we have P h (a) V ≤ a V for a ∈ V and therefore recalling that we are considering a probability measure d (y) over U.
Step 6: Having now shown the error bound, the final step of the proof involves bounding the size, depth and number of trainable parameters for DNNs in the class N defined in Step 1. To do this, we notice from (50) that δ s −1/2 N −1/2 exp(−γs 1/d ), from which, after a few steps, implies log(δ −1 ) log(s) + log(N ) + γs 1/d .