Optimal deep learning of holomorphic operators between Banach spaces

Authors: ben  adcock
nick  dexter
sebastian  moraga
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-11 23:46:25

Abstract: Operator learning problems arise in many key areas of scientific computing where Partial Differential Equations (PDEs) are used to model physical systems. In such scenarios, the operators map between Banach or Hilbert spaces. In this work, we tackle the problem of learning operators between Banach spaces, in contrast to the vast majority of past works considering only Hilbert spaces. We focus on learning holomorphic operators -- an important class of problems with many applications. We combine arbitrary approximate encoders and decoders with standard feedforward Deep Neural Network (DNN) architectures -- specifically, those with constant width exceeding the depth -- under standard $\ell^2$-loss minimization. We first identify a family of  DNNs such that the resulting Deep Learning (DL) procedure achieves optimal generalization bounds for such operators. For standard fully-connected architectures, we then show that there are uncountably many minimizers of the training problem that yield equivalent optimal performance. The DNN architectures we consider are `problem agnostic', with width and depth only depending on the amount of training data $m$ and not on regularity assumptions of the target operator. Next, we show that DL is optimal for this problem: no recovery procedure can surpass these generalization bounds up to log terms. Finally, we present numerical results demonstrating the practical performance on challenging problems including the parametric diffusion, Navier-Stokes-Brinkman and Boussinesq PDEs.

Main paper: 
 introduction: 
 Operator learning is increasingly being investigated for problems arising in computational science and engineering. These problems are often posed in terms of Partial Differential Equations (PDEs), which can be viewed as operators mapping function spaces to function spaces. Depending on the requirements for well-posedness of the PDE, both the input and solution spaces are often Hilbert, or, more generally, Banach spaces. The aim of operator learning is to efficiently capture the dynamic behavior of these operators using surrogate models, typically based on Deep Neural Networks (DNNs). Specifically, we want to learnEQUATIONwhere Y is the PDE solution space, X represents the data supplied to the PDE, i.e., possibly multiple functions describing initial and boundary conditions or forcing terms or, equivalently, a vector of parameters defining such functions.Let µ be a probability measure on X . Given noisy training dataEQUATIONwhere X 1 , . . . , X m ∼ i.i.d. µ and E i is noise, a typical operator learning methodology consists of three objects: an approximate encoder E X : X → R d X , an approximate decoder D Y : R d Y → Y and a DNN N : R d X → R d Y . It then approximates F asF ≈ F := D Y • N • E X . (1.3)The encoder and decoder are either specified by the problem, learned separately from data, or learned concurrently with N . The goal, as in all supervised learning problems, is to ensure good generalization via the learned operator F from as little training data m as possible.As noted in, e.g., [16, 66] , the theory of deep operator learning is still in its infancy. We contribute to this growth in the following ways. We consider learning classes of holomorphic operators (Assumption 2.2), with arbitrary approximate encoders E X and decoders D Y . As we explain in §2.3 (see also [52, §5.2] , [53, §3.4 ] and [41] ) these operators are relevant in many applications, notably those involving parametric PDEs. The main contributions of this work are as follows.1. We consider operators taking values in general Banach spaces. As noted, the vast majority of existing work (with the notable exception of [16] ) considers Hilbert spaces. 2. We consider standard feedforward DNN architectures (constant width, width exceeds depth) and training procedures (ℓ 2 -loss minimization). 3. (Theorem 3.1) We construct a family of DNNs such that any approximate minimizer of the corresponding training problem satisfies a generalization bound that is explicit in the various error sources: namely, an approximation error, which decays algebraically in the amount of training data m; encoding-decoding errors, which depend on the accuracy of the learned encoders and decoders; an optimization error, and; a sampling error, which depends on the noise E i in (1.2). 4. These DNN architectures are problem agnostic; they depend on m only. In particular, the architectures are completely independent on the regularity assumptions of target operator. 5. (Theorem 3.2) We show that training problems based on any family of fully-connected DNNs possess uncountably many minimizers that achieve the same generalization bounds. 6. (Theorems 3.1-3.2) We provide bounds in both the L 2 µ -and L ∞ µ -norms that hold in high probability, rather than just expectation. 7. (Theorems 4.1-4.2) We show that the generalization bound is optimal with respect to m: no learning procedure (not necessarily DL-based) can achieve better rates in m up to log terms. 8. Finally, we present a series of experiments demonstrating the efficacy of DL on challenging problems such as the parametric diffusion, Navier-Stokes-Brinkman and Boussinesq PDEs, the latter two of which involve operators whose codomains are Banach, as opposed to Hilbert, spaces.Approximating an operator between function spaces with training data obtained through numerical PDEs solves presents a formidable challenge. Nevertheless, in recent years, significant advances have been made through the development of DL techniques, leading to the field of operator learning [15, 40, 51, 53, 56, 58, 60, 62, 69, 83, 93, 103] . These approaches often leverage intricate DNN architectures to approximate the complex mappings inherent in physical modelling scenarios. Many works have also focused on the practical aspects of operator learning in real-world applications [13, 24, 35-37, 42, 45, 48, 49, 59, 61, 64, 65, 70, 73, 77, 80, 81, 84, 96, 98-101, 104, 105] .On the theoretical side, universal approximation theorems for operator learning have been developed in [50, 55, 68, 69] and elsewhere. Such bounds are typically not quantitative in the size of the DNN needed to achieve a certain error. For this, one typically either restricts to specific operators (e.g., certain PDEs) or imposes regularity conditions. One such assumption is Lipschitz regularity -see [10, 16, 55, 66, 87] and references therein. However, learning Lipschitz operators suffers from a curse of parametric complexity [54] , meaning that algebraic rates may not be achievable. Another common assumption is holomorphy. While stronger, it is, as noted, very relevant to operator learning problems involving parametric PDEs. Quantitative approximation results for holomorphic operators have been shown in [26, 31, 41, 55, 71] and elsewhere.However, these works do not consider the generalization error, i.e., the error incurred when learning the approximation (1.3) from the finite training data (1.2) . This is particularly important in applications of operator learning where data is obtained through expensive numerical PDE solves, since such problems are highly data-starved. Several works have tackled this question from the perspective of statistical learning theory and nonparametric estimation [16, 55, 66] , but only for Lipschitz operators.As observed in [6, §9.5] , this approach generally leads to a best O(m -1/2 ) decay of the L 2 µ -norm error with respect to m. Theorem 4.1 shows that such a rate is strictly suboptimal for learning the classes of holomorphic operators we consider. Our generalization bounds in Theorems 3.1-3.2 do not use such techniques, and yield near-optimal rates in both the L 2 µ -and L ∞ µ -norms. See also [30] for some related work in this direction for reduced-order modelling with convolutional autoencoders.Our work is inspired by recent research on learning holomorphic, Banach-valued functions [2, 5, 6] . We extend both these works, in particular, [5] , to learning holomorphic operators. We also significantly improve the error decay rates in [5] with respect to m and show they can be achieved using substantially smaller DNNs with standard training (i.e., ℓ 2 -loss minimization). See Remarks C.1-C.2. Our theoretical guarantees fall into the category of encoder-decoder-nets [52] , which includes the well-known PCA-Net [10] and DeepONet [68] frameworks. As in other recent works [16, 30, 55, 66] , in Theorems 3.1-3.2 we assume the encoder-decoder pair (E X , D Y ) in (1.3) have been learned, and focus on the generalization error when training the DNN N . , respectively, for 1 ≤ p < ∞, and analogously for p = ∞ (see, e.g., [8, 44] ). Notice that |||F ||| L p µ (X ;Y) ≤ ∥F ∥ L p µ (X ;Y) for 1 ≤ p < ∞, while |||F ||| L ∞ µ (X ;Y) = ∥F ∥ L ∞ µ (X ;Y) . Throughout this work, ℓ p (N), 0 < p ≤ ∞ denotes the standard ℓ p space with (quasi-)norm ∥ • ∥ p . We also define the monotone ℓ p space ℓ p M (N) as the space of all sequences z = (z i ) ∞ i=1 ∈ R N whose minimal monotone majorant z ∈ ℓ p (N). Here z = (z i ) ∞ i=1 is defined as zi = sup j≥i |z j |. Given a (componentwise) activation function σ, we consider feedforward DNNs of the formEQUATIONwhere A l : R N l → R N l+1 are affine maps, and N 0 = n and N L+2 = k. We define width(N ) = max{N 1 , . . . , N L+1 } and depth(N ) = L. We denote a class of DNNs of the form (2.1) with a fixed architecture (i.e., fixed activation function, depth and widths) as N , and write width(N ) = max{N 1 , . . . , N L+1 } and depth(N ) = L. 
assumptions and setup: 
 Let F : X → Y be the unknown operator we seek to learn andE X : X → R d X , D X : R d X → X , E Y : Y → R d Y , D Y : R d Y → Ybe approximate encoders and decoders for X and Y, respectively. As mentioned, we assume that these maps have already been learned, and focus on the training of the DNN N in (1.3). Our main results allow for arbitrary encoders and decoders (subject to the assumptions detailed below), and provide generalization bounds that are explicit in these terms: specifically, they depend on how well each encoder-decoder pair approximates the respective identity map on X or Y.In order to formulate the precise notion holomorphy for the operator F , we require the following. Let D = [-1, 1] N and ϱ be the uniform probability measure on D. Given ρ > 1, we define the Bernstein ellipse E(ρ) = (x + x -1 )/2 : x ∈ C, 1 ≤ |x| ≤ ρ ⊂ C, and, for convenience, we let E(1) = [-1, 1]. Next, for ρ = (ρ i ) i∈N ≥ 1, we define the Bernstein polyellipse as the productE(ρ) = E(ρ 1 ) × E(ρ 2 ) × • • • ⊂ C N . Definition 2.1 (Holomorphic map). Let ε > 0, b ∈ ℓ 1 (N) with b ≥ 0. A Banach-valued function f : D → Y is (b, ε)-holomorphic if it is holomorphic in the region R(b, ε) = E(ρ) : ρ ≥ 1, ∞ j=1 (ρ j + ρ -1 j )/2 -1 b j ≤ ε ⊂ C N , b = (b j ) j∈N . (2.2)See, e.g., [17, 85] . As noted in [7] we can, by rescaling b, assume that ε = 1. For convenience, we define the following unit ball, consisting of all such functions of norm at most one over R(b, 1):H(b) = {f : D → V (b, 1)-holomorphic : ∥f (x)∥ Y ≤ 1, ∀x ∈ R(b, 1)}. (2.3)Assumption 2.2. Let D = [-1, 1] N and ϱ be the uniform probability measure on D.(A.I) There is a measurable mapping ι : X → R N such that pushforward ς := ι♯µ is a quasi-uniform measure supported on D and ι| supp(µ) : X → ℓ ∞ (N) is Lipschitz with constant L ι ≥ 0.(A.II) The operator F has the form F = f • ι, where f ∈ H(b) for some b ∈ ℓ p M (N) and 0 < p < 1.(A.III) The map E X := ι d X • D X • E X is measurable (here ι d X : X → R d X is the restriction of ι, i.e., ι d X (X) = (ι(X) i ) d X i=1) and the pushforward ς := E X ♯µ is absolutely continuous with respect to ϱ. (A.IV) The maps D Y and E Y are linear and bounded. Now let X 1 , . . . , X m ∼ i.i.d. µ and consider the training data{(X i , Y i )} m i=1 ⊂ (X × Y) m , where Y i = F (X i ) + E i ∈ Y (2.4)and E i ∈ Y represents noise. Let N be a class of DNNs N : R d X → R d Y , and defineF ≈ F := D Y • N • E X , where N ∈ argmin N ∈N 1 m m i=1 ∥Y i -D Y • N • E X (X i )∥ 2 Y .(2.5)We now discuss (A.I)-(A.IV). In §6 we describe future work on relaxing these assumptions.(A.I) is a weak assumption. It asserts that there is a Lipschitz map ι under which the pushforward of µ is a quasi-uniform measure supported in D. As we discuss in Example 2.3, this is notably the case when µ is the law of some random field with an affine parametrization involving bounded random variables -a situation that occurs frequently in parametric and stochastic PDE problems. (A.II) describes the specific holomorphy of the operator F -see Remark 2.4 for details. Note that we require b ∈ ℓ p M (N), not just b ∈ ℓ p (N). It is known [7] that one cannot learn holomorphic functions (and hence operators) from finite data if b ∈ ℓ p (N) only. (A.III) is a relatively weak assumption. In view of (A.I), we expect it to hold as long as the D X • E X ≈ I X sufficiently well. Finally, (A.IV) is a standard assumption, which holds for instance in the case of PCA-Net and DeepONet. The former also enforces the learned encoder E X to be linear, which is not needed in our setup. Moreover, both approaches usually only deal with the case where both X and Y are Hilbert spaces. Problems such as (2.6) are ubiquitous in scientific computing, with many applications in engineering, biology, physics, finance and beyond. In many such applications, it is common to assume that the measure µ on X is the law of a random fieldEQUATIONfor functions a 0 , ϕ i ∈ X , where the x i are random variables and c i ≥ 0 are scalars that ensure that a ∈ L ∞ (Ω). Under some mild assumptions, (2.8) is then the Karhunen-Loève (KL) expansion of the measure µ. See, e.g., [91] (see also [55, §3.5.1] ). The x i are typically independent. While in some settings, they may have infinite support, it is also common in practice to assume they range between finite maxima and minima. After rescaling, one may therefore assume thatx ∈ D = [-1, 1] N .Problems of this type fits into our framework. Suppose that x = (x 1 , x 2 , . . .) ∼ ϱ. The measure µ is then given as the pushforward µ = a♯ϱ and f :D → Y is the parametric solution map f : x ∈ D → u(•; a(x)) ∈ Y.If needed, the map ι can be defined in a number of different ways. Suppose, for instance, that X is a Hilbert space, e.g., X = L 2 (Ω), and {ϕ i } ∞ i=1 is a Riesz system (this holds, for instance, in the case of a KL expansion, in which case {ϕ i } ∞ i=1 is an orthonormal basis). Then {ϕ i } ∞ i=1 has a unique biorthogonal dual Riesz system {ψ i } ∞ i=1 . We may therefore defineι : a → ⟨a -a 0 , ψ i ⟩ L 2 (Ω) /c i ∞ i=1 . Notice that ι is a bounded linear map and F (X) = f • ι(X) = f (x) for X = a(x) ∼ µ.However, evaluating ι is often not required for computations (see §A.1). This example considers an affine parametrization (2.8) inducing the measure µ. Note that other parametrizations can be considered. Common examples include the quadratic a(z;x) = a 0 (z) + ( ∞ i=1 c i x i ϕ i (z)) 2 and log-transformed a(z, x) = exp ( ∞ i=1 c i x i ϕ i (z)) parametrizations [17].Remark 2.4 (Holomorphy assumption) In the previous example, the operator F stems from the solution map f : D → Y of a parametric PDE. The regularity of solution maps of parametric PDEs has been intensively studied, and it is known that many such maps are (b, ε)-holomorphic (hence the resulting operator satisfies (A.II)). Consider, for instance, the affine diffusion problem (2.7)-(2.8). Under a mild uniform ellipticity condition, the solution map of the standard weak form of the PDE f :x ∈ D → u(a(•; x)) ∈ H 1 0 (Ω) is (b, ε)-holomorphic with b = (b i ) ∞ i=1 and b i = c i ∥ϕ i ∥ L ∞ (Ω). See, e.g., [3, Prop. 4 .9], as well as §B.3. Similar results are known for other parametric PDEs. This includes parabolic PDEs, various types of nonlinear, elliptic PDEs, PDEs over parametrized domains, parametric hyperbolic problems and parametric control problems. See [19] or [3, Chpt. 4] for reviews. 
main results i: upper bounds: 
 We now present our first two main results. In these results, given an optimization problem min t f (t), we say that t is a τ -approximate minimizer for some τ ≥ 0 if f ( t) ≤ min t f (t) + τ 2 . Theorem 3.1 (Existence of good DNN architectures). Let m ≥ 3, δ > 0, 0 < ϵ < 1 and L = L(m, ϵ) = log 4 (m) + log(1/ϵ). Then there exists a class N of hyperbolic tangent (tanh) DNNs N : R d X → R d Y depending on m and ϵ only withEQUATIONsuch that following holds. Suppose that Assumption 2.2 holds andEQUATIONwhere I X : X → X is the identity map and c > 0 is a universal constant. Let X 1 , . . . , X m ∼ i.i.d. µ and consider the noisy training data (2.4) with arbitrary noise E i ∈ Y. Then, with probability at least 1 -ϵ, every τ -minimizer N of (2.5), where τ ≥ 0 is arbitrary, yields an approximation F that satisfiesEQUATIONEQUATIONand, if Y is a Hilbert space,∥F -F ∥ L 2 µ (X ;Y) ≲ E app,2 + E X ,2 + E Y,2 + E opt,2 + E samp,2 . (3.5)Here, the approximation error terms E app,q , q = 2, ∞, are given by EQUATIONEQUATIONwhere ) approximate the identity maps on X and Y, respectively. Observe that these terms are increasing in m for fixed encoders and decoders. Therefore, as one expects, the accuracy of the encoder-decoder approximationsEQUATIOND X • E X ≈ I X and D Y • E Y ≈ I Y shouldincrease with increasing m to ensure decay to zero of the generalization error as m → ∞. The specific terms in (3.7) (for q = 2) are quite standard in operator learning. See, e.g., [53, 55] . When the encoders and decoders are computed via PCA, as in PCA-Net, standard bounds can be derived for these terms [53] . For similar analysis in the case of DeepONets, see [55] . Finally, the error (3.3)-(3.5) involves an optimization error E opt , which primarily depends on how accurately the optimization problem (2.5) is solved (i.e., the term τ ), and a sampling error E samp , which depends on the error in the training data (2.4).Theorem 3.1 allows Y to be a Banach or a Hilbert space. Overall, when Y is only a Banach space, we obtain a weaker L 2 µ -norm bound involving the Pettis norm (3.3) and, moreover, the approximation error E app,q is worse by a factor of 1/2 than when Y is a Hilbert space. (Note that one can establish a bound for the Bochner L 2 µ -norm error when Y is a Banach space via (3.4) and the inequality∥ • ∥ L 2 µ (X ;Y) ≤ ∥ • ∥ L ∞ µ (X ;Y). However, we do not believe the resulting bound is sharp). As we discuss in Remark D.18, the discrepancies between the two cases stem from the lack of an inner product structure and, in particular, the absence of Parseval's identity when Y is a Banach space.Observe that the DNN architecture in Theorem 3.1 is independent of the smoothness of the operator being learned. We term such an architecture problem agnostic. This theorem considers tanh activations only. However, as we discuss in Remark D.11, other activations can be readily used instead. Other key facets of Theorem 3.1 are the width and depth bounds (3.1). Qualitatively, these agree with empirical practice: namely, better performing DNNs tend to be wider than they are deep, and relatively shallow DNNs perform well in practice (see [24, 25] and references therein). We also see this later in §5.On the other hand, the family N is not fully connected. As we describe in §C.2.1, while the weights on the final layer can be arbitrary real numbers, the weights and biases in the hidden layers come from a finite (but large) set: they are handcrafted to approximately emulate certain multivariate orthogonal polynomials. Since fully-connected DNNs are typically used in practice, Theorem 3.1 is essentially a theoretical contribution. In our next result, we consider the more practical scenario of fully-connected DNNs. Theorem 3.2 (Fully-connected DNN architectures are good). There are universal constants c 1 , c 2 , c 3 , c 4 ≥ 1 such that the following holds. Let m, δ, ϵ and L be as in Theorem 3.1,d X ≥ c 1 (m + log(1/ϵ)), L ι • ∥I X -D X • E X ∥ L 2 µ (X ;X ) ≤ c(δ) • (m + log(1/ϵ)) -1/2 , (3.9)where c(δ) > 0 depends on δ only, consider any class N of fully-connected DNNs satisfying(n 0 , n L+2 ) = (d X , d Y ), N 1 , . . . , N L+1 ≥ c 2 • (m + log(1/ϵ)) • (m/L) δ , L ≥ c 3 • log(m/L).(3.10) Suppose that Assumption 2.2 holds and that the pushforward ς in (A.I) is the tensor-product of a univariate probability distribution with mean zero and variance ω ≳ 1. Let X 1 , . . . , X m ∼ i.i.d. µ and consider (2.4) with arbitrary E i ∈ Y. Then the following hold with probability at least 1 -ϵ.(A) Uncountably many 'good' minimizers. The problem (2.5) has uncountably many minimizers that satisfy (3.3) with τ = 0 or (3.5) with τ = 0 if Y is a Hilbert space. They also satisfy (3.4) with τ = 0 and the modified right-hand side√ LE app,∞ +LE X ,∞ + √ LE Y,∞ +E opt,∞ + √ LE samp,∞ .(B) Good minimizers are stable. Suppose that E X ∈ L ∞ µ (X ; R d X ) and let τ o > 0 be arbitrary. Then there is a neighbourhood of DNN parameters around the parameters of each minimizer in (A) for which the approximation corresponding to any parameters in this neigbourhood also satisfies the same bounds as in (A) with τ = τ o .(C) Good minimizers can be far apart in parameter space. For sufficiently large m, there are at least (m/(c 4 L)) 2δm minimizers satisfying the bounds in (A) such that, for any two such minimizers, their parameters satisfy ∥θ ′ ∥ = ∥θ∥ and ∥θ ′ -θ∥ ≳ 1.This theorem states that DL with fully-connected DNN architectures of sufficient width and depth (3.10) can succeed, since there are minimizers that yield the optimal bounds of Theorem 3.1. Such minimizers are uncountably many in number (A), stable to perturbations (B) and many of them (exponentially in m) have sufficiently distinct and nonvanishing/nonexploding parameters (C). This theorem does not imply that all minimizers are 'good' -an issue we discuss further in §6 -but our numerical results in §5 suggest that (approximate) minimizers obtained through training do, at least for the experiments considered, achieve the rates specified in Theorem 3.1. 
main results ii: lower bounds: 
 We now show that the various approximation errors are nearly optimal. For this, we ignore the encoding-decoding, optimization and sampling errors and proceed as follows. Let C(X ; Y) be the Banach space of continuous operators. We term an (adaptive) sampling map as any mapEQUATIONwhere X 1 ∈ X , X 2 = X 2 (F (X 1 )) ∈ X potentially depends on the previous evaluation F (X 1 ), X 3 = X 3 (F (X 1 ), F (X 2 )) ∈ X , and so forth. Next, we term a reconstruction map as any map R :EQUATIONwhere the infimum is taken over all such L and R. In other words, θ m (b) measures how well one can learn holomorphic operators using arbitrary training data and an arbitrary reconstruction procedure. Theorem 4.1 (Optimal L 2 error rates). Suppose that (A.I) holds. Then, for any 0 < p < 1 there is a constant c(p) > 0 such that the following hold.(i) For each m ∈ N, there is a b ∈ ℓ p M (N), b ≥ 0, ∥b∥ p,M = 1 such that θ m (b) ≥ c(p) • m 1/2-1/p . (ii) There is a b ∈ ℓ p M (N), b ≥ 0, ∥b∥ p,M = 1 such that θ m (b) ≥ c(p) • m 1/2-1/p log 2/p (2m) , ∀m ∈ N.This theorem shows that the error E app,2 in Theorems 3.1-3.2 is optimal, up to log terms, whenever Y is a Hilbert space: there does not exist a reconstruction map surpasses the rate m 1/2-1/p for learning holomorphic operators. Note that this result applies not only to DL-based procedures, but any procedure that learns such operators from m samples. Another consequence of this result is that adaptive sampling, i.e., active learning, is of no benefit. As shown by Theorems 3.1-3.2, the optimal rate m 1/2-1/p can, up to log terms, be achieved through inactive learning, i.e., i.i.d. sampling from µ. EQUATIONwhere the infimum is taken over all reconstruction maps R : (X × Y) m → L ∞ µ (X ; Y) only. Theorem 4.2 (Optimal L ∞ error rates). Suppose that (A.I) holds and that the pushforward ς is the tensor-product of a univariate probability distribution with mean zero and variance ω ≳ 1. Then, for any 0 < p < 1 there is a constant c(p) > 0 such that the following hold.(i) For each m ∈ N, there is a b ∈ ℓ p M (N), b ≥ 0, ∥b∥ p,M = 1 such that θm (b) ≥ c(p) • m 1-1/p log(m) . (ii) There is a b ∈ ℓ p M (N), b ≥ 0, ∥b∥ p,M = 1 such that θm (b) ≥ c(p) • m 1-1/p log 2/p+1 (2m) , ∀m ∈ N.As with the previous theorem, this result asserts that the rate m 1-1/p is optimal in the L ∞ -norm when Y is a Hilbert space. However, it is strictly weaker than Theorem 4. Since the main objective of this work is to examine the approximation error, we follow a standard setup and fix the encoder and decoders for each experiment, so that E X and D Y in (1.3) do not change for different choices of N . We also set up our experiments so that encoding-decoding (3.7) and sampling (3.8) errors are zero. We do this in a standard way. To ensure that E X ,q = 0, we truncate the parametric expansions (2.8) after d terms (henceforth termed the parametric dimension) and define the encoder E X accordingly. This means we effectively consider a parametric PDE depending on finitely-many parameters. We use Finite Element Methods (FEMs) to both solve the PDE (for generating training and testing data) and define the decoder D Y (see (A.2)). To ensure that E Y,q = 0, we compute errors with respect to the Bochner L 2 µ (X ; Y)-norm, where Y = D Y (R d Y ) is the FEM discretization of Y. In other words, we use the same FEM code to generate test data and compute the errors as we do to construct the operator approximation F . See §A.1 for further details.The DNNs in our experiments are fully-connected and of the form (2.1). We denote by σ L × N DNN a DNN N with activation function σ, width N and depth L. To solve (2.5) we use Adam [47] with early stopping and an exponentially decaying learning rate. We train our DNN architectures for 60,000 epochs and results are averaged over a number of trials. See §A.2 for further details.Parametric elliptic diffusion equation. Our first example is the parametric elliptic diffusion equation (2.7) . This PDE arises in many scientific computing applications, such as groundwater flow modelling, see, e.g., [95] . We describe the full PDE and its FE discretization in §B.3. In our experiments, we consider both affine (B.1) and log-transformed (B.2) diffusion coefficients. The latter is particularly useful in the groundwater flow problem as the permeability of various layers of sediment can vary on logarithmic scales. Differing from most prior work, we consider a novel mixed variational formulation [32] of (2.7), which has a number of key practical benefits (see §B.3.1). In this case, Y = L 2 (Ω) is a Hilbert space. Fig. 1 compares the error versus the amount of training data m for various DNN architectures for learning the solution map of this PDE in d = 4 and d = 8 parametric dimensions with these two diffusion coefficients. We observe that architectures with the Exponential Linear Unit (ELU) or hyperbolic tangent (tanh) activation generally outperform similar architectures with the Rectified Linear Unit (ReLU) activation (as we discuss in Remark D.11, this difference is in agreement with our theoretical analysis). Overall, the best performing DNNs appear to roughly match the plotted rate m -1 . As we explain further in §B.3.2, this rate is precisely that predicted by our theory. In particular, the parametric solution map (recall Remark 2.4) is (b, ε)-holomorphic with b ∈ ℓ p M (N) for any p < 2/3, giving an effective convergent rate m 1/2-1/p that is arbitrarily close to m -1 . Another important fact that we observe is that despite the parametric dimension doubling from 4 to 8, there is little change in the error behaviour.Parametric Navier-Stokes-Brinkman equations. We next consider the parametric Navier-Stokes-Brinkman (NSB) equations. See §B.4 and (B.14) for the full definition. Here the solution is a pair (u, p), where u is the velocity field and p is the pressure. These equations describe the dynamics of a viscous fluid flowing through porous media with random viscosity. See, e.g., [28, 43, 46, 94] . We use a mixed variational formulation [34] to discretize the PDE. This formulation is more sophisticated that standard variational formulations, but conveys various practical advantages. Unlike the previous example, it leads to Y being either Y = L 4 (Ω) for u or Y = L 2 (Ω) for p. See §B.4.1 for details. Fig. 2 compares a variety of DNN architectures for approximating the velocity field component in d = 4 and d = 8 parametric dimensions. Here again we observe the ELU and tanh DNN architectures outperform similar sized ReLU architectures. We also observe a rate close to m -1 . Note that it is currently unknown whether this or the next example possess the same (b, ε)-holomorphy guarantee as that of the previous example. Yet we observe the same rate, and therefore conjecture that such a property does indeed hold in these cases. Similar to the previous example, there is also no deterioration of the rate when moving from d = 4 to d = 8.Parametric stationary Boussinesq equation. Our final example is a parametric stationary Boussinesq PDE. See §B.5 and (B.16) for the full definition. Here the solution is a triplet (u, φ, p), where u is the velocity field, φ is the temperature and p is the pressure of the solution. The Boussinesq model arises in a variety of engineering, fluid dynamics and natural convection problems where changes in temperature affect the velocity of a fluid [14, 22, 39] . Similar to the previous example, we consider a fully mixed variational formulation (see §B.5.1), which leads toY = L 4 (Ω) (for u), Y = L 4 (Ω) (for φ) or Y = L 2 0(Ω) (for p). Fig. 3 provides numerical results. Our observations are in line with the previous two examples, with the ELU and the smaller tanh networks being most often the best performers in this problem. Once more, the errors roughly correspond to the rate m -1 and there is no deterioration with increasing d. 
conclusions and limitations: 
 The purpose of this work was to derive near-optimal generalization bounds for learning certain classes of holomorphic operators that arise frequently in operator learning tasks involving PDEs. Complementing and extending previous works [26, 31, 41, 55, 71] on the approximation of such operators via DNNs, we showing sharp algebraic rates of convergence in m, thus confirming that such operators can be learned efficiently and without the curse of dimensionality. It is notable that the sizes of the various DNNs in Theorems 3.1-3.2 also do not succumb to the so-called curse of parametric complexity [54] , since the width and depth bounds are at most algebraic in m.We end by discussing a number of limitations. First, assumption (A.I) may not hold in some applications. The domain D can easily be replaced by bounded hyperrectangle through rescaling and the condition that ς be quasi-uniform relaxed to quasi-ultraspherical (by considering ultraspherical polynomials). However, it is currently an open problem whether our results can be extended to the case where µ is Gaussian, in which case ς would typically be a tensor-product Gaussian measure on R N and the relevant polynomials would be the Hermite polynomials. Second, the reader may have noticed that the encoder E X defined in (A.III) and used to construct the approximation (2.5) involves the pair ( E X , D X ) and the map ι d X . This is a technical requirement -also found in other theoretical works on operator learning -needed to obtain encoding-decoding errors of the form E X ,q , q = 2, ∞. It is unknown whether it can be relaxed. It is also unknown whether the assumption on ς in (A.III) can be relaxed. We believe this can be done, at least if the L 2 µ -norm in (3.2) is replaced by the L ∞ µ -norm. Whether this is possible without modifying (3.2) is currently unknown. Third, a limitation of Theorem 3.2 is that it only asserts that some minimizers are 'good', not all. Techniques from statistical learning theory can provide stronger bounds that hold for all minimizers. Yet, as noted in §1.2, these tools typically produce slower rates of decay in m. Overcoming this limitation -e.g., by refining these tools for the holomorphic setting or showing that the 'good' minimizers can indeed be obtained via standard training -is a topic of future work.Finally, as noted, our theorems provided worse generalization bounds when Y is a Banach space than when Y is a Hilbert space. Our numerical results in Figs. 2 3 suggest that this factor is an artefact of the proofs. Whether it can be removed is an interesting open problem. 
 
 Appendix: In this section, we describe our experimental setup.We first explain how all our experiments are formulated as operator learning problems. We do this in a standard way, by first truncating the random field which generates the measure µ, and then using an FEM to discretize the output space.All our examples follow Example (2.3) and involve operators of the form (2.6) , where F a represents a different type of PDE in each case (specifically, either an elliptic diffusion, Navier-Stokes-Brinkman or Boussinesq PDE). We consider both affine (2.8) and log-transformed parametrizations of the random field a(x). Thus, in general we writea(x) = a(•; x) = g a 0 (•) + ∞ i=1 c i x i ϕ i (•) , (A.1)where g : R → R is a (measurable) map. In the affine case g(t) = 1. In the log-transformed case g(t) = exp(t).As discussed, since the main objective of this work is to examine the approximation error, we set up our experiments so that encoding-decoding errors are zero. We do this as follows. First, we fix a parametric dimension d and truncate the expansion in (A.1) after d terms, giving a mapa d : [-1, 1] d → X and measure µ = a d ♯ϱ d , where ϱ d is the uniform probability measure on [-1, 1] d .We then define the operator F asF (a d (x)) = f (x) = u(•; a d (x)), where u(•; a) is the solution of the PDE F a u = 0.In alignment with our theorems, we focus on the in-distribution performance of the learned approximation F . This means we define the encoder only on supp(µ), asE X (X) = x when X = a d (x) with x ∈ [-1, 1] d .As a result, the encoding-decoding error E X ,q in (3.7) satisfies E X ,q = 0.To perform our simulations, we use FEMs to solve the PDE and discretize the output space Y. Let {φ i } K i=1 ⊂ Y be a FEM basis and set d Y = K. Then we define the decoder asD Y : R K → Y, D Y (c) = K i=1 c i φ i (A.2) and set Y = D Y (R K ) as the discretization of Y.With this in hand, we now describe the simulation of training data in general terms. First, we draw x 1 , . . . , x m ∼ i.i.d. ϱ d . Then, for each training sample x i , we compute Y i by using the FEM to solve the PDE with parameterX i = a d (x i ). Notice that Y i ∈ Y in this setup.We consider a DNN architecture with input dimension n 0 = d and output dimension n L+2 = K.After training, we evaluate the learned approximationF = D Y • N • E X over supp(µ) as F (X) = D Y • N • E X (X) = D Y • N (x) for X = a d (x) with x ∈ [-1, 1] d .Finally, as noted, we us the same FEM discretization to generate testing data, which allows us to measure the error with respect to the L 2 µ (X ; Y)-norm. This effectively means that the encoding-decoding error E Y,q in (3.7) satisfies E Y,q = 0 as well.In this work, we investigate the trade-off between the accuracy of the learned operator and the number of samples m used in training. Our methodology is summarized as follows.(i) Implementation. Trained models were not retained due to space limitations on the cluster. (iii) Choice of architectures and initialization. Based on the strategies in [1] , we fix the number of nodes per layer N and depth L such that the ratio β := L/N is β = 0.5. In addition, we initialize the weights and biases using the HeUniform initializer from keras setting the seed to the trial number. We consider the Rectified Linear Unit (ReLU)σ 1 (z) := max{0, z}, hyperbolic tangent (tanh) σ 2 (z) := e z -e -z e z + e -z , or Exponential Linear Unit (ELU) σ 3 (z) = z z > 0, e z -1 z ≤ 0activation functions in our experiments. (iv) Optimizers for training and parametrization. To train the DNNs, we use the Adam optimizer [47] , incorporating an exponentially-decaying learning rate. We train our models for 60,000 epochs or until converging to a tolerance level of ϵ tol = 5 • 10 -7 in single precision. In light of the nonmonotonic convergence behavior observed during the minimization of the nonconvex loss (see, e.g., [1, 2] ), we implement early stopping. More precisely, we save the weights and biases of the partially trained network once the ratio between the current loss and the last checkpoint loss is reduced below 1/8, or if the current weights and biases produce the best loss value observed in training. We then restore these weights after training only if the loss value of the current weights is larger than that of the saved checkpoint. (v) Training data and design of experiments. First, we define a 'trial' as a complete training run for a DNN approximating a specific function, initialized as mentioned above.Following the setup of §A.1, we run several trials solving the problem:Given training data {(X i , Y i )} m i=1 ⊂ (X × Y) m , X i ∼ i.i.d. µ, Y i = F (X i ) + E i ∈ Y, approximate F ∈ L 2µ (X ; Y). We generate the measurements Y i using mixed variational formulations of the parametric elliptic, Navier-Stokes-Brinkman and Boussinesq PDEs discretized using FEniCS with input data X i . The noise E i ∈ Y encompasses the discretization errors from numerical solution. Further details of the discretization can be found in §B. Each of our architectures is trained across a range of datasets with increasing sizes. This involves using a set of training data consisting of values {(X i , Y i ))} m i=1 , where m denotes the size of the training data and belongs to the set {10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 200, 300, 400, 500}. After training we calculate the testing error for each trial and run statistics across all trials for each dataset. (vi) Testing data and error metric. The testing data is generated similarly to the training data, obtaining solutions at different points X i ∈ X for i = 1, . . . , m test . However, the testing data{(X i , Y i = F (X i )+E i )} mtest i=1is generated using a deterministic high-order sparse grid collocation method [76] . In particular, we use sparse grid quadrature rules to compute approximations to the Bochner norms∥F ∥ L 2 µ (X ; Y) = X ∥F (X)∥ 2 Y dµ(X)) 1/2 ≈ mtest i=1 ∥F (X i )∥ 2 Y w i 1/2, where µ = a d ♯ϱ d is the pushforward measure defined in (A.1) and w i , i = 1, . . . , m test are the quadrature weights. We use these approximations to compute the relative L 2 µ (X ; Y) errore test F = mtest i=1 ∥F (X i ) -F (X i )∥ 2 Y w i 1/2 mtest i=1 ∥F (X i )∥ 2 Y w i 1/2.We use a high order isotropic Clenshaw Curtis sparse grid quadrature rule to evaluate e test F , as described in [2] . This method shows superior convergence over Monte Carlo integration to evaluate the global Bochner error. The sparse grid rule gives m test points at a level ℓ for d dimensions. We rely on the TASMANIAN sparse grid toolkit [88] [89] [90] for the generation of the isotropic rule to study the generalization performance of the DNN. (vii) Visualization. The graphs in Figs. 1 2 3 show the geometric mean (the main curve) and plus/minus one (geometric) standard deviation (the shaded region). We use the geometric mean because our errors are plotted in logarithmic scale on the y-axis. See [3, Sec. A.1] for further discussion about this choice.In this section, we provide full details of the parametric PDEs considered in our numerical experiments. We also describe their variational formulations and numerical solution using FEM.We consider two parametric coefficients of the form (A.1) in our numerical experiments. Our first coefficient isa 1 (z, x) = 2.62 + ∞ j=1 x j sin(πz 1 j) j 3/2 , ∀z ∈ Ω, (B.1)where x j ∈ [-1, 1], ∀j. Our second example involves a log-transformed coefficient, which is a rescaling of an example from [76] of a diffusion coefficient with one-dimensional (layered) spatial dependence given bya 2 (z, x) = exp   1 + x 1 √ πβ 2 1/2 + ∞ j=2 ζ j ϑ j (z)x j   , ∀z ∈ Ω, ζ j := ( √ πβ) 1/2 exp -(⌊j/2⌋πβ) 2 8 , ϑ j (z) := sin (⌊j/2⌋πz 1 /β p ) if j is even cos (⌊j/2⌋πz 1 /β p ) if j is odd , (B.2)wherex j ∈ [-1, 1], ∀j.Here we let β c = 1/8, andβ p = max{1, 2β c }, β = β c /β p .In both cases we consider truncation of the expansion after d terms, giving the map a j,d : [-1, 1] d → X , with j = 1 corresponding to (B.1) and j = 2 corresponding to (B.2). Our input samples are then{X i = a j,d (x i )} m i=1 ⊂ X with x i ∈ [-1, 1]d drawn identically and independently from ϱ d and j ∈ {1, 2}. Note for the Boussinesq problem we also consider an additional parametric dependence in the tensor K describing the thermal conductivity of the fluid. See §B.5 and (B.17).We require several function space definitions for the development of the mixed variation formulations of the Poisson, Navier-Stokes-Brinkman and Boussinesq PDEs. Let Ω ⊂ R n , n ∈ {2, 3}, be the physical domain of a PDE. We write L p (Ω), 1 ≤ p ≤ ∞, for the L p -space of scalar-valued functions with respect to the Lebesgue measure (to avoid confusion with the Bochner space L 2 µ (X ; Y)). We denote the standard Sobolev spaces as W s,p (Ω) for s ∈ R and p > 1, and write H k (Ω) when p = 2 and s = k. Additionally, we consider the space of traces of functions in H 1 (Ω), denoted by H 1/2 (∂Ω), and its dual, H -1/2 (∂Ω) (see, e.g., [11, Sec. 1.2] for further details). We also define the following closed subspace of H 1 (Ω):H 1 0 (Ω) := C ∞ 0 (Ω) ∥•∥1,Ω.HereC ∞ 0 (Ω) ∥•∥1,Ωdenotes the closure of C ∞ 0 (Ω) (i.e., the space of C ∞ (Ω) functions with compact support) with respect to the norm ∥ • ∥ 1,Ω , which, for any v ∈ H 1 (Ω), is given by∥v∥ 1,Ω := |v| 2 1,Ω + ∥v∥ 2 L 2 (Ω) 1/2, where |v| 1,Ω := ∥∇v∥ L 2 (Ω) .For scalar functions u and vector fields v, we use ∇u and div(v) to denote their gradient and divergence, respectively. For tensor fields σ and τ , represented by (σ i,j ) n i,j=1 and (τ i,j ) n i,j=1 , respectively, we define div(σ) as the divergence operator div acting along the rows of σ, and we define the trace and the tensor inner-product astr(σ) = n i=1 σ i,i , and τ : σ = n i,j=1 τ i,j σ i,j ,respectively. Furthermore, we introduce the notation L p (Ω) and L p (Ω) to represent the vectorial and tensorial counterparts of L p (Ω), respectively, and H 1 (Ω) and H 1 (Ω) for the vectorial and tensorial counterparts of H 1 (Ω), respectively. Keeping this in mind, we introduce the Banach spacesH(div q ; Ω) := v ∈ L 2 (Ω) : div(v) ∈ L q (Ω) , H(div q ; Ω) := τ ∈ L 2 (Ω) : div(τ ) ∈ L q (Ω) (B.3) with norms ∥v∥ H(divq;Ω) := ∥v∥ L 2 (Ω) + ∥div(v)∥ L q (Ω) , ∥τ ∥ H(divq;Ω) := ∥τ ∥ L 2 (Ω) + ∥div(τ )∥ L q (Ω) .The cases of q = 4/3 and q = 2 appear in the mixed variational formulations of the considered PDEs, and for the latter we simply write H(div; Ω).Often, under certain conditions, such as incompressibility conditions [33, eq.(2.4)], it is convenient to define variants of these spaces. For example, we defineL 2 tr (Ω) := τ ∈ L 2 (Ω) : tr(τ ) = 0 , (B.4)which represents the space of integrable functions with zero trace over Ω. Furthermore, given the decomposition (see, e.g., [32] )H(div 4/3 ; Ω) = H 0 (div 4/3 ; Ω) ⊕ R I , (B.5)we may also considerH 0 (div 4/3 ; Ω) := τ ∈ H(div 4/3 ; Ω) : Ω tr(τ ) = 0 , (B.6)as the space of elements in H(div 4/3 ; Ω) with zero mean trace. Finally, we defineL 2 skew (Ω) = {η ∈ L 2 (Ω) : η + η t = 0}, (B.7)and the space of L 2 (Ω) functions with zero integral over Ω asL 2 0 (Ω) = ν ∈ L 2 (Ω) : Ω ν = 0 . (B.8)We now describe our first example, which is the parametric elliptic diffusion equation. Let Ω ⊂ R 2 be a bounded Lipschitz domain, ∂Ω be the boundary ofΩ, f ∈ L 2 (Ω) and g ∈ H 1/2 (∂Ω). Given x ∈ [-1, 1] d , we consider the PDE -div(a(x)∇u(x)) = f, in Ω (B.9) u(x) = g, on ∂Ω.Here a(x) = a(•, x) ∈ L ∞ (Ω) =: X is the parametric diffusion coefficient. The terms f and g are not parametric.Our first step in precisely defining the problem is to identify sufficient conditions for the solution map x → u(x) to be well-defined. To do this, we diverge from the standard variational formulation involving the space Y = H 1 0 (Ω) (see Remark B.3) and instead consider a mixed variational formulation of (B.9). Using a mixed formulation to study the solution of PDEs offers several benefits over the standard variational formulation. One key advantage is that it allows us to introduce additional variables that can be of physical interest. Additionally, mixed formulations can naturally accommodate different types of boundary conditions and introduce Dirichlet boundary conditions directly into the formulation rather than imposing them on the search space. For further details on mixed formulations, we refer to [32] and references within.Assume that there exists r, M > 0 such that, for allx ∈ [-1, 1] d , 0 < r ≤ essinf z∈Ω a(z, x) =: a min (x) and a max (x) := esssup z∈Ω a(z, x) ≤ M. (B.10)Then the problem can be recast as a first-order system: givenx ∈ [-1, 1] d , find (σ, u)(x) ∈ H(div; Ω) × L 2 (Ω) such that d a(x) (σ(x), τ ) + b(τ , u(x)) = G(τ ), ∀τ ∈ H(div; Ω), b(σ, v) = F (v), ∀v ∈ L 2 (Ω). (B.11)Here d and b are the bilinear forms defined byd a(x) (σ, τ ) = Ω σ • τ a(x) , ∀(τ , σ) ∈ H(div; Ω) × H(div; Ω), b(τ , v) = Ω div(τ )v, ∀(τ , v) ∈ H(div; Ω) × L 2 (Ω)and the functionals G ∈ (H(div; Ω)) ′ and J ∈ (L 2 (Ω)) ′ are defined byJ(v) = - Ω f v, ∀v ∈ L 2 (Ω) and G(τ ) = ⟨γ(τ ) • n, g⟩ 1/2,∂Ω , ∀τ ∈ H(div; Ω). (B.12)For the experiments in this work, we consider Ω = (0, 1) 2 and f = 10. For the boundary condition, we consider a constant value u(z, x) = 0.5 on the bottom of the boundary (0, 1) × {0}, and zero boundary conditions on the remainder of the boundary.Consider the affine parametrization (B.1). Setting M = 2.7 and observing that j∈Nx j sin(πz 1 j) j 3/5 ≤ j∈N 1 j 3/5 ≈ 2.61238 = 2.62 -r, (B.13)for some r < 0.00762, we deduce that (B.10) holds, which makes the mixed variational formulation (B.11) well defined, i.e., for each x ∈ [-1, 1] ∞ , there exists a unique solution (σ, u)(x) ∈ H(div; Ω) × L 2 (Ω). Moreover, one can show that the parametric solution map In view of this property, this example falls within our theory. Note that b ∈ ℓ p M (N) for every p < 2/3. Thus, we expect a theoretical rate of convergence with respect to the amount of training data that is arbitrarily close to m 1/2-3/2 = m -1 . This holomorphy result applies to the affine diffusion (B.1), not the log-transformed diffusion (B.2). However, we expect that it is possible to extend [74, Prop. A.3.2] to the latter case.x → (σ, u)(x) is (b, ε)-holomorphic for 0 < ε < 0.00762 and where b = (b i ) ∞ i=1 is given by b j = ∥ sin(πj•)/j 3/2 ∥ L ∞ (Ω) = j -3/2 . See [74, Prop. A.3.2].We use so-called conforming Finite Element (FE) discretizations [18, Chp. 3] . Given a number of Degrees of Freedom (DoF) K, this results in finite-dimensional spaces H K ⊆ H(div; Ω) and Q K ⊆ L 2 (Ω). Specifically, we consider a regular triangulation T K of Ω made up of triangles of minimum diameter h min = 0.0844 and maximum diameter h max = 0.1146. This corresponds to a total number of DoF K = 2622. See Fig. 4 for an illustration of the FE mesh. u : [-1, 1] d × Ω → R 2 and p : [-1, 1] d × Ω → R such that ηu -λdiv(a(x)e(u(x))) + (u(x) • ∇)u(x) + ∇p(x) = f in Ω div(u(x)) = 0 in Ω (B.14) u = u D on ∂Ω in 0 on ∂Ω wall (a(x)∇e(u) -pI)ν = 0 on ∂Ω out Ω p = 0,where λ = Re -1 and Re is the Reynolds number, a(x) = a(•; x) ∈ X := L ∞ (Ω) and a :[-1, 1] d × Ω → R +is the random viscosity of the fluid, η ∈ R + is the scaled inverse permeability of the porous media, u is the velocity of the fluid, e(u) = 1 2 (∇u + (∇u) t ) is the symmetric part of the gradient, p is the pressure of the fluid and f : Ω → R is an external force independent of the parameters. Here, the fourth condition imposes a zero normal Cauchy stress (a(x)∇e(u) -pI)ν = 0 for the output boundary ∂Ω out .In addition, the incompressibility of the fluid imposes the following compatibility condition on u D :Γ u D • n = 0 on ∂Ω in .The third condition also imposes a no-slip condition on the walls Ω wall [34, eq.( 2.3)].The analysis of the detailed mixed formulation used for this problem in the nonparametric case can be found in [34] . Over the last decade, many works have used a mixed formulation employing a Banach space framework, allowing one to solve different PDEs in continuum mechanics in suitable Banach spaces. The advantage of this formulation is that no augmentation is required, the spaces are simpler and closer to the original model, and it allows one to obtain more direct approximations of the variables of physical interest [34, Sec. 1] .Based on the analysis in [34] , the mixed variational formulation of the parametric NSB equations in (B.14) becomes: givenx ∈ [-1, 1] d , find (u, t, σ, γ)(x) ∈ L 4 (Ω) × L 2 tr (Ω) × H 0 (div 4/3 ; Ω) × L 2 skew (Ω) such that λ Ω a i (x)t(x) : s - Ω s : σ(x) - Ω (u ⊗ u)(x) : s = 0, Ω t(x) : τ + Ω γ(x) : τ + Ω u(x) • div(τ ) = ⟨τ n, u D ⟩ ∂Ωin , Ω δ : σ(x) + Ω v • div(σ(x)) - Ω ηu(x) • v = Ω f • v, (B.15) for all (v, s, τ , δ) ∈ L 4 (Ω) × L 2 tr (Ω) × H 0 (div 4/3 ; Ω)× L 2 skew (Ω).Numerically, the skew-symmetry of γ is imposed by searching for γ ∈ L 2 (Ω) and settingγ = 0 γ -γ 0 .Moreover, we impose the Neumann boundary condition via a Nietsche method as in [34, Sec. 5.2] . Specifically, we add κ⟨(σ + u ⊗ u)n), τ n⟩ ∂Ωout = 0 to the second equation where κ ≫ 1 is a large constant (e.g., κ = 10 4 ). As usual in this formulation, the pressure p ∈ L 2 (Ω) can be computed according to the post-processing formula Note that above we omitted the term x for simplicity.p = - 1 2 tr(σ + (u ⊗ u)).In our experiments, we consider approximating solutions to the parametric NSB problem with λ = 0.1, a scaled inverse permeability of η = 10 + z 2 1 + z 2 2 , an external force f = (0, -1) ⊤ and random viscosity a j,d as in (B.1)-(B.2) with j ∈ {1, 2}.We consider the unit square Ω = (0, 1) 2 as the domain, an inlet boundary ∂Ω in = (0, 1) × {1}, an outlet boundary ∂Ω out = {1} × (0, 1) and walls ∂Ω wall = {0} × (0, 1) ∪ (0, 1) × {0}. For simplicity, we use the same mesh as that of the previous example. See Fig. 4 . On the Neumann boundary ∂Ω out we consider a zero normal Cauchy stress, a Dirichlet condition u D = (0.0625) -1 ((z 2 -0.5)(1z 2 ), 0) on ∂Ω in and a no-slip velocity on ∂Ω wall .Fig. 6 provides a comparison between a reference solution of the vector field u and pressure p computed by the FEniCS FEM solver and the approximation generated by a ELU 4 × 40 DNN.We report the performance of the DNNs approximating (u, p)(x) ∈ L 4 (Ω) × L 2 (Ω). Note that any solver based on the above formulation outputs several other variables, e.g., (t, σ, γ)(x) ∈ L 2 tr (Ω)×H 0 (div 4/3 ; Ω)×L 2 skew (Ω). One could also approximate these auxiliary variables using DNNs. However, we restrict our experiments to (u, p) as these are the primary variables of interest in the problem.To conclude this discussion, in Fig. 7 we plot the numerical results for the approximation of the pressure p in the above problem. This complements Fig. 2 , which showed results for the velocity field u. We once more observe similar results: ELU and tanh DNNs outperform ReLU DNNs, the rate of convergence appears to be close to O(m -1 ) and there is no degradation with increasing parametric dimension d.To recap, in our first example, we considered a mixed formulation of a parametric diffusion equation that provably satisfies the (b, ε)-holomorphy assumption. Using this formulation, we considered problems with nonzero Dirichlet boundary conditions, whereas previous works [2, 19, 27] study the more restrictive case of homogeneous Dirichlet boundary conditions, where u ∈ H 1 0 (Ω). In our next example, we studied a more complicated parametric PDE, namely, the parametric NSB equations. While this example currently lacks a holomorphy guarantee, we observe a convergence rate that aligns with what we expect. We conjecture that the m -1 rate holds both for this and for even more complicated problems. To illustrate this claim with an example, we now consider a parametric coupled partial differential equation in three dimensions (Ω ⊂ R 3 ) with two random coefficients affecting different parts of the coupled problem. The nonparametric version of this problem is based on [20] . Specifically, we consider the Boussinesq formulation in [20] that combines a parametric incompressible Navier-Stokes equation with a parametric heat equation. The parametric dependence affects both equations. The Navier-Stokes equation is affected by a parametric variable multiplying the temperature-dependent viscosity, and the equation for heat flow is affected directly by the thermal conductivity of the fluid. To be more precise, given x ∈ [-1, 1] d , our goal is to find the velocity u :[-1, 1] d × Ω → R 2 , pressure p : [-1, 1] d × Ω → R and temperature φ : [-1, 1] d × Ω → R of a fluid such that -div(2a(x)ϖ(φ(x))e(u(x))) + (u(x) • ∇)u(x) + ∇p(x) = φ(x)g in Ω, div(u(x)) = 0 in Ω, (B.16) -div(K(x)∇φ(x)) + u(x) • ∇φ(x) = 0 in Ω, u = u D on ∂Ω, φ = φ D on ∂Ω, Ω p(x) = 0.Here g = (0, 0, -1) ⊤ is a gravitational force andK(x) = K(•; x) ∈ L ∞ (Ω), where K : [-1, 1] d × Ω → R 3×3is a parametric uniformly-positive tensor describing the thermal conductivity of the fluid. It is given explicitly byK(z, x) =   1.89 + j∈N x j sin(πz 3 j) j 9/5   exp(-z 1 ) 0 0 0 exp(-z 2 ) 0 0 0 exp(-z 3 ), ∀z ∈ Ω, (B.17)for x ∈ [-1, 1] d .The term ϖ : R → R + is a temperature-dependent viscosity given by ϖ(φ) = 0.1 + exp(-φ) and the term a(x) = a(•; x) ∈ L ∞ (Ω), where a : [-1, 1] d × Ω → R is a parametric variable affecting the viscosity of the fluid. As in the previous example e(u) is the symmetric part of ∇u. Note that in this case, we have (a(x), K(x)) ∈ X := L ∞ (Ω) × L ∞ (Ω).The complete derivation of a fully-mixed variational formulation for the non-parametric Boussinesq equation in Banach spaces can be found in [20, Sec. 3.1] . To make the presentation simpler we rewrite it for the parametric case. Givenx ∈ [-1, 1] d , find (u, t, σ, φ, t, σ)(x) ∈ L 4 (Ω) × L 2 tr (Ω) × H 0 (div 4/3 ; Ω) × L 4 (Ω) × L 2 (Ω) × H(div 4/3 ; Ω) such that - Ω v • div(σ(x)) + 1 2 Ω t(x)u(x) • v - Ω φ(x)g • v = 0 Ω 2a(x)ϖ(φ(x))t sym (x) : s - 1 2 Ω (u ⊗ u)(x) : s = Ω σ(x) : s Ω τ : t(x) + Ω u(x) • div(τ ) = ⟨τ ν, u D ⟩ Γ - Ω ψ div( σ(x)) + 1 2 Ω ψ(x)u(x) • t = 0 Ω K(x) t(x) • s - 1 2 Ω φ(x)u(x) • s = Ω σ(x) • s Ω τ • t(x) + Ω φ(x) div( τ ) = ⟨ τ • ν, φ D ⟩ Γ Ω tr(2σ + u ⊗ u)(x) = 0, (B.18) for all (v, s, τ , ψ, σ, τ ) ∈ L 4 (Ω) × L 2 tr (Ω) × H 0 (div 4/3 ; Ω) × L 4 (Ω) × L 2 (Ω) × H(div 4/3 ; Ω). Here p ∈ L 2 0 (Ω) can be recovered as p = - 1 6 tr(2σ + 2cI + u ⊗ u), where c = - 1 6|Ω| Ω tr(u ⊗ u). (B.19)As in the previous example, we omitted the term x from this equation for simplicity. For further details on this formulation we refer to [20] and references within. 18 ) by using DNNs and study the approximation capabilities as we increase the number of training samples m. As in the previous example (see Remark B.1), we do not aim to approximate the other variables (t, σ, t, σ)(x) ∈ L 2 tr (Ω) × H 0 (div 4/3 ; Ω) × L 2 (Ω) × H(div 4/3 ; Ω). In our experiments, we consider the unit cube Ω = (0, 1) 3 as the domain in R 3 . We consider a nonzero boundary condition u D = (1, 1, 0) on the bottom face of the cube ∂Ω bottom = (0, 1) × (0, 1) × {0}, and zero on the other faces. We set φ D = exp(4(-(z 1 -0.5) 2 -(z 2 -0.5) 2 )) on ∂Ω bottom and zero otherwise. For simplicity, we consider the same parametric coefficients a 1,d and a 2,d given by (B.1) and (B.2), respectively. See Fig. 8 for an example of the solution (u, p, φ)(x) for a given x ∈ [-1, 1] 4 .Given x ∈ [-1, 1] d , we approximate the solution (u, p, φ)(x) ∈ (L 4 (Ω) × L 2 0 (Ω) × L 4 (Ω)) of (B.To conclude this section, we provide a comparison of the performance of the DNN architectures in approximating the velocity field u and pressure p for the Boussinesq PDE. This complements Fig. 3 , which showed results for the temperature φ. As with φ, the convergence rate for p agrees roughly with the rate m -1 , and does not appear to deteriorate with the parametric dimension d. On the other hand, the convergence rate for the velocity field u is somewhat slower.In this section, we first introduce additional notation that is needed for the proofs of the main results. We then give a brief overview of the proofs.Let (X , ∥ • ∥ X ) and (Y, ∥ • ∥ Y ) be Banach spaces, G : X → Y and B ⊆ X . We define the Lipschitz constant as L = Lip(G; B, Y) as the smallest constant L ≥ 0 such that ∥G(X ′ ) -G(X)∥ Y ≤ L∥X ′ -X∥ X , ∀X, X ′ ∈ B.We require some notation for sequences. Let (Z, ∥ • ∥ Z ) be a Banach space, d = N ∪ {∞} and write ν = (ν k ) d k=1 for an arbitrary multi-index in N d 0 . If Λ ⊆ N d 0 is a finite or countable set of multi-indices and 0 < p ≤ ∞ we define the space ℓ p (Λ; Z) as the set of all Z-valued sequences c = (c ν ) ν∈Λ for which ∥c∥ p;Z < ∞, where∥c∥ p;Z = ν∈Λ ∥c ν ∥ p Z 1/p 0 < p < ∞, sup ν∈Λ ∥c ν ∥ Z p = ∞. When (Z, ∥ • ∥ Z ) = (R, | • |), we just write ℓ p (Λ) and ∥ • ∥ p . We also define the following:|||c||| p;Z = sup z * ∈B(Z * ) ∥z * (c)∥ p , ∀c ∈ ℓ p (Λ; Z), where z * (c) = (z * (c ν )) ν∈Λ ∈ R |Λ| for c = (c ν ) ν∈Λ . Notice that |||c||| p;Z ≤ ∥c∥ p;Z .Given a sequence c = (c ν ) ν∈Λ , we writesupp(c) = {ν ∈ Λ : c ν ̸ = 0} ⊆ Λ.For d ∈ N ∪ {∞}, we write e j , j = 1, . . . , d, for the canonical basis vectors in R d . We also write 0 for the multi-index of zeros and 1 for the multi-index of ones. Next, for a Banach space (Z, ∥ • ∥ Z ) and 0 < p ≤ 2, we define the weighted ℓ p w (Λ; Z) space as the space of Z-valued sequences c = (c ν ) ν∈Λ for which ∥c∥ p,w;Z < ∞, where∥c∥ p,u;Z = ν∈Λ w 2-p ν ∥c ν ∥ p Z 1/p .Notice that ∥ • ∥ 2,u;Z coincides with the unweighted norm ∥ • ∥ 2;Z .We write {P n } n∈N0 for the classical Legendre polynomials on [-1, 1] with normalization P n (1) = 1. Since 1 -1 |P n (x)| 2 dx = (n + 1/2) -1 , we define the orthonormal (with respect to the uniform probability measure) Legendre polynomials asψ n (x) = √ 2n + 1P n (x), x ∈ [-1, 1], n ∈ N 0 . Let D = [-1, 1] N as before, F = ν = (ν i ) ∞ i=1 ∈ N N 0 : |supp(ν)| < ∞ bethe set of multi-indices with finitely-many nonzero terms and define the multivariate Legendre polynomials asΨ ν (x) = i∈N ψ νi (y i ) ≡ i∈supp(ν) ψ νi (x i ), ∀x ∈ D, ν = (ν i ) ∞ i=1 ∈ F.Here, the second equality follows from the fact that ψ 0 = 1. Then it is known that the set{Ψ ν } ν∈F ⊂ L 2 ϱ (D) (C.2)constitutes an orthonormal basis for L 2 ϱ (D) [19, §3] . For later use, we also define the sequencesu = (u ν ) ν∈F , where u ν = ∥Ψ ν ∥ L ∞ ϱ (D) = k∈N √ 2ν k + 1, ∀ν ∈ F, (C.3) and v = (v ν ) ν∈F , where v ν = u 5+ξ ν , ∀ν ∈ F.(C.4) Here ξ ≥ 0 will be chosen later in the proof.Given an optimization problem min t f (t), we say that t is a (σ, τ )-approximate minimizer for some σ ≥ 1 and τ ≥ 0 if f ( t) ≤ σ 2 min t f (t) + τ 2 . (In the main paper, we consider only σ = 1, but for the proofs it is useful to allow σ > 1).Finally, for convenience, given X 1 , . . . , X m ∈ X , we define the semi-norm|||G||| disc,µ = 1 m m i=1 ∥G(X i )∥ 2 Yof an operator G : X → Y and Since our goal in this work is 'agnostic' DNN architectures (i.e., independent of the smoothness of the underlying operator), in step (a) we first define a nonlinear set (D.2) spanned by Legendre polynomials with nonzero indices in certain sets of bounded weighted cardinality. This is effectively a form of sparse polynomial approximation, and the analysis of the resulting learning procedure (D.3) relies heavily on techniques from compressed sensing. In order to bound the encoding-decoding error, we also require several results on Lipschitz continuity (Lemma D.2) and norm equivalences (Lemma D.4) for multivariate polynomials.EQUATIONStep (b) relies on what have now become fairly standard results in DNN approximation theory: namely, the approximate emulation of orthogonal polynomials via DNNs of given width and depth. We present such a result in Lemma D.9, then use this to define the DNN family N in (D. 19 ).We then analyze the DNN training problem (2.5) in step (c). Using the emulation result of step (b), we first show that any approximate minimizer N of (2.5) yields a polynomial that is also an approximate minimizer of the polynomial training problem (D.3) (Lemma D.12). We may then apply the results shown in Step (a) to prove a generalization bound (Theorem D.13). Up to this point, we have not used the holomorphy assumption. We now use this assumption to bound the various best polynomial approximation errors that arise in the previously-derived generalization bound ( §D.7). Finally, in §D.8 we put all these estimates together to complete the proof.Remark C.1 Theorem 3.1 is a generalization and improvement of [5, Thms. 4.1 & 4.2], which deals with the case of learning Banach-valued functions rather than operators. Specifically, the setting of [5] can be considered a special case of this paper where X = ℓ ∞ (N) and the encoding error E X ,q = 0. Moreover, Theorem 3.1 improves the main results of [5] in three key ways. First, the the DNN architectures bounds are much narrower: width(N ) ≲ (m/L) 1+δ versus width(N ) ≲ m 3+log 2 (m) in the latter. Second, Theorem 3.1 considers standard training, i.e., ℓ 2 -loss minimization, whereas [5, Thms. 4.1 & 4.2] requires regularization. Third, for Banach spaces, the error decay rate with respect Both Λ and k will be chosen later in the proof. For any Banach space (Z, ∥ • ∥ Z ) define the spaceP S;Z = ν∈S c ν Ψ ν : c ν ∈ Z, S ∈ S . (D.2)Then, given the training data (2.4), consider the problemmin p∈P S; Y 1 m m i=1 ∥Y i -p • E X (X i )∥ 2 Y , (D.3)whereY = D Y (R d Y ).Here and throughout the proofs, we slightly abuse notation: the polynomial p : R N → R, whereas E X has codomain R d X . However, by construction, p is independent of all but the first d X variables. Hence we may consider p as a function R d X → R. This aside, if p is an approximate minimizer of (D.3), then we define the approximation to F as F ≈ F = p • E X . (D.4) Notice that p • E X = D Y • P • E X , where P : R d X → R d Y isWe require several lemmas. The first relates the L ∞ ϱ -norm of a polynomial to its Pettis L 2 ϱ -norm. Lemma D.1 (Nikolskii inequality for polynomials). Let (Z, ∥ • ∥ Z ) be any Banach space and p = ν∈S c ν Ψ ν for some finite set S ⊂ F, where c ν ∈ Z. Then∥p∥ L ∞ ϱ (D;Z) ≤ |S| u |||p||| L 2 ϱ (D;Z) .Proof. By definition∥p∥ L ∞ ϱ (D;Z) = |||p||| L ∞ ϱ (D;Z) = sup z * ∈B(Z * ) ∥z * (p)∥ L ∞ ϱ (D) .Fix z * ∈ B(Z * ) and write z * (p) = ν∈S z * (c ν )Ψ ν . By the triangle inequality and the definition of the weights (C.3), we have∥z * (p)∥ L ∞ ϱ (D) ≤ ν∈S |z * (c ν )|∥Ψ ν ∥ L ∞ ϱ (D) = ν∈S u ν |z * (c ν )|.We now apply the Cauchy-Schwarz inequality, (C.1) and Parseval's identity to get∥z * (p)∥ L ∞ ϱ (D) ≤ |S| u ν∈S |z * (c ν )| 2 = |S| u ∥z * (p)∥ L 2 ϱ (D) .Since z * was arbitrary, we deduce that∥p∥ L ∞ ϱ (D;Z) ≤ |S| u sup z * ∈B(Z * ) ∥z * (p)∥ L 2 ϱ (D) = |S| u |||p||| L 2 ϱ (D;Z) ,as required.Next, we require the following bound on the Lipschitz constant of a multivariate polynomial. Lemma D.2 (Lipschitz continuity for polynomials). Let (Z, ∥ • ∥ Z ) be any Banach space and suppose that p = ν∈S c ν Ψ ν for some finite S ⊂ F, where c ν ∈ Z. Then satisfiesLip(p; B ∞ (N), Z) ≤ 1 2 |S| v • |||p||| L 2 ϱ (D;Z) , where B ∞ (N) = {x ∈ R N : ∥x∥ ∞ ≤ 1} is the unit ball of ℓ ∞ (N).Proof. Fix z * ∈ B(Z * ) and let p = z * (p) = ν∈S z * (c ν )Ψ ν =: ν∈S cν Ψ ν be the corresponding scalar-valued polynomial. Let x, x ′ ∈ D. Then the mean value theorem gives thatp(x ′ ) -p(x) = ∞ i=1 (x ′ i -x i ) ∂ p ∂x i (tx + (1 -t)x ′ )for some 0 ≤ t ≤ 1. Since B ∞ (N) ≡ D, this and the fact that D is convex give that|p(x ′ ) -p(x)| ≤ ∥x ′ -x∥ ∞ ∞ i=1 ∂ p ∂x i L ∞ ϱ (D).(D.5)We now consider the terms in the sum separately. Using the definition of the Legendre polynomials (see §C.1.4), we have∂ p ∂x i = ν∈S cν u ν P ′ νi (x i ) j̸ =i P νj (x j ).The unnormalized Legendre polynomials satisy1 = P n (1) = ∥P n ∥ L ∞ ([-1,1]) and n(n + 1)/2 = P ′ n (1) = ∥P ′ n ∥ L ∞ ([-1,1]) . Hence ∂ p ∂x i L ∞ ϱ (D) ≤ ν∈S |c ν |u ν ν i (ν i + 1)/2.We deduce that∞ i=1 ∂ p ∂x i L ∞ ϱ (D) ≤ ν∈S |c ν |u ν ∞ i=1 ν i (ν i + 1)/2.Now, by definition of the weights u ν and v ν (see (C.3) and (C.4)), we have(ν i + 1) ≤ j∈N (2ν j + 1) = u 2 ν and ∞ i=1 ν i ≤ j∈N (2ν j + 1) = u 2 ν . Hence ∞ i=1 u ν ν i (ν i + 1) ≤ u 5 ν ≤ v ν .Here we also used the fact that u ≥ 1. We now apply this, the Cauchy-Schwarz inequality and Parseval's identity to get∞ i=1 ∂ p ∂x i L ∞ ϱ (D) ≤ 1 2 ∥p∥ L 2 ϱ (D) |S| v .Substituting this into (D.5) now gives|p(x ′ ) -p(x)| ≤ 1 2 ∥p∥ L 2 ϱ (D) |S| v ∥x ′ -x∥ ∞ .We now recall that p = z * (p) and z * ∈ B(Z * ) was arbitrary to get∥p(x ′ ) -p(x)∥ Z = sup z * ∈B(Z * ) |z * (p(x ′ ) -p(x))| ≤ 1 2 sup z * ∈B(Z * ) ∥z * (p)∥ L 2 ϱ (D) |S| v ∥x ′ -x∥ ∞ = 1 2 |||p||| L 2 ϱ (D;Z) |S| v ∥x ′ -x∥ ∞ ,as required.We now show that this result can be used to imply a norm equivalence for polynomials. For this we first require the following lemma. Note that in this and subsequent results, we abuse notation and write ς for both the measure on [-1, 1] d X defined in (A.III) and the measure on D = [-1, 1] N defined by tensoring this measure (corresponding to the first d X variables x 1 , . . . , x d X ) with the uniform measure on D\[-1, 1] d X (corresponding to the remaining variables x d X +1 , x d X +2 , . . .). Lemma D.3 (Closeness of L 2 norms). Let (Z, ∥ • ∥ Z ) be any Banach space, ς be the measure defined in (A.I) and ς be the measure defined in (A.III). Suppose that f ∈ L 2 ς (D; Z) is Lipschitz continuous with constant L = Lip(f ; B ∞ (N), Z) < ∞ and that f depends only on its first d X variables. Thenf ∈ L 2 ς (D; Z) and |||f ||| L 2 ς (D;Z) -δ ≤ |||f ||| L 2 ς (D;Z) ≤ |||f ||| L 2 ς (D;Z) + δ, where δ = L • ∥ι d X -ι d X • D X • E X ∥ L 2 µ (X ;ℓ ∞ (N)) . Proof. Fix z * ∈ B(Z * ), let g = z * (f ) ∈ L 2 ς (D) and set a = ∥g∥ L 2 ς (D) , b = ∥g∥ L 2 ς (D) . Notice that |g(x) -g(x ′ )| ≤ ∥z * ∥ Z * ∥f (x) -f (x ′ )∥ Z ≤ L∥x -x ′ ∥ ∞ , ∀x, x ′ ∈ B ∞ (N).Now, with slight abuse of notation, g(ι(X)) = g(ι d X (X)) due to the assumption on f . Therefore, using this and the Cauchy-Schwarz inequality,a 2 -b 2 = D (g(ι d X • D X • E X (X))) 2 -(g(ι d X (X))) 2 dµ(X) ≤ L D ∥ι d X (X) -ι d X • D X • E X (X)∥ ∞ |g(ι d X (X))| + |g(ι d X • D X • E X )| dµ(X) ≤ L∥ι d X -ι d X • D X • E X ∥ L 2 µ (X ;ℓ ∞ (N)) (a + b) .We deduce that a -b ≤ δ. Since z * was arbitrary, we get|||f ||| L 2 ς (D;Z) = sup z * ∈B(Z * ) ∥z * (f )∥ L 2 ς (D) ≤ sup z * ∈B(Z * ) ∥z * (f )∥ L 2 ς (D) + δ = |||f ||| L 2 ς (D;Z) + δ,which gives the upper bound. The same argument applied to b 2 -a 2 also gives the lower bound.Lemma D.4 (Norm equivalences for polynomials). Let (Z, ∥ • ∥ Z ) be any Banach space, ς be the measure defined in (A.I), ς be the measure defined in (A.III) and S ⊂ F with supp(ν) ∈ {1, . . . , d X }, ∀ν ∈ S. Suppose that|S| v • ∥ι d X -ι d X • D X • E X ∥ L 2 µ (X ;ℓ ∞ (N)) ≤ c, (D.6)for some sufficiently small universal constant c > 0. Then the norm equivalence|||p||| L 2 ϱ (D;Z) ≲ |||p||| L 2 ς (D;Z) ≲ |||p||| L 2 ϱ (D;Z)holds for all p = ν∈S c ν Ψ ν , where c ν ∈ Z.Proof. Combining Lemmas D.2 and D.3, we see that|||p||| L 2 ς (D;Z) -δ ≤ |||p||| L 2 ς (D;Z) ≤ |||p||| L 2 ς (D;Z) + δ, where δ = 1 2 |||p||| L 2 ϱ (D;Z) |S| v ∥ι d X -ι d X • D X • E X ∥ L 2 µ (X ;ℓ ∞ (N)) ≤ c|||p||| L 2 ϱ (D;Z) /2. By (A.I), there are constants c 1 ≥ c 2 > 0 such that c 1 |||p||| L 2 ϱ (D;Z) ≤ |||p||| L 2 ς (D;Z) ≤ c 2 |||p||| L 2 ϱ (D;Z) . We now take c = c 1 /2.We now analyze (D.3). Our analysis relies on the following result, which shows an error bound subject to a certain discrete metric inequality (D.7).Q : Y → Y = D Y (R d Y ) be a bounded linear operator and π Q = ∥Q∥ Y→Y . Suppose that ∥p -q∥ disc,ς ≥ α max{|||p -q||| L 2 ϱ (D;Y) , |||p -q||| L 2 ς (D;Y) }, ∀p, q ∈ P S; Y (D.7)for some α > 0. Then, for any F ∈ L ∞ µ (X ; Y), p ∈ P S; Y and q ∈ P S;Y , we have|||F -p • E X ||| L 2 µ (X ;Y) ≤ |||F -Q • F ||| L 2 µ (X ;Y) + α -1 ∥p -Q • q∥ disc,ς + π Q |||F -q • E X ||| L 2 µ (X ;Y) and |||F -p • E X ||| L ∞ µ (X ;Y) ≤ |||F -Q • F ||| L ∞ µ (X ;Y) + √ 2k/α∥p -Q • q∥ disc,ς + π Q ∥F -q • E X ∥ L ∞ µ (X ;Y) .Proof. By the triangle inequality and properties of Q, we have|||F -p • E X ||| L 2 µ (X ;Y) ≤ |||F -Q • F ||| L 2 µ (X ;Y) + |||p • E X -Q • q • E X ||| L 2 µ (X ;Y) + |||Q • F -Q • q • E X ||| L 2 µ (X ;Y) ≤ |||F -Q • F ||| L 2 µ (X ;Y) + |||p • E X -Q • q • E X ||| L 2 µ (X ;Y) + π Q |||F -q • E X ||| L 2 µ (X ;Y) . Now, since p, Q • q ∈ P S; Y , the second term can be bounded by |||p • E X -Q • q • E X ||| L 2 µ (X ;Y) = |||p -Q • q||| L 2 ς (D;Y) ≤ α -1 ∥p -Q • q∥ disc,ς . (D.8)This yields the first result.For the second result, we once more write|||F -p • E X ||| L ∞ µ (X ;Y) ≤ |||F -Q • F ||| L ∞ µ (X ;Y) + |||p • E X -Q • q • E X ||| L ∞ µ (X ;Y) + π Q |||F -q • E X ||| L ∞ µ (X ;Y). For the second term, we use (A.III) to write|||p • E X -Q • q • E X ||| L ∞ µ (X ;Y) = |||p -Q • q||| L ∞ ς (D;Y) ≤ ∥p -Q • q∥ L ∞ ϱ (D;Y) . Notice that p -Q • q is apolynomial supported in a set S with |S| u ≤ 2k. We now apply Lemma D.1 and (D.7) to obtain|||p • E X -Q • q • E X ||| L ∞ µ (X ;Y) ≤ √ 2k|||p -Q • q||| L 2 ϱ (D;Y) ≤ √ 2k/α|||p -Q • q||| disc,ς . (D.9)This gives the result.Q : Y → Y = D Y (R d Y )be a bounded linear operator, π Q = ∥Q∥ Y→Y and suppose that (D.7) holds. Let F be as in (D.4) for some (σ, τ )-approximate minimizer p of (D.3). Then|||F -F ||| L 2 µ (X ;Y) ≤ |||F -Q • F ||| L 2 µ (X ;Y) + σ + 1 α |||F -Q • F ||| disc,µ + π Q |||F -q • E X ||| L 2 µ (X ;Y) + σ + 1 α |||F -q • E X ||| disc,µ + τ α + σ + 1 α √ m ∥E∥ 2;Y and ∥F -F ∥ L ∞ µ (X ;Y) ≤ ∥F -Q • F ∥ L ∞ µ (X ;Y) + √ 2k(σ + 1) α |||F -Q • F ||| disc,µ + π Q ∥F -q • E X ∥ L ∞ µ (X ;Y) + √ 2k(σ + 1) α |||F -q • E X ||| disc,µ + √ 2kτ α + √ 2k(σ + 1) α √ m ∥E∥ 2;Yfor all q ∈ P S;Y , where E = (E i ) m i=1 ∈ Y m is the (Banach-valued) vector of noise terms.Proof. We apply the previous lemma with p = p. This gives|||F -F ||| L 2 µ (X ;Y) ≤ |||F -Q • F ||| L 2 µ (X ;Y) + α -1 ∥p -Q • q∥ disc,ς + π Q |||F -q • E X ||| L 2 µ (X ;Y) ∥F -F ∥ L ∞ µ (X ;Y) ≤ ∥F -Q • F ∥ L ∞ µ (X ;Y) + √ 2k/α∥p -Q • q∥ disc,ς + π Q ∥F -q • E X ∥ L ∞ µ (X ;Y) . (D.10)Consider the second term. We have∥p -Q • q∥ disc,ς = |||p • E X -Q • q • E X ||| disc,µ ≤ |||F -Q • q • E X ||| disc,µ + |||F -p • E X ||| disc,µ .Consider the second term of this expression. By the triangle inequality and the facts that p is a (σ, τ )-minimizer and Q • q is feasible, we obtain|||F -p • E X ||| disc,µ ≤ 1 m m i=1 ∥Y i -p • E X ∥ 2 Y + 1 √ m ∥E∥ 2;Y ≤ σ 1 m m i=1 ∥Y i -Q • q • E X ∥ 2 Y + τ + 1 √ m ∥E∥ 2;Y ≤ σ|||F -Q • q • E X ||| disc,µ + τ + σ + 1 √ m ∥E∥ 2;Y Therefore, we get ∥p -Q • q∥ disc,ς ≤ (σ + 1)|||F -Q • q • E X ||| disc,µ + τ + σ + 1 √ m ∥E∥ 2;Y .We now estimate the first term in this expression as follows:|||F -Q • q • E X ||| disc,µ ≤ |||F -Q • F ||| disc,µ + |||Q • F -Q • q • E X ||| disc,µ ≤ |||F -Q • F ||| disc,µ + π Q |||F -q • E X ||| disc,µ .Therefore, we conclude that∥p -Q • q∥ disc,ς ≤ (σ+1)|||F -Q • F ||| disc,µ +(σ+1)π Q |||F -q • E X ||| disc,µ +τ + σ + 1 √ m ∥E∥ 2;Y .Combining this with (D.10) now gives the result.For the proof of the next lemma and subsequent steps of the proof, we let Λ = {ν 1 , . . . , ν N } and define the matrixA = Ψ νj (E X (X i )) √ m m,N i,j=1 ∈ R m×N . Lemma D.7. Let 0 < ϵ < 1, 0 < δ < 1, k > 0 and suppose that m ≥ c 0 • δ -2 • k • (log(eN ) • log 2 (k/δ) + log(2/ϵ)) (D.11)for some universal constant c 0 > 0. Let T = {c ∈ R N : ∥c∥ 2 = 1, ∥c∥ 0,v ≤ k} andθ + = sup E∥Ac∥ 2 2 : c ∈ T , θ -= inf E∥Ac∥ 2 2 : c ∈ T .Then ∥Ac∥2 2 ≥ (θ --(1 + θ + )c 1 δ) ∥c∥ 2 2 , ∀c ∈ R N , ∥c∥ 0,v ≤ k. with probability at least 1 -ϵ, where c 1 > 0 is a universal constant.Proof. The proof is based on [12, Thm. 2.13 ]. Since ∥c∥ 1,v ≤ ∥c∥ 0,v ∥c∥ 2 ≤ √ k, we see thatT ⊆ c ∈ R N : ∥c∥ 1,v ≤ √ k . Define the random vector X ∈ R N by X = (Ψ νi (E X (X))) N i=1 for X ∼ µ. Observe that |⟨X, e j ⟩| = |Ψ νi (E X (X))| ≤ u νi ≤ v νialmost surely, by (A.III) and the definition of the weights. Therefore, by [12, Thm. 2.13] , ifEQUATION2 ).Therefore, we have shown that∥Ac∥ 2 2 ≥ θ --(1 + θ + )c 1 δ, ∀c ∈ T, with probability at least 1 -2 exp(-c 2 δ 2 m/k), provided m satisfies (D.12). To conclude the result, we observe that (D.11) implies (D.12), up to a possible change in the universal constant c 0 . Moreover, it also implies that 2 exp(-c 2 δ 2 m/k) ≤ ϵ.Hence we obtain the result.Lemma D.8. There exist universal constants c 0 , c 1 , c 2 > 0 such that the following holds. Suppose thatm ≥ c 0 • k • (log(eN ) • log 2 (k) + log(2/ϵ)) (D.13) and √ k • ∥ι d X -ι d X • D X • E X ∥ L 2 µ (X ;ℓ ∞ (N)) ≤ c 1 . (D.14)Then (D.7) holds with probability at least 1 -ϵ and constant α ≥ c 2 .Proof. We shall apply the previous lemma with k replaced by 2k. First, we estimate θ + and θ -. Let c ∈ T , with T as defined therein with 2k in place of k. Write p = ν∈Λ c ν Ψ ν for the corresponding (scalar-valued) polynomial. Then∥Ac∥ 2 2 = 1 m m i=1 |p • E X (X i )| 2 ,and therefore E∥Ac∥2 2 = ∥p∥ 2 L 2 ς (D) .Since ∥c∥ 0,v ≤ 2k, (D.14) implies (D.6). We now apply Lemma D.4 and the fact that ∥p∥L 2 ϱ (D) = ∥c∥ 2 = 1 to get c 3 ≤ E∥Ac∥ 2 2 ≤ c 4 for universal constants c 4 ≥ c 3 > 0. Since c ∈ T was arbitrary to deduce that c 3 ≤ θ -≤ θ + ≤ c 4 .(D.15)We now show that (D.7) holds with the desired probability. Let p, q ∈ P S; Y be arbitrary. Then their difference h = p -q can be expressed as 16 )h = ν∈Λ d ν Ψ ν (D.where d = (d ν ) ν∈Λ ∈ Y N satisfies ∥d∥ 0,u ≤ ∥d∥ 0,v ≤ 2k.Now, this, (D.14) and Lemma D.4 imply that|||h||| L 2 ς (D;Y) ≤ c 4 |||h||| L 2 ϱ (D;Y) . Therefore, in order to prove (D.7), it suffices to show that |||h||| L 2 ϱ (D;Y) ≲ ∥h∥ disc,ς . Observe that ∥h∥ disc,ς = ∥Ad∥ 2;Y ≥ |||Ad||| 2;Y = sup y * ∈B(Y * ) ∥y * (Ad)∥ 2 ,where we recall that for a vector z= (z i ) N i=1 ∈ Y N , we write y * (z) = (y * (z i )) N i=1 ∈ R N .By linearity, we have y * (Ad) = Ay * (d). Therefore, by Lemma D.7,∥h∥ 2 disc,ς = ∥Ad∥ 2 2;Y ≥ sup y * ∈B(Y * ) ∥Ay * (d)∥ 2 2 ≥ sup y * ∈B(Y * ) (θ --(1 + θ + )c 5 δ) ∥y * (d)∥ 2 2 ≥ (c 3 -(1 + c 4 )c 5 δ) sup y * ∈B(Y * ) ∥y * (h)∥ 2 L 2 ϱ (D) = (c 3 -(1 + c 4 )c 5 δ) |||h||| 2 L 2ϱ (D;Y) . for some universal constant c 5 > 0, provided m satisfies (D.11). We now set δ = c 3 /(2(1+ c 4 )c 5 ) to get ∥h∥ 2 disc,ς ≥ c 3 /2|||h||| 2 L 2ϱ (D;Y) , provided m satisfies (D.11) with this value of δ. However, this is implied by the condition (D.13). We deduce the result.Recall that Λ ⊂ F, |Λ| = N is an arbitrary set and S is defined by (D.1). In this and what follows, we slightly abuse notation and consider a DNN N : R d → R as a function R N → R which depends on only the first d variables. Lemma D.9 (Approximating Legendre polynomials with tanh DNNs). Let Γ ⊂ Λ with supp(ν) ⊆ {1, . . . , d}, ∀ν ∈ Γ, and m(Γ) = max ν∈Γ ∥ν∥ 1 < ∞. There exists a fully-connected family N o of tanh DNNs N : R d → R withwidth(N o ) ≲ m(Γ), depth(N o ) ≲ log(m(Γ)),such that, for any 0 < δ < 1 and ν ∈ Γ, there is a DNNN ν ∈ N o with ∥N ν -Ψ ν ∥ L ∞ ϱ (D) ≤ δ. (D.17)Moreover, the zero network 0 : x → 0 also belongs to N o (trivially, since tanh(0) = 0).Proof. We follow the argument of [5, Thm. 7.4] , which is based on [79, Prop. 2.6]. Let ν ∈ Γ. Then via the fundamental theorem of algebra we can write Ψ ν (x) as a product of ∥ν∥ 1 numbers asΨ ν (x) = i∈supp(ν) νi j=1 d i (x i -r ij ).Here {r ij } νi j=1 are the roots of the univariate Legendre polynomial P νi . We append m(Γ) -∥ν∥ 1 ones and write Ψ ν (x) as a product of exactly m(Γ) numbers. Now define the affine mapA ν : R d → R m(Γ) , (D.18)so that A ν (x) is the vector consisting of the values d i (x i -r ij ) for j = 1, . . . , ν i and i ∈ supp(ν) and 1 otherwise. To complete the proof, we need to construct a tanh DNN mapping R m(Γ) → R that approximately multiplies these numbers. To do this, we argue as in the proof of [5, Thm. 7.4] to see that there is a tanh DNN N ν withwidth(N ν ) ≲ m(Γ), depth(N ν ) ≲ log(m(Γ))that satisfies the desired bound (D.17). Since these width and depth bounds are independent of ν, we deduce the result.Fix δ > 0, let Γ = ∪ S∈S S, d = d X and consider the corresponding family N o and DNNs N ν , ν ∈ Γ, asserted by this lemma. For any ν ∈ Γ, we have ν ∈ S for some S ∈ S, and any such S satisfies |S| v ≤ k. Therefore u 2(5+ξ) ν = v 2 ν ≤ k for any ν ∈ S and we deduce that∥ν∥ 1 ≤ d j=1 (2ν j + 1) = u 2 ν ≤ k 1/(5+ξ) .This implies that m(Γ) ≤ k 1/(5+ξ) and thereforewidth(N o ) ≲ k 1/(5+ξ) , depth(N o ) ≲ log(k).With this in hand, we now define the family N of DNNs N : Rd X → R d Y by N =                  N = C          N ν1 . . . N ν |S| 0 . . . 0          : S = {ν 1 , . . . , ν |S| } ∈ S, C ∈ R d Y ×⌊k⌋                  . (D.19)Here we also use the fact that |S| ≤ |S| v ≤ k. Notice that this family satisfieswidth(N ) ≤ k 1+1/(5+ξ) , depth(N ) ≲ log(k), (D.20) due to the bounds for N o . Now let N ∈ N and write C = [c 1 | • • • |c ⌊k⌋ ], where c i ∈ R d Y . Then D Y • N = |S| i=1 D Y (c i )N νi = |S| i=1 c i N νi , where c i ∈ Y.Therefore, we can associate N with the space of functions P S; Y , where, for any arbitrary Banach space (Z, ∥ • ∥ Z ),P S;Z = ν∈S c ν N ν : c ν ∈ Z, S ∈ S .We now require the following lemma, which relates the distance between a function in P S;Z and the corresponding polynomial in P S;Z .Lemma D.10 (Discrete norms of polynomials and their approximating DNNs). Let S ⊂ F, (Z, ∥ • ∥ Z ) be any Banach space and suppose that p = ν∈S c ν Ψ ν , where c ν ∈ Z. Definep = ν∈S c ν N ν . Then ∥p -p∥ disc,ς ≤ ∥p -p∥ L ∞ ϱ (D;Z) ≤ δ |S||||p||| L 2 ϱ (D;Z) . (D.21)Moreover, if Z = Y, S ∈ S and (D.7) holds with α satisfying δ |S| u /α < 1 then|||p||| L 2 ϱ (D;Y) ≤ 1 α -δ |S| ∥p∥ disc,ς .Proof. By (A.III) and the definition of the N ν we have ∥p -p∥ disc,ς ≤ ∥p -p∥ L ∞ ϱ (D;Z) = sup z * ∈B(Z * ) ∥z * (p -p)∥ L ∞ ϱ (D) ≤ sup z * ∈B(Z * ) ν∈S |z * (c ν )|∥Ψ ν -N ν ∥ L ∞ ϱ (D)≤ δ|||c||| 1;Z .We now apply the Cauchy-Schwarz inequality and Parseval's identity to obtain∥p -p∥ disc,ς ≤ δ |S||||c||| 2;Z = δ |S| u |||p||| L 2 ϱ (D;Z), which gives the first result.For the second result, we apply (D.7) with q = 0 to get 21 ) and the fact that δ |S|/α < 1 now gives the result.|||p||| L 2 ϱ (D;Y) ≤ α -1 ∥p∥ disc,ς ≤ α -1 ∥p -p∥ disc,ς + ∥p∥ disc,ς . Using (D.Remark D.11 (Other activation functions) As seen in this section, a key step in our proofs is emulating the polynomials via DNNs of quantifiable width and depth. There is an extensive literature on this topic. See, e.g., [5, 21, 23, 25, 38, 57, 63, 67, 72, 78, 79, 85, 86, 92, 97] and references therein. The proof of Lemma D.9 reduces this to the task of emulating the multiplication operation(x 1 , . . . , x d ) ∈ R d → x 1 • • • x d ∈ R via a DNN.As shown in the proof of [5, Thm. 7.4] (which is based on [79, Prop. 2.6]), this can in turn be achieved using a binary tree of ⌈log 2 (d)⌉ DNNs that approximately compute the multiplication of two numbers (x, y) ∈ R 2 → xy ∈ R. Further, this task can be achieved via the identity xy = ((x + y) 2 -(x -y) 2 )/4 by using a DNN that approximately computes the squaring function x ∈ R → x 2 ∈ R. To summarize, provide a DNN of quantifiable width and depth can approximately compute the squaring function, it can also approximately emulate the multivariate Legendre polynomials.In view of this, we can adapt our main theorems to various other activation functions without change. This includes Rectified Polynomial Units (RePUs), where the emulation is, in fact, exact (see, e.g., [57, Lem. 2.1] ). It also includes the Exponential Linear Unit (ELU) used in our numerical experiments and many others. See, for instance, Proposition 4.7 of [38] and the ensuing discussion. Rectified Linear Units (ReLUs) are slightly different, as in this case the depth of the DNN that performs the approximation multiplication depends on the desired accuracy (see, e.g., [79, Prop. 2.6] ). One could modify Theorem 3.1 to consider ReLU DNNs, with the result being a worse depth bound than that presented for tanh DNNs. D.6 Analysis of (2.5) Lemma D.12 (Approximate minimizers of (2.5) yield approximate minimizers of (D.3)). Suppose that (D.7) holds, let N be the family of DNNs defined in §D.5 and N be any (σ, τ )-approximate minimizer of (2.5). LetD Y • N = ν∈S ĉν N ν ∈ P S; Y ,where S ∈ S and ĉν ∈ Y, and definep = ν∈S ĉν Ψ ν .Then p is a (σ ′ , τ ′ )-approximate minimizer of (D.3), whereσ ′ ≤ σ(1 + δ √ k/α) τ ′ ≤ τ + σδ √ k/α ∥F ∥ L ∞ µ (X ;Y) + 1 √ m ∥E∥ 2;Y + δ √ k|||p||| L 2 ϱ (D;Y) . (D.22)Proof. Let p = ν∈S c ν Ψ ν ∈ P S; Y be arbitrary and N ∈ N be the corresponding DNN so thatD Y • N = ν∈S c ν N ν .Then by the triangle inequality and the fact that N is an approximate minimizer, we have1 m m i=1 ∥Y i -p • E X (X i )∥ 2 Y ≤ 1 m m i=1 ∥Y i -D Y • N • E X (X i )∥ 2 Y + ∥p -D Y • N ∥ disc,ς ≤ σ 1 m m i=1 ∥Y i -D Y • N • E X (X i )∥ 2 Y + τ + ∥p -D Y • N ∥ disc,ς ≤ σ 1 m m i=1 ∥Y i -p • E X (X i )∥ 2 Y + τ + σ∥p -D Y • N ∥ disc,ς + ∥p -D Y • N ∥ disc,ς .We now apply Lemma D.10 to the last two terms, noting that |S| ≤ |S| v ≤ k since S ∈ S, to obtain1 m m i=1 ∥Y i -p • E X (X i )∥ 2 Y ≤ σ 1 m m i=1 ∥Y i -p • E X (X i )∥ 2 Y + τ + σδ √ k|||p||| L 2 ϱ (D;Y) + δ √ k|||p||| L 2 ϱ (D;Y). Consider the third term. We first apply (D.7) with q = 0 to get|||p||| L 2 ϱ (D;Y) ≤ α -1 ∥p∥ disc,ς = α -1 |||p • E X ||| disc,µ .We then use the triangle inequality to get|||p • E X ||| disc,µ ≤ |||F -p • E X ||| disc,µ + |||F ||| disc,µ ≤ 1 m m i=1 ∥Y i -p • E X (X i )∥ 2 Y + ∥F ∥ L ∞ µ (X ;Y) + 1 √ m ∥E∥ 2;Y .Therefore, we obtain1 m m i=1 ∥Y i -p • E X (X i )∥ 2 Y ≤ σ(1 + δ √ k/α) 1 m m i=1 ∥Y i -p • E X (X i )∥ 2 Y + τ + σδ √ k/α ∥F ∥ L ∞ µ (X ;Y) + 1 √ m ∥E∥ 2;Y + δ √ k|||p||| L 2 ϱ (D;Y) .Since p ∈ P S; Y was arbitrary we get the result.Q : Y → Y = D Y (R d Y )be a bounded linear operator, π Q = ∥Q∥ Y→Y and suppose that (D.7) holds with α ≥ c 0 and α -δ √ k ≥ c 1 for suitable universal constants c 0 , c 1 > 0. Let N be the family of DNNs defined in §D.5 and N be any (σ, τ )-approximate minimizer of (2.5). Then the approximationF = D Y • N • E X satisfies |||F -F ||| L 2 µ (X ;Y) ≲ |||F -Q • F ||| L 2 µ (X ;Y) + σ|||F -Q • F ||| disc,µ + π Q |||F -q • E X ||| L 2 µ (X ;Y) + σ|||F -q • E X ||| disc,µ + τ + σδ √ k∥F ∥ L ∞ µ (X ;Y) + σ √ m ∥E∥ 2;Y and ∥F -F ∥ L ∞ µ (X ;Y) ≲ ∥F -Q • F ∥ L ∞ µ (X ;Y) + √ kσ|||F -Q • F ||| disc,µ + π Q ∥F -q • E X ∥ L ∞ µ (X ;Y) + √ kσ|||F -q • E X ||| disc,µ + √ kτ + σδk∥F ∥ L ∞ µ (X ;Y) + √ kσ √ m ∥E∥ 2;Yfor all q ∈ P S;Y .Proof. Let p be as in Lemma D.12. Then|||F -F ||| L 2 µ (X ;Y) ≤ |||F -p • E X ||| L 2 µ (X ;Y) + |||p • E X -D Y • N • E X ||| L 2 µ (X ;Y) ∥F -F ∥ L ∞ µ (X ;Y) ≤ ∥F -p • E X ∥ L ∞ µ (X ;Y) + ∥p • E X -D Y • N • E X ∥ L ∞ µ (X ;Y) . (D.23)For the second term, we have, by (A.III) and Lemma D.10,|||p • E X -D Y • N • E X ||| L 2 µ (X ;Y) ≤ ∥p • E X -D Y • N • E X ∥ L ∞ µ (X ;Y) = ∥p -D Y • N ∥ L ∞ ς (D;Y) ≤ ∥p -D Y • N ∥ L ∞ ϱ (D;Y) ≤ δ √ k|||p||| L 2 ϱ (D;Y). We now apply this, Lemma D.12, Theorem D.6 and the facts that α -1 ≲ 1 and σ ′ ≥ 1 to obtain|||F -F ||| L 2 µ (X ;Y) ≲ |||F -Q • F ||| L 2 µ (X ;Y) + σ ′ |||F -Q • F ||| disc,µ + π Q |||F -q • E X ||| L 2 µ (X ;Y) + σ ′ |||F -q • E X ||| disc,µ + τ ′ + σ ′ √ m ∥E∥ 2;Y and ∥F -F ∥ L ∞ µ (X ;Y) ≲ ∥F -Q • F ∥ L ∞ µ (X ;Y) + √ kσ ′ |||F -Q • F ||| disc,µ + π Q ∥F -q • E X ∥ L ∞ µ (X ;Y) + √ kσ ′ |||F -q • E X ||| disc,µ + √ kτ ′ + √ kσ ′ √ m ∥E∥ 2;Yfor any q ∈ P S;Y , where σ ′ and τ ′ are as in (D.22) . Due to the various assumptions, these satisfyσ ′ ≲ σ, τ ′ ≲ τ + σδ √ k∥F ∥ L ∞ µ (X ;Y) + σ √ m ∥E∥ 2;Y + δ √ k|||p||| L 2 ϱ (D;Y) . We deduce that |||F -F ||| L 2 µ (X ;Y) ≲ |||F -Q • F ||| L 2 µ (X ;Y) + σ|||F -Q • F ||| disc,µ + π Q |||F -q • E X ||| L 2 µ (X ;Y) + σ|||F -q • E X ||| disc,µ + τ + σδ √ k∥F ∥ L ∞ µ (X ;Y) + σ √ m ∥E∥ 2;Y + δ √ k|||p||| L 2 ϱ (D;Y) and ∥F -F ∥ L ∞ µ (X ;Y) ≲ ∥F -Q • F ∥ L ∞ µ (X ;Y) + √ kσ|||F -Q • F ||| disc,µ + π Q ∥F -q • E X ∥ L ∞ µ (X ;Y) + √ kσ|||F -q • E X ||| disc,µ + √ kτ + σδk∥F ∥ L ∞ µ (X ;Y) + √ kσ √ m ∥E∥ 2;Y + δk|||p||| L 2 ϱ (D;Y) .We now bound the final term. Using Lemma D.10 once more, in combination with the fact that α -δ √ k ≳ 1, we see that|||p||| L 2 ϱ (D;Y) ≤ 1 α -δ √ k ∥D Y • N ∥ disc,ς ≲ ∥D Y • N ∥ disc,ς . ≤ 1 m m i=1 ∥Y i -D Y • N • E X (X i )∥ 2 Y + 1 √ m ∥Y ∥ 2;Y ≤ σ 1 m m i=1 ∥Y i -0∥ 2 Y + τ + 1 √ m ∥Y ∥ 2;Y ≤ (1 + σ) ∥F ∥ L ∞ µ (X ;Y) + 1 √ m ∥E∥ 2;Y + τ.Here, we also used the fact that N is an approximate minimizer in the fourth step, as well as the facts that the zero network 0 ∈ N and that D Y is a linear map. Plugging this into the previous expressions now gives the result.When σ = 1 (as will be the case when we come to prove Theorem 3.1), the error bounds in Theorem D.13 involve best polynomial approximation error terms of the form|||F -q • E X ||| L 2 µ (X ;Y) + |||F -q • E X ||| disc,µ , ∥F -q • E X ∥ L ∞ µ (X ;Y) + √ k|||F -q • E X ||| disc,µfor arbitrary q ∈ P S;Y . By the triangle inequality, these are bounded byE 2 (F, q) := |||F -q • ι||| L 2 µ (X ;Y) + |||F -q • ι||| disc,µ + |||q • ι -q • E X ||| L 2 µ (X ;Y) + |||q • ι -q • E X ||| disc,µ , E ∞ (F, q) := ∥F -q • ι∥ L ∞ µ (X ;Y) + √ k|||F -q • ι||| disc,µ + ∥q • ι -q • E X ∥ L ∞ µ (X ;Y) + √ k|||q • ι -q • E X ||| disc,µ . (D.24)In this section, we construct a suitable polynomial q in the case where (A.II) holds and thereby derive a bound for these term. We first require the following lemma. Lemma D.14. Let G ∈ L ∞ µ (X ; Y) and 0 < ϵ < 1. Then the following hold.(a) With probability at least 1 -ϵ on the draw of the X i , we have|||G||| disc,µ ≤ ∥G∥ L 2 µ (X ;Y) / √ ϵ.(b) Suppose that m ≥ 2r log(2/ϵ) for some r > 0. Then, with probability at least 1 -ϵ on the draw of the X i , we have|||G||| disc,µ ≤ √ 2 ∥G∥ L ∞ µ (X ;Y) / √ r + ∥G∥ L 2 µ (X ;Y) .Proof. Observe that the random variable|||G||| 2 disc,µ satisfies E|||G||| 2 disc,µ = ∥G∥ 2 L 2 µ (X ;Y) .For (a), we use Markov's inequality to getP |||G||| disc,µ ≥ ∥G∥ L 2 µ (X ;Y) / √ ϵ ≤ E|||G||| 2 disc,µ ∥G∥ 2 L 2 µ (X ;Y) /ϵ = ϵ,as required.The proof of (b) is based on [3, Lem. 7.11 ]. We repeat it here for convenience. Define the random variable Z i = ∥G(X i )∥Y and observe thatE(Z i ) = E X∼µ ∥G(X)∥ 2 Y = ∥G∥ 2 L 2 µ (X ;Y) =: a. Let X i = Z i -E(Z i ) so that |||G||| 2 disc,µ = 1 m m i=1 Z i = 1 m m i=1 X i + a. Now let b = ∥G∥ 2 L ∞ µ (X ;Y) and observe that X i ≤ Z i ≤ b, -X i ≤ E(Z i ) ≤ b, a.s..We also havem i=1 E(X 2 i ) ≤ m i=1 E(Z 2 i ) ≤ b m i=1 E(Z i ) = abm.Therefore, Bernstein's inequality for bounded random variables (see, e.g., [29, Cor. 7.31] ) implies thatP 1 m m i=1 X i ≥ t ≤ 2 exp - t 2 m/2 ab + bt/3for any t > 0. We now set t = a + b/r and notice that t 2 m/2 ab+bt/3 ≥ 3m 5r ≥ log(2/ϵ). Therefore,1 m m i=1 X i < a + b/rwith probability at least 1 -ϵ. It follows that|||G||| disc,µ ≤ 2a + b/r ≤ √ 2 √ a + b/rwith the same probability. Substituting the values for a and b now gives the result.Lemma D.15. Let F ∈ L ∞ µ (X ; Y), q ∈ P S;Y be arbitrary and m ≥ 2r log(6/ϵ) for some r > 0 and 0 < ϵ < 1. ThenE 2 (F, q) ≲ ∥F -q • ι∥ L ∞ µ (X ;Y) / √ r + ∥F -q • ι∥ L 2 µ (X ;Y) + k/ϵ|||q||| L 2 ϱ (D;Y) ∥ι d X -ι d X • D X • E X ∥ L 2 µ (X ;ℓ ∞ (R d X )). andE ∞ (F, q) ≲(1 + k/r)∥F -q • ι∥ L ∞ µ (X ;Y) + √ k∥F -q • ι∥ L 2 µ (X ;Y) + k|||q||| L 2 ϱ (D;Y) ∥ι d X -ι d X • D X • E X ∥ L 2 µ (X ;ℓ ∞ (R d X )) + √ k(1 + k/r)|||q||| L 2 ϱ (D;Y) ∥ι d X -ι d X • D X • E X ∥ L ∞ µ (X ;ℓ ∞ (R d X ))with probability at least 1 -ϵ.Using (C.3)-(C.4), we see thatv ν = k∈supp(ν) (2ν k + 1) (5+ξ)/2and therefored ν ≤ k∈I(ν,ρ) ρ -ν k +1 k (ρ k -1) 2 (2ν k + 1) γfor all ρ satisfying (D.25) and some γ = γ(p, ξ) ≥ 0. We now use [4, Lem. 5.4] to deduce that d = (d ν ) ν∈F ∈ ℓ p (N) with ∥d∥ p ≤ C(b, p, ξ) for some C(b, p, ξ) ≥ 0 depending on b, p and ξ only. Returning to c, this gives∥c∥ p p,v;Y = ν∈F v 2-p ν ∥c ν ∥ p Y ≤ ν∈F |d ν | p ≤ C(b, p, ξ) p .Hence c ∈ ℓ p v (F; Y) with ∥c∥ p,v;Y ≤ C(b, p, ξ), as required. The second step involves the application of the weighted Stechkin's inequality (see [3, Lem. 3.12] ). This gives that min ∥c -c S ∥ q,v;Y : S ⊂ F, |S| v ≤ k =: σ k (c) q,v;Y ≤ ∥c∥ p,v;Y k 1/q-1/p , for any q ∈ (p, 2] and k > 0. Applying this result with q = 2 implies the existence of the set S 1 (recall that ∥ • ∥ 2,v;Y = ∥ • ∥ 2;Y ) and applying it with q = 1 implies the existence of the set S 2 .We now define the setΛ HCI n =    ν = (ν k ) ∞ k=1 ∈ F : k:ν k ̸ =0 (ν k + 1) ≤ n, ν k = 0, k > n    ⊂ F. (D.26)Notice that Λ HCI n is isomorphic to an index set in N n 0 by the natural restriction map. Lemma D.17 where f S = ν∈S c ν Ψ ν . Moreover, if Y is a Hilbert space, then we also have that∥f -f S ∥ L 2 ϱ (D;Y) ≤ C(b, p, ξ) • k 1/2-1/p + n 1/2-1/p . (D.28)Here C(b, p, ξ) > 0 is a constant depending on b, p and ξ only.Proof of Lemma D.17. The previous lemma implies that there exist index sets S 1 , S 2 ⊂ F with This follows from the same approach and [5, Cor. 8.2] once more.|S 1 | v , |S 2 | v ≤ k/2 such that ∥c -c S1 ∥ 2;Y ≤ C(b, p, ξ) • k 1/2-1/p , ∥c -c S2 ∥ 1,v;Y ≤ C(b, p, ξ) • k 1-1/p . Now define S = S 1 ∪ S 2 ∩ Λ and notice that S ⊆ Λ and |S| v ≤ |S 1 | v + |S 2 | v ≤ k. Hence S ∈ S. Since v ν ≥ u ν = ∥Ψ ν ∥ L ∞ ϱ (D) , we have (using [5, Lem. 5.1]) ∥f -f S ∥ L ∞ ϱ (D;Y) ≤ ∥c -c S ∥ 1,u;Y ≤ ∥c -c S2 ∥ 1,v;Y + ∥c -c Λ ∥ 1,u;Y . Hence,Step 4: Bounding the polynomial error terms. It remains to bound the error terms F -q • ι in (D.37). Using (A.I), (A.II) and Lemma D.17, we now notice that∥F -q • ι∥ L ∞ µ (X ;Y) = ∥f -q∥ L ∞ ς (D;Y) ≲ ∥f -q∥ L ∞ ϱ (D;Y) ≲ C(b, p, ξ) • (k 1-1/p + n 1-1/p). Recall that n ≥ k. Therefore, using this, (D.31) and (3.6), we see thata Y ∥F -q • ι∥ L ∞ µ (X ;Y) / m/L + ∥F -q • ι∥ L 2 µ (X ;Y) ≲ a Y ∥F -q • ι∥ L ∞ µ (X ;Y) ≲ a Y • C(b, p, ξ) • (m/L) 1-1/p = E app,2and likewise for the L ∞ µ -norm bound. Combining this with (D.37) now completes the proof.Proof of Theorem 3.1 when Y is a Hilbert space. The two differences in theorem statement when Y is a Hilbert space are: (i) the L 2 -norm error is with respect to the stronger Bochner norm, and (ii) the approximation error terms E app,q , q = 2, ∞, are smaller by a factor of 1/2 in the exponent (recall (3.6)). We treat both issues separately.For the (i), we commence with the supporting results in §D.2. First, we note that Lemmas D.1 and D.2 also hold in the Bochner L 2 -norm, since these results already give upper bounds involving the Pettis L 2 -norm. Next, we observe that the proof of Lemma D.3 is readily adapted to yield an equivalent result in the Bochner L 2 -norm with the same constant δ. The same therefore applies to Lemma D.4.We next consider the analysis of (D. Up to this point, we have not used the fact that Y is a Hilbert space. We now need this property. As in §D.4, the next step is to establish that (D.38) holds with high probability. Lemma D.7 is unchanged, therefore our focus is on Lemma D.8. We now describe the steps needed to modify the proof of this lemma to assert (D. i=1 . Since b ij = 0, ∀j, whenever c νi = 0, this vector also satisfies ∥b j ∥ 0,v ≤ 2k. Hence, by Lemma D.7 and Parseval's identity twice, ∥h∥ 2 disc,ς ≥ d Y j=1 (θ --(1 + θ + c 5 δ)∥b j ∥ 2 2 = (θ --(1 + θ + c 5 δ) N i=1 ∥c νi ∥ 2 2 = (θ --(1 + θ + c 5 δ)∥h∥where ⌊k⌋+r) . Finally, we define the approximationN = C          N ν1 . . . N ν |S∪Γo | 0 . . . 0          (E.4) and C ∈ R d Y ×(F = D Y • N • E X (E.5)and, for convenience,F = p • E X . E.2 Estimation of the DNN minimizer Lemma E.1 ( N is a minimizer). If B is full rank, then F (X i ) = Y i , ∀i = 1, . . . , m.Therefore, N is a minimizer of (2.5).Proof. Observe thatF (X i ) = p • E X (X i ) = p • E X (X i ) + √ r r i=1 (B) ij (B † e + yz) j = p • E X (X i ) + √ r(B(B † e + yz)) i = p • E X (X i ) + √ r(e) i = Y i .Here, in the penultimate step we use the facts that Byz = yBz = 0 since z ∈ N (B) and BB † = I since r ≥ m and B is full rank by assumption. This gives the first result.For the second result, we recall thatY i is a closest point to Y i from Y = D Y (R d Y ). Therefore, for any DNN N , 1 m m i=1 ∥Y i -D Y • N • E X (X i )∥ 2 Y ≥ 1 m m i=1 ∥Y i -Y i ∥ 2 Y = 1 m m i=1 ∥Y i -D Y • N • E X (X i )∥ 2 Yas required.Lemma E.2 (Bounding F in terms of F ). Suppose that (D.7) holds with α ≥ c 0 and α -δ √ k ≥ c 1 , and also that√ r∥ι d X -ι d X • D X • E X ∥ L 2 µ (X ;ℓ ∞ (N)) ≤ c 2 , where c 0 , c 1 , c 2 > 0 are suitable universal constants. Then the approximation F satisfies |||F -F ||| L 2 µ (X ;Y) ≲ |||F -F ||| L 2 µ (X ;Y) + (1 + δ √ r)δ √ k 1 + √ m √ rσ min (B) ∥F ∥ L ∞ µ (X ;Y) + 1 √ m ∥E∥ 2;Y + √ m(1 + δ √ r) √ rσ min (B)   1 m m i=1 ∥ Y i -F (X i )∥ 2 Y + 1 √ m ∥E∥ 2;Y   + (1 + δ √ r)∥z∥ 2 and ∥F -F ∥ L ∞ µ (X ;Y) ≲ ∥F -F ∥ L ∞ µ (X ;Y) + (1 + δ) √ rδ √ k 1 + √ m √ rσ min (B) ∥F ∥ L ∞ µ (X ;Y) + 1 √ m ∥E∥ 2;Y + √ m(1 + δ) σ min (B)   1 m m i=1 ∥ Y i -F (X i )∥ 2 Y + 1 √ m ∥E∥ 2;Y   + (1 + δ) √ r∥z∥ 2 .Proof. By the triangle inequality,|||F -F ||| L 2 µ (X ;Y) ≤ |||F -F ||| L 2 µ (X ;Y) + ||| F -F ||| L 2 µ (X ;Y) ∥F -F ∥ L ∞ µ (X ;Y) ≤ ∥F -F ∥ L ∞ µ (X ;Y) + ∥ F -F ∥ L ∞ µ (X ;Y) . (E.6)Consider the second term. We have||| F -F ||| L 2 µ (X ;Y) = |||p -D Y • N ||| L 2 ς (D;Y) ≤ |||p -p||| L 2 ς (D;Y) + |||q||| L 2 ς (D;Y) ∥ F -F ∥ L ∞ µ (X ;Y) = ∥p -D Y • N ∥ L ∞ ς (D;Y) ≤ ∥p -p∥ L ∞ ς (D;Y) + ∥q∥ L ∞ ς (D;Y) , (E.7)where q = p -p = r i=1 (B † e + yz) i N ei . Lemma D.10 and (A.III) give that|||p -p||| L 2 ς (D;Y) = ∥p -p∥ L ∞ ς (D;Y) ≤ δ √ k|||p||| L 2 ϱ (D;Y) . (E.8)Now consider the other term in (E.7). Define q = r i=1 (B † e + yz) i Ψ ei . Then Lemma D.10 and (A.III) once more give that|||q||| L 2 ς (D;Y) ≤ |||q||| L 2 ς (D;Y) + |||q -q||| L 2 ς (D;Y) ≤ |||q||| L 2 ς (D;Y) + δ √ r|||q||| L 2 ϱ (D;Y) ∥q∥ L ∞ ς (D;Y) ≤ ∥q∥ L ∞ ς (D;Y) + ∥q -q∥ L ∞ ς (D;Y) ≤ ∥q∥ L ∞ ς (D;Y) + δ √ r|||q||| L 2 ϱ (D;Y).Here, we also used the fact that |Γ o | = r. We now apply Lemmas D.1 and D.4 and (A.III) once more to get|||q||| L 2 ς (D;Y) ≲ (1 + δ √ r)|||q||| L 2 ϱ (D;Y) , ∥q∥ L ∞ ς (D;Y) ≲ (1 + δ) √ r|||q||| L 2 ϱ (D;Y) . (E.9)We next analyze the term |||q||| We now use the definition of e and the inequality |||e||| 2;Y ≤ ∥e∥ 2;Y to obtain|||q||| L 2 ϱ (D;Y) ≤ √ m √ rσ min (B) 1 m m i=1 ∥ Y i -p • E X (X i )∥ 2 Y + ∥z∥ 2 .We next apply the triangle inequality and Lemma D.10 once more to get|||q||| L 2 ϱ (D;Y) ≤ √ m √ rσ min (B)   1 m m i=1 ∥ Y i -p • E X (X i )∥ 2 Y + ∥p -p∥ disc,ς   + ∥z∥ 2 ≤ √ m √ rσ min (B)   1 m m i=1 ∥ Y i -p • E X (X i )∥ 2 Y + δ √ k|||p||| L 2 ϱ (D;Y)   + ∥z∥ 2 .Combining this with (E.7) and (E.9), we deduce that||| F -F ||| L 2 µ (X ;Y) ≲ (1 + δ √ r) δ √ k 1 + √ m √ rσ min (B) |||p||| L 2 ϱ (D;Y) + ∥z∥ 2 + √ m √ rσ min (B)   1 m m i=1 ∥ Y i -p • E X (X i )∥ 2 Y + 1 √ m ∥E∥ 2;Y   and ∥ F -F ∥ L ∞ µ (X ;Y) ≲ (1 + δ) √ r δ √ k 1 + √ m √ rσ min (B) |||p||| L 2 ϱ (D;Y) + ∥z∥ 2 + √ m √ rσ min (B)   1 m m i=1 ∥ Y i -p • E X (X i )∥ 2 Y + 1 √ m ∥E∥ 2;Y   .It remains to bound the term |||p||| L 2 ϱ (D;Y) . Using (D.7) and the fact that p is a minimizer of (D.3) and that the zero polynomial is feasible for (D.3), we get EQUATIONP σ min ≤ 1 -c 1 n/m -t ≤ 2 exp(-c 2 mt 2 ), where c 1 , c 2 > 0 are universal constants. Lemma E.4 (Bounding σ min (B)). Suppose that √ mδ ≤ √ ω/8, 3 (5+ξ)/2 √ r 2 ∥ι d X -ι d X • D X • E X ∥ L ∞ µ (X ,ℓ ∞ (N)) ≤ √ ω 8 , (E.10)where ω is the variance of the univariate probability measure as specified in Theorem 3.2 andd X ≥ r ≥ c (m + log(2/ϵ)) (E.11)for some universal constant c ≥ 1. Then, with probability at least 1 -ϵ, the matrix B is full rank andσ min (B) ≥ √ ω/4.Proof. Define the matricesB ′ = 1 √ r Ψ ej • E X (X i ) m,r i,j=1 , B ′′ = 1 √ r Ψ ej • ι d X (X i ) m,r i,j=1 . Then, since r ≥ m, σ min (B) = inf ∥B ⊤ d∥ 2 : d ∈ C m , ∥d∥ 2 = 1 ≥ σ min (B ′ ) -∥(B -B ′ ) ⊤ ∥ 2 = σ min (B ′ ) -∥B -B ′ ∥ 2 ≥ σ min (B ′′ ) -∥B -B ′ ∥ 2 -∥B ′ -B ′′ ∥ 2 . Now, for any c ∈ C r , ∥(B -B ′ )c∥ 2 2 = 1 r m i=1   r j=1 N ej (x i ) -Ψ ej (x i ) c j   2 ≤ 1 r m i=1 δ 2 ∥c∥ 2 1 ≤ mδ 2 ∥c∥ 2 2 .We deduce that ∥B -B ′ ∥ 2 ≤ √ mδ ≤ √ ω/8. Hence σ min (B) ≥ σ min (B ′′ ) - √ ω/8 -∥B ′ -B ′′ ∥ 2 .Now let c ∈ C r and p = r i=1 c i Ψ ei be the corresponding polynomial. Then, by (A.I), (A.III) and Lemma D.2,∥(B ′ -B ′′ )c∥ 2 = 1 r m i=1 p • ι d X (X i ) -p • ι d X • D X • E X (X i ) 2 ≤ Lip(p, B ∞ (N), R)∥ι d X -ι d X • D X • E X ∥ L ∞ µ (X ,ℓ ∞ (N)) ≤ 3 (5+ξ)/2 √ r 2 ∥p∥ L 2 ϱ (D) ∥ι d X -ι d X • D X • E X ∥ L ∞ µ (X ,ℓ ∞ (N)) ≤ √ ω 8 ∥c∥ 2 .Here, in the third step we used the fact that|Γ o | v = 3 5+ξ r, since u ei = √ 3. We deduce that ∥B ′ -B ′′ ∥ 2 ≤ √ ω/8 and therefore σ min (B) ≥ σ min (B ′′ ) - √ ω/4.It remains to show that σ min (B ′′ ) ≥ √ ω/2 with high probability. By construction, Ψ ej (x) = √ 3x j . Now recall that the pushforward ς is a tensor-product of a univariate probability measure supported in [-1, 1] with mean zero and variance ω > 0. Therefore (B ′′ ) ⊤ = √ 3 √ ω √ r A, where A = (ι(X j ) i / √ ω) r,m i,j=1 ∈ R r×m .P(σ min (B ′′ ) ≤ √ ω/2) ≤ P(σ min (A/ √ r) ≤ 1/(2 √ 3)) ≤ ϵ.This gives the result.We are now ready to complete the proof of Theorem 3.2.Proof of Theorem 3.2, Statement (A). We divide the proof into a series of steps.Step 1: Setup. Let m, δ, ϵ and L be as in the theorem statement. We once more assume without loss of generality that δ ≤ 1/5. Let n be as in (D.29) and Λ = Λ HCI n , let ξ be as in (D.30) andk = m c 1 L , (E.12)where c 1 ≥ 1 is a constant that will be chosen in the next step. Letδ = min 2 -m /r 2 , √ ω/(8 √ m) ,where ω is the variance of the univariate probability measure andr = ⌈c 2 (m + log(1/ϵ))⌉, (E.13)where c 2 ≥ 2 will also be chosen in the next step. Note that d X ≥ r > m by assumption. Finally, let S be as in (D.1), Γ be as in (E.1) and {N ν } ν∈Γ be the corresponding family of tanh DNNs ensured by Lemma D.9. Finally, let F be given by (E.5), with DNN N as in (E.4).Step 2: N is a minimizer. By construction, widthEQUATIONUsing the values of k and r, we see that width( N ) ≲ (m + log(1/ϵ))(m/L) δ , depth( N ) ≲ log(m/L).Therefore, N ∈ N is feasible for (2.5). Lemma E.1 now implies that N is a minimizer, provided B is full rank. We will show that this holds in the next step.Step c 0 • k • (log(eN ) • log 2 (k) + log(6/ϵ)) ≲ k • log 2 (m) log 2 (m) + log(1/ϵ) + log(1/ϵ) ≲ k • log 4 (m) + log(1/ϵ) . = k • L(m, ϵ).Hence, due to the definition of k, we get that (D.13) holds once more with probability at least 1 -ϵ/4. The remainder of this step is identical to Step 2 of the proof of Theorem 3.1.Step 5: Error analysis. We now apply Lemma E.2 to the approximation F . Since σ min (B) ≳ √ ω ≳ 1, δ ≤ δ √ r ≤ 1 and r ≥ m, we deduce that|||F -F ||| L 2 µ (X ;Y) ≲ |||F -F ||| L 2 µ (X ;Y) + δ √ k ∥F ∥ L ∞ µ (X ;Y) + 1 √ m ∥E∥ 2;Y + ∥z∥ 2 + 1 m m i=1 ∥ Y i -F (X i )∥ 2 Y + 1 √ m ∥E∥ 2;Y and ∥F -F ∥ L ∞ µ (X ;Y) ≲ ∥F -F ∥ L ∞ µ (X ;Y) + √ rδ √ k ∥F ∥ L ∞ µ (X ;Y) + 1 √ m ∥E∥ 2;Y + √ r∥z∥ 2 + √ m   1 m m i=1 ∥ Y i -F (X i )∥ 2 Y + 1 √ m ∥E∥ 2;Y  with probability at least 1 -ϵ/2, where F = p • E X and p is the corresponding minimizer of (D.3). We now appeal to Theorem D.6 with σ = 1 and τ = 0, recalling that α ≳ 1 and δ√ k ≲ 1, to get that |||F -F ||| L 2 µ (X ;Y) ≲ |||F -Q • F ||| L 2 µ (X ;Y) + |||F -Q • F ||| disc,µ + π Q |||F -q • E X ||| L 2 µ (X ;Y) + |||F -q • E X ||| disc,µ + 1 √ m ∥E∥ 2;Y + δ √ k∥F ∥ L ∞ µ (X ;Y) + ∥z∥ 2 + 1 m m i=1 ∥ Y i -F (X i )∥ 2 Yand, since δ √ r ≲ 1, ∥F -F ∥ L ∞ µ (X ;Y) ≲ ∥F -Q • F ∥ L ∞ µ (X ;Y) + √ k|||F -Q • F ||| disc,µ + π Q ∥F -q • E X ∥ L ∞ µ (X ;Y) + √ k|||F -q • E X ||| disc,µ + ∥E∥ 2;Y + √ rδ √ k∥F ∥ L ∞ µ (X ;Y) + √ r∥z∥ 2 + √ m 1 m m i=1 ∥ Y i -F (X i )∥∥ Y i -Y i ∥ Y ≤ ∥F (X i ) -Q • F (X i )∥ Y + ∥E i ∥ Y .We now use this, the fact that p is a minimizer and Q • q is feasible in combination with triangle inequality to get1 m m i=1 ∥ Y i -F (X i )∥ 2 Y ≤ 1 m m i=1 ∥Y i -Q • q • E X (X i )∥ 2 Y + |||F -Q • F ||| disc,µ + 1 √ m ∥E∥ 2;Y ≤ |||F -Q • q • E X ||| disc,µ + |||F -Q • F ||| disc,µ + 2 √ m ∥E∥ 2;Y ≤ π Q |||F -q • E X ||| disc,µ + 2|||F -Q • F ||| disc,µ + 2 √ m ∥E∥ 2;Y .Now let f ∈ H(b) be the function asserted by (A.II) and q = f S be the polynomial asserted by Lemma D.17. Substituting this into the previous expression, recalling (D.35) and using the definition of δ, we deduce that|||F -F ||| L 2 µ (X ;Y) ≲ |||F -Q • F ||| L 2 µ (X ;Y) + |||F -Q • F ||| disc,µ + π Q E 2 (F, q) + ∥z∥ 2 + 2 -m + E samp,2 ∥F -F ∥ L ∞ µ (X ;Y) ≲ ∥F -Q • F ∥ L ∞ µ (X ;Y) + √ m|||F -Q • F ||| disc,µ + π Q E ∞ (F, q) + √ r∥z∥ 2 + 2 -m + E ′ samp,∞with probability at least 1 -ϵ/2. Here E 2 (F, q) is as in (D. 24 ), E ∞ (F, q) is as in (D. 24|||F -F ||| L 2 µ (X ;Y) ≲ a Y ∥F -q • ι∥ L ∞ µ (X ;Y) / m/L + ∥F -q • ι∥ L 2 µ (X ;Y) + E X ,2 + E Y,2 + ∥z∥ 2 + 2 -m + E samp,2 .(E.14)The second bound above is identical to the corresponding bound in (D.36), except with E opt,∞ replaced by √ r∥z∥ 2 +2 -m , E samp,∞ and E ∞ (F, q) replaced by E ′ samp,∞ and E ∞ (F, q), respectively, and withk replaced by √ m. We once more follow the same arguments as in Step 3, with these changes. This yields the corresponding version of (D.37), which is∥F -F ∥ L ∞ µ (X ;Y) ≲ a Y √ L∥F -q • ι∥ L ∞ µ (X ;Y) + √ m∥F -q • ι∥ L 2 µ (X ;Y) + E ′ X ,∞ + E ′ Y,∞ + √ r∥z∥ 2 + 2 -m + E ′ samp,∞ .(E.15)Here E ′ X ,∞ = LE X ,∞ and E ′ Y,∞ = √ LE Y,∞ .Lemma D.9 as the composition of an affine map defined by the fundamental theorem of algebra and a tanh DNN that approximately multiplies m(Γ) numbers. In this specific case, we haveΨ ei (x) = x i .Therefore, the corresponding affine map (D.18) A ei : R d X → R m(Γ) has a bias vector that is all ones, except for a single entry that has value zero. We deduce that the bias vector b of the full DNN (E.4) contains a subblock of size rm(Γ) that is all ones, except for r zeroes. There are As noted, the proof of this theorem is based on [7, Thm. 4.4] . We recap the details as they will also be needed in the proof of the next result. First, we recall some basic definitions. See [82] or [29, Ch. 10] for more details. Let K be a subset of a normed space (X , ∥ • ∥ X ). Then its Gelfand m-width is Proof of Theorem G.1. We proceed in a similar series of steps to those of the last proof.Step 1: Setup. and notice that we may writeA = [A 1 • • • A r ] ,where A l = (ι(X i )) π((l-1)(m+1)+j) m,m+1 i,j=1 . Let σ be the one-dimensional probability measure associated with ς. Then notice that each A l is a random matrix whose entries are drawn i.i.d. from σ. Since σ is supported in [-1, 1], we deduce that the A l are independent subgaussian random matrices with the same distribution. Write γ for this distribution. Let t 1 , . . . , t r > 0 and E l,t l be the event E l,t l = ∃u ∈ N (A l ) : ∥u∥ 2 = 1, ∥u∥ ∞ < t l /(m + 1) , l = 1, . . . , r.(G.1)We will make a suitable choice of t 1 , . . . , t r later. For a > 0, the function x → max{x, a -x} is minimized at x = a/2 and takes value a/2 there. We deduce that sup (G.3)Step 3: Bounding the expected error. Fix l = 1, . . . , r and suppose the event E l,t l defined (G.1) occurs. Let u l be the corresponding vector and define c l = b π(l(m+1)) ∥u l ∥ ∞ u l , l = 1, . . . , r. 


Decision: The paper presents some results on learning holomorphic operators between Banach spaces using deep neural networks. Theoretical results about generalization bounds and optimality of these bounds are given along with some numerical experiments relating to several different parametric PDEs. Based on the reviews, it is agreed that the strengths of the paper are its powerful, well-formulated theoretical results, its rigorous presentation, and its originality. All reviews furthermore indicated no major weaknesses in the paper, indicating that the results are technically sound and interesting to the community.

Regarding the rebuttal, as all reviews were quite positive, there was little to discuss besides superficial changes (figure sizes, minor clarifications/rewordings). The authors responded to all questions and criticisms of the reviewers sufficiently and the reviewers participated in the discussion as necessary.

Since all reviewers were in agreement with respect to the quality of the submission, I did not see it necessary to have a post-rebuttal discussion.

I recommend accepting the paper.

Reviews: review: paper_summary: paper_summary: The paper addresses the challenge of learning operators between Banach spaces, focusing on holomorphic operators which have many applications. The authors employ standard (DNN) architectures with constant width greater than depth and $\ell_2$-loss minimization. They identify a family of DNNs that achieve optimal generalization bounds for such operators. For fully-connected architectures, they show there are infinitely many solutions that offer optimal performance. The results demonstrate that deep learning achieves the best possible generalization bounds, and no method can surpass these bounds, up to logarithmic terms.
strength_weakness: strength_weakness: I think some of the figures of the experiments results are too small, making it a bit difficult to recognize the results. It would be good if the authors could review this issue.
questions: questions: No questions.
limitations: limitations: Theorem 3.2 only asserts that some minimizers are ‘good’, not all.
paper_summary: paper_summary: Authors provide theoretical results on sample efficiency for learning holomorphic neural operators $\mathcal{X}\rightarrow\mathcal{Y}$ (results are available for both Hilbert and Banach spaces) with a neural network of a particular structure. The structure of the network is encoder $\mathcal{X}\rightarrow\mathbb{R}^{d_1}$, decoder $\mathcal{Y}\rightarrow\mathbb{R}^{d_2}$ and DNN $\mathbb{R}^{d_1}\rightarrow \mathbb{R}^{d_2}$, with both encoder and decoder considered to be given and DNN trained on $l_2$ loss in a usual supervised manner.

Besides the upper bound on approximation error, authors also provide lower bound and conclude that active learning (deliberate collection of data used to train neural operator) is not justified, since non-adaptive training is already optimal up to log terms.

Authors illustrate their theoretical results (relative test error vs number of samples for selected DNN configurations) with several numerical experiments, including diffusion equation in $D=2$, Navier-Stokes-Brinkman equation in $D=2$ and Boussinesq equation in $D=3$.
strength_weakness: strength_weakness: Overall, in my opinion, the research is of excellent quality, so I can not point to any major weaknesses.

The two main issues that can be relevant are: (i) the requirements for the encoder and decoder are not sufficiently discussed, and I find restricting the assumption that they are already available, and (ii) the significance of the numerical confirmation is not apparent.
questions: questions: 1. After line 102 there is a definition of Bochner and Pettis $L^p$ norms. Presumably, the later expression is a definition of $|||F|||$ which is not apparent.
2. I want to clarify several statements on the encoder and decoder. 
   1. In (1.3) authors introduce encoder $\mathcal{E}$ and decoder $\mathcal{D}$, and later approximate encoder and decoder $\widetilde{\mathcal{E}}$, $\widetilde{\mathcal{D}}$ after line 114. What is the difference between those?
   2. In particular, in the discussion (lines 145, 146) of assumption A.IV (line 131) that requires both $\mathcal{E}$ and $\mathcal{D}$ to be linear and bounded, authors claim that $\widetilde{\mathcal{E}}$ is not required to be linear. It is not clear to me how this agrees with the fact that $\mathcal{E}_{\mathcal{X}}$ is defined in A.III as a composition of approximate encoder and decoder. Can the authors please clarify that?
   3. According to (1.3) encoder and decoder are defined for the approximate operator $\hat{F}$. This means, as I understand, they are not rigidly fixed by $F$ since presumably $F$ can be approximated with many different encoders and decoders. Approximate encoders and decoders are not formally defined, so it is not clear what they approximate. Does it all mean that authors consider not a continuous operator learning problem but suppose that discretization (of the space of parameters and $\mathcal{Y}$) is already performed ($\mathcal{E}$ and $\mathcal{D}$ are selected) and $\widetilde{\mathcal{E}}$ and $\widetilde{\mathcal{D}}$ approximate $\mathcal{E}$ and $\mathcal{D}$?
   4. Why there are $\widetilde{\mathcal{E}}\_{\mathcal{X}}$ and $\widetilde{\mathcal{D}}\_{\mathcal{X}}$ but not the approximate encoder and decoder for $\mathcal{Y}$?
3. Several questions on Theorem 3.1:
   1. Why $a_{\mathcal{Y}}$ is called approximation error? What is the meaning of this term (it looks like some kind of volume)? In the ideal setting $\mathcal{D}\_{\mathcal{Y}}\circ\mathcal{E}\_{\mathcal{Y}} = I$, so we have a perfect approximation that has nonzero approximation error.
   2. Since $\lim_{m\rightarrow\infty}m \big/ L(m) = \infty$ many terms in the upper bound will diverge unless the encoder and decoder are of sufficient quality. More specifically, the encoder and decoder should improve their accuracy as $\sim m^{-1/2}$ for u.b. to be useful in case of $2$ norm and as $\sim m^{-1}$ in case of Chebyshev norm. Can the authors please discuss this consequence of their result in more detail?
   3. Can the authors also explain the consequence of the optimization error term? It seems that $\tau$ that appears in it is to some extent beyond control and completely independent of $m$.
4. Questions about experiments:
   1. In experiments authors strive to confirm theoretical scaling $m^{-1}$. Since the encoder and decoder are perfect, and noise is absent, the only two sources of error are approximation error and optimization error. In my opinion, one can clearly see the signs of saturation for larger $m$. This is also clear from operator learning literature that relative error rarely fell below $10^{-3}$. Can the authors please comment on that? Is it because optimization error prevails for large $m$?
   2. Dimensions $d=4,8$ in Figures 1, 2, 3 and in the text (e.g., lines 296, 297: "PDE in $d=4$ and $d=8$ dimensions") are slightly misleading because they refer to the dimension of the parameter space and can be confused with the number of spatial dimensions.
   3. All results are given for perfect encoder and decoder. Can the authors share some thoughts on whether it is possible to test their upper bound with an imperfect decoder and encoder? What are those experiments might be? Is it possible to create a controllable setting to test terms in ub that contain encoder and decoder errors?
limitations: limitations: yes
paper_summary: paper_summary: This work examines learning holomorphic operators between Banach spaces. It shows that these learning tasks can be tackled using encoder-decoder approaches with deep neural networks (DNNs), achieving optimal approximation rates up to log factors. The architectures used are problem-agnostic and depend only on the amount of available training data. The study also delves into theoretical aspects of the associated training problem for standard feed-forward networks and presents numerical findings to support its claims.
strength_weakness: strength_weakness: Overall I do not have any serious concerns about this work and consider it to be a valuable contribution. Some questions are listed below.
questions: questions: * Perhaps I missed it, but I could not locate a precise reference discussing the holomorphicity of the Navier-Stokes-Brinkman and Boussinesq equations, similar to Section B.3.2. for the Diffusion problem. It would be helpful if such a discussion could be included in the manuscript.

* Could you briefly summarize to what extent Theorems 4.1 and 4.2 require different techniques compared to those used in [7]? How does the use of a general Banach space $\mathcal{X}$ complicate the analysis?

* Regarding the experimental results, do you have any insights into why different activation functions perform differently across individual tasks? Can this behavior be explained using the results from Theorems 3.1 and 3.2?

* Do the authors think the results from Theorem 3.2 could be extended to include information on optimization procedures used for solving (2.5) (e.g., stochastic gradient descent and variants thereof)? Such an extension would enhance the practical relevance of the result. I am only seeking an educated guess or a high-level opinion, without expecting detailed analysis.

* I would appreciate it if the authors could briefly summarize the differences between the Banach space and Hilbert space cases for $\mathcal{Y}$. Specifically, could you indicate precisely where the inner product structure is utilized in the proofs to achieve better rates?  Also please double check formulas (3.3) and (3.5) in this context, currently it seems that they are the same.
limitations: limitations: Limitations and social impact have been discussed.

Comments: comment: Thanks for your response!
Thanks for your second positive assessment and further comments.

Your comment about tanh networks is a good one. We agree this is a good starting point for further investigations and will take it on board in our future work.

Thanks also for flagging the point about the Pettis norm. As another also referee pointed out, this was a typo in the equation between lines 102 and 103. We will add this notation to that equation in the revision.
I am pleased with the explanations provided by the authors and am willing to increase my score to 8. I am also looking forward to follow-up work that may address my first and fourth points (as mentioned, I did not expect detailed responses here). Here are a few additional comments:

* Regarding the tanh network: One possible reason for its behavior could be the vanishing gradient property of the tanh function, which may lead to suboptimal performance of the optimization algorithm in larger (deeper) networks. However, this is merely a rough guess, and further investigation is needed since there could be many other contributing factors, as the authors have noted.
* Regarding the Pettis norm: It might be helpful for the reader to introduce the associated notation ($|||.|||$) in the equation between 102 and 103. I believe this would address the source of my confusion.