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.

§ 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 learn

\begin{equation}
\label{eq:operator_equation}
F: \cX \to \cY, \qquad X \in \cX \mapsto F(X) \in \cY,
\end{equation}

where $\cY$ 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 $\mu$ be a probability measure on $\cX$. Given noisy training data

{ (X_i,F(X_i) + E_i) }^m_i=1

where $X_1,\ldots,X_m \sim_{\mathrm{i.i.d.}} \mu$ and $E_i$ is noise, a typical operator learning methodology consists of three objects: an approximate encoder $\cE_{\cX} : \cX \rightarrow \bbR^{d_{\cX}}$, an approximate decoder $\cD_{\cY} : \bbR^{d_{\cY}} \rightarrow \cY$ and a DNN $\widehat{N} : \bbR^{d_{\cX}} \rightarrow \bbR^{d_{\cY}}$. It then approximates $F$ as

F ≈F : = _ ∘N ∘_.

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 $\widehat{F}$ from as little training data $m$ as possible.

§.§ Contributions

As observed in, e.g., [Chen et al., 2023, Liu et al., 2024], 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 <ref>), with arbitrary approximate encoders $\cE_{\cX}$ and $\cD_{\cY}$. As we explain in <ref> (see also <cit.>, <cit.> and [Herrmann et al., 2022]) these operators are relevant in many applications of operator learning, notably those involving parametric PDEs. The main contributions of this work are as follows.

* We consider operators taking values in general Banach spaces. As noted, the vast majority of existing work (with the notable exception of [Chen et al., 2023]) considers Hilbert spaces.

* We consider standard feedforward DNN architectures (constant width, width exceeds depth) and training procedures ($\ell^2$-loss minimization).

* (Theorem <ref>) 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 training-data-op.

* These DNN architectures are problem agnostic; they just depend on $m$ only. In particular, the architectures are completely independent on the regularity assumptions of target operator.

* (Theorem <ref>) We then show that training problems based on any family of fully-connected DNNs possess uncountably many minimizers that achieves the same generalization bounds.

* (Theorems <ref>-<ref>) We provide bounds in both the $L^2_{\mu}$- and $L^{\infty}_{\mu}$-norms that hold in high probability, rather than just expectation.

* (Theorems <ref>-<ref>) 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.

* Finally, we present a series of experiments demonstrating the efficacy of DL on challenging problem such as the parametric diffusion, Navier-Stokes-Brinkman and Boussinesq PDEs.

§.§ Relation to previous work

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 [Cao et al., 2023, Kovachki et al., 2021, Li et al., 2020, Li et al., 2020, Kovachki et al., 2023, Lanthaler, 2023, Lu et al., 2021, Lanthaler et al., 2023, Li et al., 2023, Hao et al., 2023, Raonić et al., 2023, Tripura and Chakraborty, 2023, Zhang et al., 2023].
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 [de Hoop et al., 2022, Lin et al., 2021, Goswami et al., 2022, Howard et al., 2022, Goswami et al., 2022, Cai et al., 2021, Osorio et al., 2022, Lin et al., 2021, Oommen et al., 2022, Zhu et al., 2023, Lu et al., 2022, Karniadakis et al., 2021, Yin et al., 2022, Zou et al., 2022, Yin et al., 2022, Kontolati et al., 2023, Goswami et al., 2023, Li et al., 2020, Li et al., 2022, Wu et al., 2023, Michałowska et al., 2023, Yeung et al., 2019, Koric et al., 2023, Zappala et al., 2023, Renn et al., 2023, Peng et al., 2023].

On the theoretical side, universal approximation theorems for operator learning have been developed in [Lu et al., 2021, Lu et al., 2019, Lanthaler et al., 2022, Kovachki et al., 2021] 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 [Bhattacharya et al., 2021, Liu et al., 2024, Lanthaler et al., 2022, Schwab et al., 2023, Chen et al., 2023] and references therein. However, learning Lipschitz operators suffers from a curse of parametric complexity [Lanthaler and Stuart, 2023, Lanthaler and Stuart, 2024], 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 [Lanthaler et al., 2022, Franco et al., 2023, Marcati and Schwab, 2023, Herrmann et al., 2022, Deng et al., 2022] and elsewhere.

However, these works do not consider the generalization error, i.e., the error incurred when learning the approximation F-approx from the finite training data training-data-op.
This is particularly important in applications of operator learning where data is obtained through expensive numerical PDE solves. Such problems are highly data-starved. Several works have tackled this question from the perspective of statistical learning theory and nonparametric estimation [Lanthaler et al., 2022, Chen et al., 2023, Liu et al., 2024], but only for Lipschitz operators. As observed in <cit.>, this approach generally leads to a best $\ord{m^{-1/2}}$ decay of the $L^2_{\mu}$-norm error with respect to $m$. Theorem <ref> shows that such a rate is strictly suboptimal. Our generalization bounds in Theorems <ref>-<ref> do not use such techniques, and yield near-optimal rates in both the $L^2_{\mu}$- and $L^{\infty}_{\mu}$-norms. See also [Franco and Brugiapaglia, 2024] for some related work in this direction for reduced-order modelling with convolutional autoencoders.

Our work is inspired by recent work on learning holomorphic, Banach-valued functions [Adcock et al., 2021, Adcock et al., 2022, Adcock et al., 2024], as opposed to operators. We extend both these works, in particular, [Adcock et al., 2022], to learning holomorphic operators. We also significantly improve the error decay rates in [Adcock et al., 2022] with respect to $m$ and show they can be achieved using substantially smaller DNNs with standard training (i.e., $\ell^2$-loss minimization). See Remarks <ref>-<ref>. Our theoretical guarantees fall into the category of encoder-decoder-nets [Kovachki et al., 2024], which includes the well-known PCA-Net [Bhattacharya et al., 2021] and DeepONet [Lu et al., 2019] frameworks. As in other recent works [Lanthaler et al., 2022, Liu et al., 2024, Chen et al., 2023, Franco and Brugiapaglia, 2024], in Theorems <ref>-<ref> we assume the encoder/decoder pair $(\cE_{\cX},\cD_{\cY})$ in F-approx have been learned, and focus on the generalization error when training the DNN $\widehat{N}$.

§ NOTATION, ASSUMPTIONS, SETUP AND EXAMPLES

§.§ Notation

Let $(\cX,\nm{\cdot}_{\cX})$ and $(\cY,\nm{\cdot}_{\cY})$ be Banach spaces and $\mu$ be a probability measure on $\cX$. Let $(\cY^*,\nm{\cdot}_{\cY^*})$ be the dual of $\cY$ and $B(\cY^*)$ be its unit ball. The Bochner and Pettis $L^p$-norms of a (strongly and weakly, respectively) measurable operator $F : \cX \rightarrow \cY$ are defined as

F_L^p_μ(;) = ( ∫_ F(X)^p_ μ(X) )^1/p  and  sup_y^* ∈B(^*) ( ∫_ | y^*(F) |^p μ(X) )^1/p,

respectively, for $1 \leq p < \infty$, and analogously for $p = \infty$ (see, e.g., [Aliprantis and Border, 2006, Hytönen et al., 2016]).
Notice that $\tnm{F}_{L^p_{\mu}(\cX ; \cY)} \leq \nm{F}_{L^p_{\mu}(\cX;\cY)}$ for $1 \leq p < \infty$, while $\tnm{F}_{L^{\infty}_{\mu}(\cX ; \cY)} = \nm{F}_{L^{\infty}_{\mu}(\cX ; \cY)}$.

Throughout this work, $\ell^p(\bbN)$, $0 < p \leq \infty$ denotes the standard $\ell^p$ space with (quasi-)norm $\nm{\cdot}_p$. We also define the monotone $\ell^p$ space $\ell^p_{\mathsf{M}}(\bbN)$ as the space of all sequences $\bm{z} = (z_i)^{\infty}_{i=1} \in \bbR^{\bbN}$ whose minimal monotone majorant $\tilde{\bm{z}}  \in \ell^p(\bbN)$. Here $\bm{z} = (\tilde{z}_i)^{\infty}_{i=1}$ is defined as $\tilde{z}_i = \sup_{j \geq i} | z_j|$.

Given a (componentwise) activation function $\sigma$, we consider feedforward DNNs of the form

N : ^n →^k, z ↦N(z) = _L+1 ( σ( _L ( σ( ⋯σ( _0 (z) ) ⋯) ) ) ),

where $\cA_l : \bbR^{N_{l}} \rightarrow \bbR^{N_{l+1}}$ are affine maps, and $N_0 = n$ and $N_{L+2} = k$. We define $\mathrm{width}(N) = \max \{ N_1,\ldots, N_{L+1} \}$ and $\mathrm{depth}(N) = L$.
We denote a class of DNNs of the form (<ref>) with a fixed architecture (i.e., fixed activation function, depth and widths) as $\cN$, and write $\mathrm{width}(\cN) = \max \{ N_1,\ldots, N_{L+1} \}$ and $\mathrm{depth}(\cN) = L$.

§.§ Assumptions and setup

Let $F : \cX \rightarrow \cY$ be the unknown operator we seek to learn and

_ : →^d_, _ : ^d_ →,    _ : →^d_, _ : ^d_ →
be approximate encoders and decoders for $\cX$ and $\cY$, respectively. As mentioned, we assume that these maps have already been learned, and focus on the training of the DNN $\widehat{N}$ in F-approx.

In order to define the precise notion holomorphy for the operator $F$, we require the following. Let $D = [-1,1]^{\bbN}$ and $\varrho$ be the uniform probability measure on $D$. Given $\rho > 1$, we define the Bernstein ellipse $\cE(\rho) = \left\{ (x+x^{-1})/2 : x \in \bbC,\ 1 \leq | x | \leq \rho \right \} \subset \bbC$ and let $\cE(1) = [-1,1]$. Given $\bm{\rho} = (\rho_i)_{i \in \bbN} \geq \bm{1}$, we define the Bernstein polyellipse as $\cE(\bm{\rho}) = \cE(\rho_1) \times \cE(\rho_2) \times \cdots \subset \bbC^{\bbN}$.
[Holomorphic map]
Let $\varepsilon > 0$, $\bm{b} \in \ell^1(\bbN)$ with $\bm{b} \geq \bm{0}$. A Banach-valued function $f : D \rightarrow \cY$ is $(\bm{b},\varepsilon)$-holomorphic if it is holomorphic in the region

\begin{equation}\label{def_R_b_e}
\cR({\bm{b},\varepsilon}) =  \bigcup \Big\lbrace \cE(\bm{\rho}): \bm{\rho} \geq \bm{1},\  \sum_{j=1}^{\infty} \left ( (\rho_j+\rho_j^{-1})/2 -1 \right) b_j \leq  \varepsilon \Big\rbrace \subset \bbC^{\bbN},\quad \bm{b} =(b_j)_{j \in \bbN}.
\end{equation}

See, e.g., [Chkifa et al., 2015, Schwab and Zech, 2019]. As noted in [Adcock et al., 2024] we can, by rescaling $\bm{b}$, assume that $\varepsilon = 1$. For convenience, we define the following unit ball, consisting of all such functions of norm at most one over $\cR(\bm{b},1)$:

(b) = { f : D → $(\bm{b},1)$-holomorphic : f(x)_ ≤1, ∀x ∈(b,1) }.

Let $D = [-1,1]^{\bbN}$ and $\varrho$ be the uniform probability measure.

(A.I) There is a measurable mapping $\iota : \cX \rightarrow \bbR^{\bbN}$ such that pushforward $\varsigma : = \iota \sharp \mu$ is a quasi-uniform measure supported on $D$ and $\iota |_{\mathrm{supp}(\mu)} : \cX \rightarrow \ell^{\infty}(\bbN)$ is Lipschitz with constant $L_{\iota} \geq 0$.

(A.II) The function $F$ has the form $F = f \circ \iota$, where $f \in \cH(\bm{b})$ for some $\bm{b} \in \ell^p_{\mathsf{M}}(\bbN)$ and $0 < p < 1$.

(A.III) The map $\cE_{\cX} : = \iota_{d_{\cX}} \circ \widetilde{\cD}_{\cX} \circ \widetilde{\cE}_{\cX}$ is measurable (here $\iota_{d_{\cX}} : \cX \rightarrow \bbR^{d_{\cX}}$ is the restriction of $\iota$, i.e., $\iota_{d_{\cX}}(X) = (\iota(X)_i)^{d_{\cX}}_{i=1}$) and the pushforward $\tilde\varsigma : = \cE_{\cX} \sharp \mu$ is absolutely continuous with respect to $\varrho$.

(A.IV) The maps $\cD_{\cY}$ and $\cE_{\cY}$ are linear and bounded.

Now let $X_1,\ldots,X_m \sim_{\mathrm{i.i.d.}} \mu$ and consider the training data

{ (X_i , Y_i ) }^m_i=1 ⊂(×)^m,  where Y_i = F(X_i) + E_i ∈
and $E_i \in \cY$ represents noise. Let $\cN$ be a class of DNNs $N : \bbR^{d_{\cX}} \rightarrow \bbR^{d_{\cY}}$, and define

F ≈F :=_ ∘N ∘_,  where N ∈N ∈ 1/m ∑^m_i=1 Y_i -_ ∘N ∘_(X_i) ^2_.

§.§ Discussion of assumptions

We now discuss (A.I)-(A.IV). In <ref> we describe future work on relaxing these assumptions.

(A.I) is a weak assumption. It asserts that there is a Lipschitz map $\iota$ under which the pushforward of $\mu$ is a quasi-uniform measure supported in a hypercube. As we discuss in Example <ref>, this is notably the case when $\mu$ is the law of some random field with an affine parametrization involving  bounded random variables. This situation occurs frequently in parametric and stochastic PDE problems.
(A.II) describes the specific holomorphy of the operator $F$. See Remark <ref> for details. Note that we require $\bm{b} \in \ell^p_{\mathsf{M}}(\bbN)$, not just $\bm{b} \in \ell^p(\bbN)$. It is known [Adcock et al., 2024] that one cannot learn holomorphic functions (and hence holomorphic operators) from finite data if $\bm{b} \in \ell^p(\bbN)$ only.
(A.III) is a relatively weak assumption. In view of (A.I), we expect it to hold as long as the $\widetilde{\cD}_{\cX} \circ \widetilde{\cE}_{\cX} \approx \cI_{\cX}$ sufficiently well.
Finally, (A.IV) is a standard assumption, which holds for instance in the case of PCA-Net and DeepONet. The former also requires the learned encoder $\widetilde{\cE}_{\cX}$ to be linear, which is not needed in our setup. Moreover, both usually deal with the case where both $\cX$ and $\cY$ are Hilbert spaces.

[Parametric PDEs]

A common operator learning problem involves learning the map

F : a ∈↦u(a) ∈,   where $u(a)$ satisfies $\cF_a u = 0$

and $\cF_{a}$ specifies a certain PDE depending on a parameter or function $a$.
A standard example is the elliptic diffusion equation over a domain $\Omega \subset \bbR^n$.
Here $a = a(\bm{x}) \in \rL^{\infty}(\Omega) = : \cX$ is the diffusion coefficient and $u = u(\cdot ; a)$ is the solution of the PDE

- ∇·( a ∇u(z ; a)) = g, z ∈Ω,  u(z ; a) = 0, z ∈∂Ω.

Problems such as PDE-operator are ubiquitous in scientific computing, with many applications engineering, biology, physics, finance and beyond. In many such applications, it is common to assume that the measure $\mu$ on $\cX$ is the law of a random field

a (x) = a(·; x) = a_0(·) + ∑^∞_i=1 c_i x_i φ_i(·),

for functions $a_0,\varphi_i \in \cX$, where the $x_i$ are random variables and $c_i \geq 0$ are scalars that ensure that $a \in \rL^{\infty}(\Omega)$. Under some mild assumptions, a-affine is then the Karhunen–Loève (KL) expansion of the measure $\mu$.
See, e.g., [Stuart, 2010] (see also <cit.>). 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 that $\bm{x} \in D = [-1,1]^{\bbN}$.

Problems of this type fits into our framework. Suppose that $\bm{x} = (x_1,x_2,\ldots) \sim \varrho$. The measure $\mu$ is then given as the pushforward $\mu = a \sharp \varrho$ and $f : D \rightarrow \cY$ is the parametric solution map $f : \bm{x} \in D \mapsto u(\cdot ; a(\bm{x})) \in \cY$. If needed, the map $\iota$ can be defined in a number of different ways. Suppose, for instance, that $\cX$ is a Hilbert space, e.g., $\cX = \rL^2(\Omega)$, and $\{ \varphi_i \}^{\infty}_{i=1}$ is a Riesz system (this holds, for instance, in the case of a KL expansion, in which case $\{ \varphi_i \}^{\infty}_{i=1}$ is an orthonormal basis). Then $\{ \varphi_i \}^{\infty}_{i=1}$ has a unique biorthogonal dual Riesz system $\{ \psi_i \}^{\infty}_{i=1}$. We may therefore define $\iota : a \mapsto \left ( \ip{a - a_0 }{\psi_i}_{\rL^2(\Omega)} /c_i \right )^{\infty}_{i=1}$. Notice that $\iota$ is a bounded linear map and $F(X) = f \circ \iota(X) = f(\bm{x})$ for $X = a(\bm{x}) \sim \mu$. However, $\iota$ is usually not required for computations (see <ref>).

This example considers an affine parametrization a-affine inducing the measure $\mu$. Other parametrizations can also be considered. Common examples include the quadratic $a(z ; \bm{x}) = a_0(z) + \left ( \sum^{\infty}_{i=1} c_i x_i \varphi_i(z) \right)^2$ and log-transformed $a(z,\bm{x}) = \exp \left ( \sum^{\infty}_{i=1} c_i x_i \varphi_i(z) \right )$ parametrizations [Chkifa et al., 2015].

[Holomorphy assumption]

In the previous example, the operator $F$ stems is the solution map $f : D \rightarrow \cY$ 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 $(\bm{b},\varepsilon)$-holomorphic. Consider, for instance, the affine diffusion problem stationary-diffusion-a-affine. Under a mild uniform ellipticity condition, the solution map of the standard weak form of the PDE $f : \bm{x} \in D \mapsto u(a(\cdot;\bm{x})) \in \rH^1_0(\Omega)$ is $(\bm{b},\varepsilon)$-holomorphic with $\bm{b} = (b_i)^{\infty}_{i=1}$ and $b_i = c_i \nm{\psi_i}_{\rL^{\infty}(\Omega)}$. See, e.g., <cit.>, as well as <ref>.

Similar holomorphy results have also been proved for many 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 [Cohen and DeVore, 2015] or <cit.> 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 $\hat{t}$ is a $\tau$-approximate minimizer for some $\tau \geq 0$ if $f(\hat{t}) \leq \min_t f(t) + \tau^2$.

[Existence of good DNN architectures]

Let $m \geq 3$, $\delta > 0$, $0 < \epsilon < 1$ and $L  = \log^4(m) + \log(1/\epsilon)$. Then there exists a class $\cN$ of hyperbolic tangent (tanh) DNNs $N : \bbR^{d_{\cX}} \rightarrow \bbR^{d_{\cY}}$ depending on $m$ and $\epsilon$ only with

width() ≲(m/L)^1+δ,  depth() ≲log(m/L),

such that following holds.
Suppose that Assumption <ref> holds and

d_ ≥⌈m/L⌉,    L_ι ·_ - _ ∘__L^2_μ(; )
≤c ·(m/L)^-1/2,

where $c > 0$ is a universal constant. Let $X_1,\ldots, X_m \sim_{\mathrm{i.i.d.}} \mu$ and consider the noisy training data DNN-training-data with arbitrary noise $E_i \in \cY$. Then, with probability at least $1- \epsilon$, every $\tau$-minimizer $\widehat{N}$ of DNN-training-prob, where $\tau \geq 0$ is arbitrary, yields an approximation $\widehat{F}$ that satisfies

F - F_L^2_μ(;)     ≲E_𝖺𝗉𝗉,2 + E_,2+ E_,2 + E_𝗈𝗉𝗍,2 + E_𝗌𝖺𝗆𝗉,2,

F - F_L^∞_μ(;)     ≲E_𝖺𝗉𝗉,∞ + E_,∞+ E_,∞ + E_𝗈𝗉𝗍,∞ + E_𝗌𝖺𝗆𝗉,∞,

and, if $\cY$ is a Hilbert space,

F - F_L^2_μ(;)  ≲E_𝖺𝗉𝗉,2 + E_,2+ E_,2 + E_𝗈𝗉𝗍,2 + E_𝗌𝖺𝗆𝗉,2.

Here, the approximation error terms $E_{\mathsf{app},q}$, $q = 2,\infty$, are given by

E_𝖺𝗉𝗉,q = a_ ·C(b,p,ξ) ·(m/L)^θ+1-1/q-1/p,

where $a_{\cY} = \nm{\cD_{\cY} \circ \cE_{\cY}}_{\cY \rightarrow \cY}$, $C(\bm{b},p,\xi) > 0$ depends on $\bm{b}$, $p$ and $\xi$ only and $\theta = 0$ if $\cY$ is a Hilbert space (as in error-bound-H) or $\theta = 1/2$ otherwise (as in error-bound-1-error-bound-2). The other terms are given by

E_, 2     = a_ ·L_ι ·√(m/(L ϵ)) ·_ - _ ∘__L^2_μ(; )

E_, ∞     =a_ ·L_ι ·√(m/L) ·( √(m/L) ·_ - _ ∘__L^2_μ(; )
+ _ - _ ∘__L^∞_μ(; ) )

E_,2     =  _ - _ ∘_ _L^2_F ♯μ(;)/√(ϵ)

E_,∞     =  _ - _ ∘_ _L^∞_F ♯μ(;) + √(m/L) ·_ - _ ∘_ _L^2_F ♯μ(;),

where $\cI_{\cX} : \cX \rightarrow \cX$, $\cI_{\cY} : \cY \rightarrow \cY$ are the identity maps and, if $\nm{\bm{E}}^2_{2;\cY} = \sum^{m}_{i=1} \nm{E_i}^2_{\cY}$,

E_𝗈𝗉𝗍 =  τ+ 2^-m     q = 2
√(m/L) τ+ 2^-m     q = ∞,  E_𝗌𝖺𝗆𝗉,q =  E_2; / √(m)     q = 2
E_2; / √(L)     q = ∞.

(Proofs of this and all other theorems are in <ref>-<ref> of the supplemental material.)
This theorem shows that there is a family of tanh DNNs which yield provable bounds for learning holomorphic operators. The error error-bound-1-error-bound-H decomposes into an approximation error Eapp-def, which decays algebraically in the amount of training data $m$. Later, in Theorems <ref>-<ref>, we show that these rates are optimal when $\cY$ is a Hilbert space, up to log factors. Next, are the encoding-decoding errors EX-def, which depend on the approximate encoders and decoders. These are standard terms in operator learning (for $q = 2$). See, e.g., [Lanthaler, 2023, Lanthaler et al., 2022]. When the encoders and decoders are computed via PCA, as in PCA-Net, standard bounds can be derived for these terms [Lanthaler, 2023].  For similar analysis in the case of DeepONets, see [Lanthaler et al., 2022]. Finally, the error error-bound-1-error-bound-H involves an optimization error $E_{\mathsf{opt}}$, which depends on $\tau$, and a sampling error $E_{\mathsf{samp}}$, which depends on the error in the training data DNN-training-data.

The bound error-bound-1 consider the Pettis $L^2$-norm. Since $\nm{\cdot}_{L^2_{\mu}(\cX ; \cY)} \leq \nm{\cdot}_{L^{\infty}_{\mu}(\cX ; \cY)}$, error-bound-2 provides a (non-sharp) estimate for the Bochner $L^2$-norm. This arises because $\cY$ is a general Banach space. If $\cY$ is a Hilbert space, then error-bound-H provides a better bound.
We also observe that DNN architecture in Theorem <ref> is independent of the smoothness of the operator being learned. We term such an architecture problem agnostic.
This theorem considers a DNN architecture using tanh activations. As we discuss in Remark <ref>, other activations can be readily used instead. Another key facet of this theorem is the width and depth bounds width-depth-bounds. 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 [De Ryck et al., 2021, de Hoop et al., 2022] and references therein). We also see this later in <ref>.

On the other hand, the family $\cN$ is not fully connected. As we describe in <ref>, 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 <ref> is essentially a theoretical contribution. In our next result, we consider the more practical scenario.

[Fully-connected DNN architectures are good]

There are universal constants $c_1,c_2,c_3,c_4 \geq 1$ such that the following holds.
Let $m$, $\delta$, $\epsilon$ and $L$ be as in Theorem <ref>,

d_ ≥c_1 (m + log(1/ϵ)),
L_ι ·_ -_ ∘__L^2_μ(; )
≤c(δ) ·(m+log(1/ϵ))^-1/2

where $c(\delta) > 0$ depends on $\delta$ only,
consider any class $\cN$ of fully-connected DNNs satisfying

(n_0,n_L+2) = (d_,d_),  n_1 , …, n_L ≥c_2 ·(m + log(1/ϵ)) ·(m/L)^δ,  L ≥c_3 ·log(m/L).

Suppose that Assumption <ref> holds and that the pushforward $\varsigma$ in (A.I) is a tensor-product of a univariate probability distribution with mean zero and variance $\omega \gtrsim 1$.  Let $X_1,\ldots, X_m \sim_{\mathrm{i.i.d.}} \mu$ and consider DNN-training-data with arbitrary $E_i \in \cY$. Then the following hold with probability at least $1-\epsilon$.

(A) Uncountably many `good' minimizers. The problem DNN-training-prob has uncountably many minimizers that satisfy error-bound-1 with $\tau = 0$ or error-bound-H with $\tau = 0$ if $\cY$ is a Hilbert space. They also satisfy error-bound-2 with $\tau = 0$ and the modified right-hand side $\sqrt{L} E_{\mathsf{app},\infty} + L E_{\cX,\infty} + \sqrt{L} E_{\cY,\infty} + E_{\mathsf{opt},\infty} + \sqrt{L} E_{\mathsf{samp},\infty}$.

(B)  Good minimizers are stable. Suppose that $\cE_{\cX} \in L^{\infty}_{\mu}(\cX ; \bbR^{d_{\cX}})$. Then, for any $\tau_o > 0$ there is a neighbourhood of DNN parameters around the parameters of each minimizer in (A) for which the approximations corresponding to the parameters in this neigbourhood also satisfy the same bounds as in (A) with $\tau = \tau_o$.

(C) Good minimizers can be far apart in parameter space. For sufficiently large $m$, there are at least $(m/(c_4 L))^{2 \delta m}$ minimizers satisfying the bounds in (A) such that, for any two such minimizers, their parameters satisfy $\nm{\bm{\theta}'} = \nm{\bm{\theta}}$ and $\nm{\bm{\theta}' - \bm{\theta}} \gtrsim 1$.

This theorem states that DL with fully-connected DNN architectures of sufficient width and depth width-depth-min can succeed, since there are minimizers that yield the optimal bounds of Theorem <ref>. 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 <ref> – but our numerical results in <ref> suggest that (approximate) minimizers obtained through training do, at least for the experiments considered, achieve the rates specified in Theorem <ref>.

§ 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(\cX ; \cY)$ be the Banach space of continuous operators. We term an (adaptive) sampling map as any map

: C(,) →^m,  F ↦(F) = ( F(X_i) )^m_i=1,

where $X_1\in \cX$, $X_2 = X_{2}(F(X_1)) \in \cX$ potentially depends on the previous evaluation $F(X_1)$, $X_3 = X_{3}(F(X_1),F(X_2)) \in \cX$, and so forth.
Next, we term a reconstruction map as any map $\cR : \cY^m \rightarrow L^2_{\mu}(\cX ; \cY)$.
Given this, we let $\cH(\bm{b} , \iota) = \{ F = f \circ \iota : f \in \cH(\bm{b}) \}$ and define

θ_m(b) = inf_,  sup_F ∈b , ι) F - ∘(F)_L^2_μ(; ) ,

where the infimum is taken over all such $\cL$ and $\cR$. In other words, $\theta_m(\bm{b})$ measures how well one can learn holomorphic operators using arbitrary training data and an arbitrary reconstruction procedure.

[Optimal $L^2$ error rates]

For any $0 < p < 1$ there is a constant $c(p) >0$ such that the following hold.

(i) For each $m \in \bbN$, there is a $\bm{b} \in \ell^p_{\mathsf{M}}(\bbN)$, $\bm{b} \geq \bm{0}$, $\nm{\bm{b}}_{p,\mathsf{M}} = 1$ such that $\theta_m(\bm{b}) \geq c(p) \cdot m^{1/2-1/p}$.

(ii) There is a $\bm{b} \in \ell^p_{\mathsf{M}}(\bbN)$, $\bm{b} \geq \bm{0}$, $\nm{\bm{b}}_{p,\mathsf{M}} = 1$ such that $\theta_m(\bm{b}) \geq c(p) \cdot \frac{m^{1/2-1/p}}{ \log^{2/p}(2m)}$, $\forall m \in \bbN$.

This theorem shows that the term $E_{\mathsf{app},2}$ in Theorems <ref>-<ref> is optimal, up to log terms whenever $\cY$ is a Hilbert space: there does not exist a reconstruction procedure that can surpass the rate $m^{1/2-1/p}$ for learning holomorphic operators. Note that result applies not only to DL-based procedures, but any procedure learns holomorphic functions from $m$ samples. Another consequence of this result is that adaptive sampling, i.e., active learning, is of no benefit. As shown by Theorems <ref>-<ref>, 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 $\mu$.

Theorem <ref> considers $L^2$-norm. For the $L^{\infty}$-norm, we present a somewhat weaker result. Let

θ̃_m(b) = inf_ { _X_1,…,X_m ∼μ sup_F ∈(b,ι) F - ({ X_i, F(X_i) })_L^∞_μ(;)  },

where the infimum is taken over all reconstruction maps $\cR : (\cX \times \cY)^{m} \rightarrow L^{\infty}_{\mu}(\cX ; \cY)$ only.

[Optimal $L^{\infty}$ error rates]

Suppose that the pushforward $\varsigma$ in (A.I) is a tensor-product of a univariate probability distribution. Then, for any $0 < p < 1$ there is a constant $c(p) >0$ such that the following hold.

(i) For each $m \in \bbN$, there is a $\bm{b} \in \ell^p_{\mathsf{M}}(\bbN)$, $\bm{b} \geq \bm{0}$, $\nm{\bm{b}}_{p,\mathsf{M}} = 1$ such that $\tilde\theta_m(\bm{b}) \geq c(p) \cdot  \frac{m^{1-1/p}}{ \log(m)}$.

(ii) There is a $\bm{b} \in \ell^p_{\mathsf{M}}(\bbN)$, $\bm{b} \geq \bm{0}$, $\nm{\bm{b}}_{p,\mathsf{M}} = 1$ such that $\tilde\theta_m(\bm{b}) \geq c(p) \cdot \frac{m^{1-1/p}}{ \log^{2/p+1}(2m)}$, $\forall m \in \bbN$.

As with the previous theorem, this result asserts that the rate $m^{1-1/p}$ is optimal in the $L^{\infty}$-norm when $\cY$ is a Hilbert space. However, it is strictly weaker than Theorem <ref> as it only considers i.i.d. random sampling from $\mu$, as opposed to arbitrary (adaptive) samples.

§ NUMERICAL EXPERIMENTS

We now present numerical results for DL applied to various different parametric PDE problems, as in Example <ref>. For a full description of our experimental setup, see <ref>-<ref>.

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 $\cE_{\cX}$ and $\cD_{\cY}$ in F-approx do not change for different choices of $\widehat{N}$. We also set up our experiments so that encoding-decoding (<ref>) and sampling (<ref>) errors are zero. We do this in a standard way. To ensure that $E_{\cX,q} = 0$, we truncate the parametric expansions a-affine after $d$ terms and define the encoder $\cE_{\cX}$ accordingly. This means we effectively consider a parametric PDE depending on finitely-many parameters. To perform our simulations we use Finite Element Methods (FEMs) to both solve the PDE (for generating training and testing data) and define the decoder $\cD_{\cY}$ (see DY-def).
To ensure that $E_{\cY,q} = 0$, we compute errors with respect to the Bochner $L^2_{\mu}(\cX ; \widetilde{\cY})$-norm, where $\widetilde{\cY} = \cD_{\cY}(\bbR^{d_{\cY}})$
is the FEM discretization of $\cY$. 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 $\widehat{F}$. See <ref> for further details.

The DNNs in our experiments are fully-connected and of the form (<ref>). We denote by $\sigma$ $L\times N$ DNN a DNN $\widehat{N}$ with activation function $\sigma$, width $N$ and depth $L$. To solve (<ref>) we use Adam [Kingma and Ba, 2017] 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 <ref> for further details.

Parametric elliptic diffusion equation. Our first example is the parametric elliptic diffusion equation stationary-diffusion. This PDE arises in many scientific computing applications, such as groundwater flow modeling, see, e.g., [Winter and Tartakovsky, 2002].
We describe the full PDE and its FE discretization in <ref>. In our experiments, we consider both affine (<ref>) and log-transformed (<ref>)
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 [Gatica, 2014] of stationary-diffusion, which has a number of key practical benefits (see <ref>). In this case, $\cY = \rL^2(\Omega)$ is a Hilbert space. Fig.  <ref> compares the error versus the amount of training data $m$ for various DNN architectures for learning the solution of this PDE in $d = 4$ and $d=8$ dimensions with these two diffusion coefficients.
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. Overall, the best performing DNNs appear to roughly match the plotted rate $m^{-1}$. As we explain further in <ref>, this rate is precisely that predicted by our theory. In particular, the parametric solution map (recall Remark <ref>) is $(\bm{b},\varepsilon)$-holomorphic with $\bm{b} \in \ell^p_{\mathsf{M}}(\bbN)$ 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 dimension of the parametric problem doubling from 4 to 8, there is little change in the error behaviour.

Elliptic diffusion equation. Average relative $L^2_{\mu}(\cX ; \widetilde{\cY})$-norm error versus $m$ for different DNNs approximating the solution operator for the elliptic diffusion equation (<ref>). The first two plots use the affine coefficient $a_{1,d}$ (<ref>) with $d=4,8$, respectively. The rest use the log-transformed coefficient $a_{2,d}$(<ref>).

NSB equations.
Average relative $L^2_{\mu}(\cX ; \widetilde{\cY})$-norm error versus $m$ for different DNNs approximating the velocity field $\bm{u}$ of the NSB problem in (<ref>). See Fig. <ref> for results for the pressure component $p$. The diffusion coefficients $a_{1,d},a_{2,d}$ and dimensions $d = 4,8$ are as in Fig. <ref>.

Boussinesq equation. Average relative $L^2_{\mu}(\cX ; \widetilde{\cY})$-norm error versus $m$ for different DNNs approximating the temperature $\varphi$ of the Boussinesq problem in (<ref>) (see Fig. <ref> for $\bm{u}$ and $p$). The diffusion coefficients $a_{1,d},a_{2,d}$ and dimensions $d = 4,8$ are as in Fig. <ref>. In this example, we also consider an additional parametric dependence in the tensor $\bbK = \bbK_d$ describing the thermal conductivity of the fluid. See <ref> and (<ref>).

Parametric Navier-Stokes-Brinkman equations.
We next consider the parametric Navier-Stokes-Brinkman (NSB) equations. See <ref> and eq:NSB for the full definition. Here the solution is a pair $(\bm{u},p)$, where $\bm{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., [Dutil et al., 2011, Khadra et al., 2000, Wang et al., 2010, Hwang and Advani, 2010]. We use a mixed variational formulation [Gatica et al., 2022] 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 $\cY$ being either $\cY = \Lft $ for $\bm{u}$ or $\cY = \rL^2(\Omega)$ for $p$. See <ref> for details.
Fig.  <ref> compares a variety of DNN architectures for approximating the velocity field component in $d=4$ and $d=8$ 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}$, and, similar to the previous example, there is no deterioration of the rate when moving from $d=4$ to $d=8$ dimensions.

Parametric stationary Boussinesq equation.
Our final example is a parametric stationary Boussinesq PDE. See <ref> and eq:BSQ for the full definition. Here the solution is a triplet $(\u,\varphi,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 [Guo et~al., 2002, Cao and Wu, 2013, Danaila et~al., 2014]. Similar to the previous example, we consider a fully mixed variational formulation (see \S \ref{ss:boussinesq-variational}), which leads to $= $ (for $u$), $= ^4(Ω)$ (for $φ$) or $= ^2_0(Ω)$ (for $p$).
Fig.\  \ref{Boussinesq_res} provides numerical results. Our observations are in line with the previous two examples where both ELU and the smaller tanh networks are most often the best performers in this problem. Once more, the errors roughly correspond to the rate $m^-1$.

\section{Conclusions and limitations}\label{s:conclusions}

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 [Lanthaler et~al., 2022, Franco et~al., 2023, Marcati and Schwab, 2023, Herrmann et~al., 2022, Deng et~al., 2022] 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 \textit{curse of dimensionality}. It is notable that size of the DNNs in Theorems \ref{thm:main-res-1}-\ref{thm:main-res-2} also do not succumb to the so-called \textit{curse of parametric complexity} [Lanthaler and Stuart, 2023, Lanthaler and Stuart, 2024], 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 $ς$ is Gaussian and supported on $^$. Second, it is unknown whether (A.III) can be relaxed. We believe this can be done, at least if the $L^2_μ$-norm in \ef{enc-dec-err} is replaced by the $L^∞_μ$-norm. Whether this is possible without modifying \ef{enc-dec-err} is currently unknown.

Third, a limitation of Theorem \ref{thm:main-res-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 \S \ref{ss:previous-work}, 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, when $$ is a Banach space the approximation error terms \ef{Eapp-def} are worse by a factor of $1/2$ than when $$ is a Hilbert space. This arises because of the loss of the inner product structure and, in particular, the absence of Parseval's identity. Our numerical results in Figs.\ \ref{NSB_res}-\ref{Boussinesq_res}  suggest that this factor is an artefact of the proofs. Whether it can be removed is an interesting open problem.

\newpage
\bibliographystyle{abbrvnat}
\begin{thebibliography}{107}
\providecommand{\natexlab}[1]{#1}
\providecommand{\url}[1]{\texttt{#1}}
\expandafter\ifx\csname urlstyle\endcsname\relax
\providecommand{\doi}[1]{doi: #1}\else
\providecommand{\doi}{doi: \begingroup \urlstyle{rm}\Url}\fi

[Adcock and Dexter, 2021]
B.~Adcock and N.~Dexter.
\newblock The gap between theory and practice in function approximation with
deep neural networks.
\newblock \emph{SIAM J. Math. Data Sci.}, 3\penalty0 (2):\penalty0 624--655,

[Adcock et~al., 2021]
B.~Adcock, S.~Brugiapaglia, N.~Dexter, and S.~Moraga.
\newblock Deep neural networks are effective at learning high-dimensional
{H}ilbert-valued functions from limited data.
\newblock In J.~Bruna, J.~S. Hesthaven, and L.~Zdeborov{\'a}, editors,
\emph{Proceedings of The Second Annual Conference on Mathematical and
Scientific Machine Learning}, volume 145 of \emph{Proc. Mach. Learn. Res.
(PMLR)}, pages 1--36. PMLR, 2021.

[Adcock et~al., 2022]
B.~Adcock, S.~Brugiapaglia, N.~Dexter, and S.~Moraga.
\newblock {Near-optimal learning of Banach-valued, high-dimensional functions
via deep neural networks}.
\newblock \emph{arXiv preprint arXiv:2211.12633}, 2022{\natexlab{a}}.

[Adcock et~al., 2022]
B.~Adcock, S.~Brugiapaglia, and C.~G. {Webster}.
\newblock \emph{Sparse Polynomial Approximation of High-Dimensional Functions}.
\newblock Comput. Sci. Eng. Society for Industrial and Applied Mathematics,
Philadelphia, PA, 2022{\natexlab{b}}.

[Adcock et~al., 2023]
B.~Adcock, N.~Dexter, and S.~Moraga.
\newblock Optimal approximation of infinite-dimensional holomorphic functions
{II}: recovery from i.i.d. pointwise samples.
\newblock \emph{arXiv:2310.16940}, 2023.

[Adcock et~al., 2024]
B.~Adcock, S.~Brugiapaglia, N.~Dexter, and S.~Moraga.
\newblock Learning smooth functions in high dimensions: from sparse polynomials
to deep neural networks.
\newblock \emph{arXiv:2404.03761}, 2024{\natexlab{a}}.

[Adcock et~al., 2024]
B.~Adcock, N.~Dexter, and S.~Moraga.
\newblock Optimal approximation of infinite-dimensional holomorphic functions.
\newblock \emph{Calcolo}, 61\penalty0 (1):\penalty0 12, 2024{\natexlab{b}}.

[Aliprantis and Border, 2006]
C.~D. Aliprantis and K.~C. Border.
\newblock \emph{Infinite Dimensional Analysis: A Hitchhiker's Guide}.
\newblock Springer--Verlag, Berlin, Heidelberg, 3rd edition, 2006.

[Aln{\ae}s et~al., 2015]
S.~Aln{\ae}s, J.~Blechta, J.~Hake, A.~Johansson, B.~Kehlet, A.~Logg,
C.~Richardson, J.~Ring, M.~E. Rognes, and G.~N. Wells.
\newblock {The FEniCS Project Version 1.5}.
\newblock \emph{Archive of Numerical Software}, 3\penalty0 (100), 2015.

[Bhattacharya et~al., 2021]
K.~Bhattacharya, B.~Hosseini, N.~B. Kovachki, and A.~M. Stuart.
\newblock {Model reduction and neural networks for parametric PDEs}.
\newblock \emph{The SMAI journal of computational mathematics}, 7:\penalty0
121--157, 2021.

[Boffi et~al., 2013]
D.~Boffi, F.~Brezzi, and M.~Fortin.
\newblock \emph{Mixed Finite Element Methods and Applications}.
\newblock Springer Berlin, Heidelberg, 1 edition, 2013.
\newblock \doi{https://doi.org/10.1007/978-3-642-36519-5}.

[Brugiapaglia et~al., 2021]
S.~Brugiapaglia, S.~Dirksen, H.~C. Jung, and H.~Rauhut.
\newblock Sparse recovery in bounded {R}iesz systems with applications to
numerical methods for {PDE}s.
\newblock \emph{Appl. Comput. Harmon. Anal.}, 53:\penalty0 231--269, 2021.

[Cai et~al., 2021]
S.~Cai, Z.~Wang, L.~Lu, T.~A. Zaki, and G.~E. Karniadakis.
\newblock {DeepM\&Mnet: Inferring the electroconvection multiphysics fields
based on operator approximation by neural networks}.
\newblock \emph{Journal of Computational Physics}, 436:\penalty0 110296, 2021.

[Cao and Wu, 2013]
C.~Cao and J.~Wu.
\newblock {{{Global regularity for the two-dimensional anisotropic Boussinesq
equations with vertical dissipation}}}.
\newblock \emph{Arch. Rational Mech. Anal.}, 208:\penalty0 985--1004, 2013.
\newblock \doi{https://doi.org/10.1007/s00205-013-0610-3}.

[Cao et~al., 2023]
Q.~Cao, S.~Goswami, and G.~E. Karniadakis.
\newblock {LNO: Laplace neural operator for solving differential equations}.
\newblock \emph{arXiv preprint arXiv:2303.10528}, 2023.

[Chen et~al., 2023]
K.~Chen, C.~Wang, and H.~Yang.
\newblock Deep operator learning lessens the curse of dimensionality for pdes.
\newblock \emph{arXiv preprint arXiv:2301.12227}, 2023.

[Chkifa et~al., 2015]
A.~Chkifa, A.~Cohen, and C.~Schwab.
\newblock {Breaking the curse of dimensionality in sparse polynomial
approximation of parametric PDEs}.
\newblock \emph{J. Math. Pures Appl.}, 103\penalty0 (2):\penalty0 400--428,

[Ciarlet, 2002]
P.~G. Ciarlet.
\newblock \emph{The Finite Element Method for Elliptic Problems}.
\newblock SIAM, 2002.

[Cohen and DeVore, 2015]
A.~Cohen and R.~A. DeVore.
\newblock Approximation of high-dimensional parametric {PDE}s.
\newblock \emph{Acta Numer.}, 24:\penalty0 1--159, 2015.

[Colmenares et~al., 2020]
E.~Colmenares, G.~N. Gatica, and S.~Moraga.
\newblock {{{A Banach spaces-based analysis of a new fully-mixed finite element
method for the Boussinesq problem}}}.
\newblock \emph{ESAIM Math. Model. Numer. Anal.}, 54\penalty0 (5):\penalty0
1525--1568, 2020.

[D\~{u}ng}}} et~al., 2023]
D.~{{{D\~{u}ng}}}, V.~K. Nguyen, and D.~T. Pham.
\newblock {Deep ReLU neural network approximation in Bochner spaces and
applications to parametric PDEs}.
\newblock \emph{J. Complexity}, 79:\penalty0 101779, 2023.

[Danaila et~al., 2014]
I.~Danaila, R.~Moglan, F.~Hecht, and S.~Le~Masson.
\newblock A newton method with adaptive finite elements for solving
phase-change problems with natural convection.
\newblock \emph{J. Comput. Phys.}, 274:\penalty0 826--840, 2014.

[Daws and Webster, 2019]
J.~Daws and C.~Webster.
\newblock Analysis of deep neural networks with quasi-optimal polynomial
approximation rates.
\newblock \emph{arXiv:1912.02302}, 2019.

[de~Hoop et~al., 2022]
M.~V. de~Hoop, D.~Z. Huang, E.~Qian, and A.~M. Stuart.
\newblock The cost-accuracy trade-off in operator learning with neural
\newblock \emph{arXiv preprint arXiv:2203.13181}, 2022.

[De~Ryck et~al., 2021]
T.~De~Ryck, S.~Lanthaler, and S.~Mishra.
\newblock On the approximation of functions by tanh neural networks.
\newblock \emph{Neural Networks}, 143:\penalty0 732--750, 2021.
\newblock ISSN 0893-6080.

[Deng et~al., 2022]
B.~Deng, Y.~Shin, L.~Lu, Z.~Zhang, and G.~E. Karniadakis.
\newblock {Approximation rates of DeepONets for learning operators arising from
advection--diffusion equations}.
\newblock \emph{Neural Networks}, 153:\penalty0 411--426, 2022.

[Dexter et~al., 2019]
N.~Dexter, H.~Tran, and C.~Webster.
\newblock {{{A mixed $\ell_1$ regularization approach for sparse simultaneous
approximation of parameterized PDEs}}}.
\newblock \emph{ESAIM Math. Model. Numer. Anal.}, 53:\penalty0 2025--2045,

[Dutil et~al., 2011]
Y.~Dutil, D.~R. Rousse, N.~B. Salah, S.~Lassue, and L.~Zalewski.
\newblock A review on phase-change materials: Mathematical modeling and
\newblock \emph{Renew. Sustain. Energy Rev.}, 15\penalty0 (1):\penalty0
112--130, 2011.

[Foucart and Rauhut, 2013]
S.~Foucart and H.~Rauhut.
\newblock \emph{A Mathematical Introduction to Compressive Sensing}.
\newblock Appl. Numer. Harmon. Anal. Birkh{\"a}user, New York, NY, 2013.

[Franco and Brugiapaglia, 2024]
N.~R. Franco and S.~Brugiapaglia.
\newblock A practical existence theorem for reduced order models based on
convolutional autoencoders.
\newblock \emph{arXiv:2402.00435}, 2024.

[Franco et~al., 2023]
N.~R. Franco, S.~Fresca, A.~Manzoni, and P.~Zunino.
\newblock {Approximation bounds for convolutional neural networks in operator
\newblock \emph{Neural Networks}, 161:\penalty0 129--141, 2023.

[Gatica, 2014]
G.~N. Gatica.
\newblock \emph{{{{A Simple Introduction to the Mixed Finite Element
\newblock Springer Cham, 2014.

[Gatica et~al., 2021]
G.~N. Gatica, R.~Oyarz\'ua, R.~Ruiz-Baier, and Y.~D. Sobral.
\newblock Banach spaces-based analysis of a fully-mixed finite element method
for the steady-state model of fluidized beds.
\newblock \emph{Comput. Math. Appl.}, 84:\penalty0 244--276, 2021.

[Gatica et~al., 2022]
G.~N. Gatica, N.~Nu\~{n}ez, and R.~Ruiz-Baier.
\newblock New non-augmented mixed finite element methods for the
navier--stokes--brinkman equations using banach spaces.
\newblock \emph{In review}, 2022.

[Goswami et~al., 2022]
S.~Goswami, K.~Kontolati, M.~D. Shields, and G.~E. Karniadakis.
\newblock {Deep transfer operator learning for partial differential equations
under conditional shift}.
\newblock \emph{Nature Machine Intelligence}, 4\penalty0 (12):\penalty0
1155--1164, 2022{\natexlab{a}}.

[Goswami et~al., 2022]
S.~Goswami, D.~S. Li, B.~V. Rego, M.~Latorre, J.~D. Humphrey, and G.~E.
\newblock {Neural operator learning of heterogeneous mechanobiological insults
contributing to aortic aneurysms}.
\newblock \emph{Journal of the Royal Society Interface}, 19\penalty0
(193):\penalty0 20220410, 2022{\natexlab{b}}.

[Goswami et~al., 2023]
S.~Goswami, A.~D. Jagtap, H.~Babaee, B.~T. Susi, and G.~E. Karniadakis.
\newblock Learning stiff chemical kinetics using extended deep neural
\newblock \emph{arXiv preprint arXiv:2302.12645}, 2023.

[G\"{u}hring}}} and Raslan, 2021]
I.~{{{G\"{u}hring}}} and M.~Raslan.
\newblock Approximation rates for neural networks with encodable weights in
smoothness spaces.
\newblock \emph{Neural Networks}, 134:\penalty0 107--130, 2021.

[Guo et~al., 2002]
Z.~Guo, B.~Shi, and C.~Zheng.
\newblock {{{A coupled lattice BGK model for the Boussinesq equations}}}.
\newblock \emph{Int. J. Numer. Meth. Fluids}, 39:\penalty0 325--342, 2002.
\newblock \doi{10.1002/d.337}.

[Hao et~al., 2023]
Z.~Hao, Z.~Wang, H.~Su, C.~Ying, Y.~Dong, S.~Liu, Z.~Cheng, J.~Song, and
\newblock {GNOT: A General Neural Operator Transformer for Operator Learning}.
\newblock In \emph{International Conference on Machine Learning}, pages
12556--12569. PMLR, 2023.

[Herrmann et~al., 2022]
L.~Herrmann, C.~Schwab, and J.~Zech.
\newblock Neural and gpc operator surrogates: construction and expression rate
\newblock \emph{arXiv preprint arXiv:2207.04950}, 2022.

[Howard et~al., 2022]
A.~A. Howard, M.~Perego, G.~E. Karniadakis, and P.~Stinis.
\newblock Multifidelity deep operator networks.
\newblock \emph{arXiv preprint arXiv:2204.09157}, 2022.

[Hwang and Advani, 2010]
W.~R. Hwang and S.~G. Advani.
\newblock {{{Numerical simulations of Stokes--Brinkman equations for
permeability prediction of dual scale fibrous porous media}}}.
\newblock \emph{Phys. Fluids.}, 22:\penalty0 113101, 2010.

[Hyt\"{o}nen et~al., 2016]
T.~Hyt\"{o}nen, J.~van Neerven, M.~Veraar, and L.~Weis.
\newblock \emph{Analysis in Banach Spaces, Volume I: Martingales and
Littlewood--Paley Theory}.
\newblock Ergebnisse der Mathematik und ihrer Grenzgebiete. 3. Folge. Springer,
Cham, 2016.

[Karniadakis et~al., 2021]
G.~E. Karniadakis, I.~G. Kevrekidis, L.~Lu, P.~Perdikaris, S.~Wang, and
\newblock Physics-informed machine learning.
\newblock \emph{Nature Reviews Physics}, 3\penalty0 (6):\penalty0 422--440,

[Khadra et~al., 2000]
K.~Khadra, P.~Angot, S.~Parneix, and J.~P. Caltagirone.
\newblock {{{Fictitious domain approach for numerical modelling of
Navier-Stokes equations}}}.
\newblock \emph{Int. J. Numer. Meth. Fluids}, 34:\penalty0 651--684, 2000.

[Kingma and Ba, 2017]
D.~P. Kingma and J.~Ba.
\newblock Adam: a method for stochastic optimization.
\newblock \emph{arXiv:1412.6980}, 2017.

[Kontolati et~al., 2023]
K.~Kontolati, S.~Goswami, M.~D. Shields, and G.~E. Karniadakis.
\newblock On the influence of over-parameterization in manifold based
surrogates and deep neural operators.
\newblock \emph{Journal of Computational Physics}, 479:\penalty0 112008, 2023.

[Koric et~al., 2023]
S.~Koric, A.~Viswantah, D.~W. Abueidda, N.~A. Sobh, and K.~Khan.
\newblock Deep learning operator network for plastic deformation with variable
loads and material properties.
\newblock \emph{Engineering with Computers}, pages 1--13, 2023.

[Kovachki et~al., 2021]
N.~Kovachki, S.~Lanthaler, and S.~Mishra.
\newblock {On universal approximation and error bounds for Fourier neural
\newblock \emph{The Journal of Machine Learning Research}, 22\penalty0
(1):\penalty0 13237--13312, 2021{\natexlab{a}}.

[Kovachki et~al., 2021]
N.~Kovachki, Z.~Li, B.~Liu, K.~Azizzadenesheli, K.~Bhattacharya, A.~Stuart, and
\newblock {Neural operator: Learning maps between function spaces}.
\newblock \emph{arXiv preprint arXiv:2108.08481}, 2021{\natexlab{b}}.

[Kovachki et~al., 2023]
N.~B. Kovachki, Z.~Li, B.~Liu, K.~Azizzadenesheli, K.~Bhattacharya, A.~M.
Stuart, and A.~Anandkumar.
\newblock {Neural Operator: Learning Maps Between Function Spaces With
Applications to PDEs.}
\newblock \emph{J. Mach. Learn. Res.}, 24\penalty0 (89):\penalty0 1--97, 2023.

[Kovachki et~al., 2024]
N.~B. Kovachki, S.~Lanthaler, and A.~M. Stuart.
\newblock Operator learning: algorithms and analysis.
\newblock \emph{arXiv:2402.15715}, 2024.

[Lanthaler, 2023]
\newblock {Operator learning with PCA-Net: upper and lower complexity bounds}.
\newblock \emph{arXiv preprint arXiv:2303.16317}, 2023.

[Lanthaler and Stuart, 2023]
S.~Lanthaler and A.~M. Stuart.
\newblock The curse of dimensionality in operator learning.
\newblock \emph{arXiv preprint arXiv:2306.15924}, 2023.

[Lanthaler and Stuart, 2024]
S.~Lanthaler and A.~M. Stuart.
\newblock The parametric complexity of operator learning.
\newblock \emph{arXiv:2306.15924}, 2024.

[Lanthaler et~al., 2022]
S.~Lanthaler, S.~Mishra, and G.~E. Karniadakis.
\newblock {Error estimates for DeepONets: A deep learning framework in infinite
\newblock \emph{Transactions of Mathematics and Its Applications}, 6\penalty0
(1):\penalty0 tnac001, 2022.

[Lanthaler et~al., 2023]
S.~Lanthaler, Z.~Li, and A.~M. Stuart.
\newblock {The Nonlocal Neural Operator: Universal Approximation}.
\newblock \emph{arXiv preprint arXiv:2304.13221}, 2023.

[Li et~al., 2020]
B.~Li, S.~Tang, and H.~Yu.
\newblock {{{Better approximations of high dimensional smooth functions by deep
neural networks with rectified power units}}}.
\newblock \emph{Commun. Comput. Phys.}, 27:\penalty0 379--411,

[Li et~al., 2020]
Z.~Li, N.~Kovachki, K.~Azizzadenesheli, B.~Liu, K.~Bhattacharya, A.~Stuart, and
\newblock Fourier neural operator for parametric partial differential
\newblock \emph{arXiv preprint arXiv:2010.08895}, 2020{\natexlab{b}}.

[Li et~al., 2020]
Z.~Li, N.~Kovachki, K.~Azizzadenesheli, B.~Liu, K.~Bhattacharya, A.~Stuart, and
\newblock {Neural operator: Graph kernel network for partial differential
\newblock \emph{arXiv preprint arXiv:2003.03485}, 2020{\natexlab{c}}.

[Li et~al., 2020]
Z.~Li, N.~Kovachki, K.~Azizzadenesheli, B.~Liu, A.~Stuart, K.~Bhattacharya, and
\newblock Multipole graph neural operator for parametric partial differential
\newblock \emph{Advances in Neural Information Processing Systems},
33:\penalty0 6755--6766, 2020{\natexlab{d}}.

[Li et~al., 2022]
Z.~Li, M.~Liu-Schiaffini, N.~Kovachki, K.~Azizzadenesheli, B.~Liu,
K.~Bhattacharya, A.~Stuart, and A.~Anandkumar.
\newblock Learning chaotic dynamics in dissipative systems.
\newblock \emph{Advances in Neural Information Processing Systems},
35:\penalty0 16768--16781, 2022.

[Li et~al., 2023]
Z.~Li, N.~B. Kovachki, C.~Choy, B.~Li, J.~Kossaifi, S.~P. Otta, M.~A. Nabian,
M.~Stadler, C.~Hundt, K.~Azizzadenesheli, and A.~Anandkumar.
\newblock {Geometry-Informed Neural Operator for Large-Scale 3D PDEs}.
\newblock \emph{arXiv preprint arXiv:2309.00583}, 2023.

[Liang and Srikant, 2017]
S.~Liang and R.~Srikant.
\newblock Why deep neural networks for function approximation?
\newblock In \emph{ICLR}, 2017.

[Lin et~al., 2021]
C.~Lin, Z.~Li, L.~Lu, S.~Cai, M.~Maxey, and G.~E. Karniadakis.
\newblock Operator learning for predicting multiscale bubble growth dynamics.
\newblock \emph{The Journal of Chemical Physics}, 154\penalty0 (10),

[Lin et~al., 2021]
C.~Lin, M.~Maxey, Z.~Li, and G.~E. Karniadakis.
\newblock A seamless multiscale operator neural network for inferring bubble
\newblock \emph{Journal of Fluid Mechanics}, 929:\penalty0 A18,

[Liu et~al., 2024]
H.~Liu, H.~Yang, M.~Chen, T.~Zhao, and W.~Liao.
\newblock Deep nonparametric estimation of operators between infinite
dimensional spaces.
\newblock \emph{J. Mach. Learn. Res.}, 25:\penalty0 1--67, 2024.

[Lu et~al., 2021]
J.~Lu, Z.~Shen, H.~Yang, and S.~Zhang.
\newblock Deep network approximation for smooth functions.
\newblock \emph{SIAM J. Math. Anal.}, 53\penalty0 (5):\penalty0 5465--5506,

[Lu et~al., 2019]
L.~Lu, P.~Jin, and G.~E. Karniadakis.
\newblock {DeepONet: Learning nonlinear operators for identifying differential
equations based on the universal approximation theorem of operators}.
\newblock \emph{arXiv preprint arXiv:1910.03193}, 2019.

[Lu et~al., 2021]
L.~Lu, P.~Jin, G.~Pang, Z.~Zhang, and G.~E. Karniadakis.
\newblock {Learning nonlinear operators via DeepONet based on the universal
approximation theorem of operators}.
\newblock \emph{Nature machine intelligence}, 3\penalty0 (3):\penalty0
218--229, 2021{\natexlab{b}}.

[Lu et~al., 2022]
L.~Lu, X.~Meng, S.~Cai, Z.~Mao, S.~Goswami, Z.~Zhang, and G.~E. Karniadakis.
\newblock A comprehensive and fair comparison of two neural operators (with
practical extensions) based on fair data.
\newblock \emph{Computer Methods in Applied Mechanics and Engineering},
393:\penalty0 114778, 2022.

[Marcati and Schwab, 2023]
C.~Marcati and C.~Schwab.
\newblock Exponential convergence of deep operator networks for elliptic
partial differential equations.
\newblock \emph{SIAM Journal on Numerical Analysis}, 61\penalty0 (3):\penalty0
1513--1545, 2023.

[Mhaskar, 1996]
\newblock Neural networks for optimal approximation of smooth and analytic
\newblock \emph{Neural Comput.}, 8\penalty0 (1):\penalty0 164--177, 1996.

[Micha{\l}owska et~al., 2023]
K.~Micha{\l}owska, S.~Goswami, G.~E. Karniadakis, and S.~Riemer-S{\o}rensen.
\newblock {Neural operator learning for long-time integration in dynamical
systems with recurrent neural networks}.
\newblock \emph{arXiv preprint arXiv:2303.02243}, 2023.

[Moraga, 2024]
\newblock \emph{Optimal and efficient algorithms for learning high-dimensional,
Banach-valued functions from limited samples}.
\newblock PhD thesis, Simon Fraser University (In preparation), 2024.

[Nguyen and Vu, 2018]
H.~H. Nguyen and V.~H. Vu.
\newblock Normal vector of a random hyperplane.
\newblock \emph{Int. Math. Res. Not. IMRN}, 2018\penalty0 (6):\penalty0
1754--1778, 2018.

[Nobile et~al., 2008]
F.~Nobile, R.~Tempone, and C.~G. Webster.
\newblock A sparse grid stochastic collocation method for partial differential
equations with random input data.
\newblock \emph{SIAM J. Numer. Anal.}, 46\penalty0 (5):\penalty0 2309--2345,

[Oommen et~al., 2022]
V.~Oommen, K.~Shukla, S.~Goswami, R.~Dingreville, and G.~E. Karniadakis.
\newblock Learning two-phase microstructure evolution using neural operators
and autoencoder architectures.
\newblock \emph{npj Computational Materials}, 8\penalty0 (1):\penalty0 190,

[Opschoor and Schwab, 2023]
J.~A.~A. Opschoor and C.~Schwab.
\newblock {Deep ReLU networks and high-order finite element methods II:
Chebyshev emulation}.
\newblock \emph{arXiv:2310.07261}, 2023.

[Opschoor et~al., 2022]
J.~A.~A. Opschoor, C.~Schwab, and J.~Zech.
\newblock {{{Exponential ReLU DNN expression of holomorphic maps in high
\newblock \emph{Constr. Approx.}, 55\penalty0 (1):\penalty0 537--582, 2022.
\newblock \doi{https://doi.org/10.1007/s00365-021-09542-5}.

[Osorio et~al., 2022]
J.~D. Osorio, Z.~Wang, G.~Karniadakis, S.~Cai, C.~Chryssostomidis, M.~Panwar,
and R.~Hovsapian.
\newblock Forecasting solar-thermal systems performance under transient
operation using a data-driven machine learning approach based on the deep
operator network architecture.
\newblock \emph{Energy Conversion and Management}, 252:\penalty0 115063, 2022.

[Peng et~al., 2023]
W.~Peng, Z.~Yuan, Z.~Li, and J.~Wang.
\newblock Linear attention coupled fourier neural operator for simulation of
three-dimensional turbulence.
\newblock \emph{Physics of Fluids}, 35\penalty0 (1), 2023.

[Pinkus, 1968]
\newblock \emph{$N$-widths in Approximation Theory}.
\newblock Springer--Verlag, Berlin, 1968.

[Raoni{\'c} et~al., 2023]
B.~Raoni{\'c}, R.~Molinaro, T.~Rohner, S.~Mishra, and E.~de~Bezenac.
\newblock Convolutional neural operators.
\newblock \emph{arXiv preprint arXiv:2302.01178}, 2023.

[Renn et~al., 2023]
P.~I. Renn, C.~Wang, S.~Lale, Z.~Li, A.~Anandkumar, and M.~Gharib.
\newblock Forecasting subcritical cylinder wakes with fourier neural operators.
\newblock \emph{arXiv preprint arXiv:2301.08290}, 2023.

[Schwab and Zech, 2019]
C.~Schwab and J.~Zech.
\newblock Deep learning in high dimension: neural network expression rates for
generalized polynomial chaos expansions in {UQ}.
\newblock \emph{Anal. Appl. (Singap.)}, 17\penalty0 (1):\penalty0 19--55, 2019.

[Schwab and Zech, 2023]
C.~Schwab and J.~Zech.
\newblock Deep learning in high dimension: neural network expression rates for
analytic functions in {$L^2(\mathbb{R}^d,\gamma_d)$}.
\newblock \emph{SIAM/ASA J. Uncertain. Quantif.}, 11\penalty0 (1):\penalty0
199--234, 2023.

[Schwab et~al., 2023]
C.~Schwab, A.~Stein, and J.~Zech.
\newblock {Deep Operator Network Approximation Rates for Lipschitz Operators}.
\newblock \emph{arXiv preprint arXiv:2307.09835}, 2023.

[Stoyanov, 2015]
\newblock User manual: Tasmanian sparse grids.
\newblock Technical Report ORNL/TM-2015/596, Oak Ridge National Laboratory, One
Bethel Valley Road, Oak Ridge, TN, 2015.

[Stoyanov et~al., 2013]
M.~Stoyanov, D.~Lebrun-Grandie, J.~Burkardt, and D.~Munster.
\newblock Tasmanian, 9 2013.
\newblock URL \url{https://github.com/ORNL/Tasmanian}.

[Stoyanov and Webster, 2016]
M.~K. Stoyanov and C.~G. Webster.
\newblock {A dynamically adaptive sparse grid method for quasi-optimal
interpolation of multidimensional functions}.
\newblock \emph{Comput. Math. Appl.}, 71\penalty0 (11):\penalty0 2449--2465,

[Stuart, 2010]
A.~M. Stuart.
\newblock Inverse problems: a {B}ayesian perspective.
\newblock \emph{Acta Numer.}, 19:\penalty0 451--559, 2010.

[Tang et~al., 2019]
S.~Tang, B.~Li, and Y.~Haijun.
\newblock {ChebNet}: efficient and stable constructions of deep neural networks
with rectified power units via {C}hebyshev approximation.
\newblock \emph{arXiv:1911.05467}, 2019.

[Tripura and Chakraborty, 2023]
T.~Tripura and S.~Chakraborty.
\newblock Wavelet neural operator for solving parametric partial differential
equations in computational mechanics problems.
\newblock \emph{Computer Methods in Applied Mechanics and Engineering},
404:\penalty0 115783, 2023.

[Wang et~al., 2010]
S.~Wang, A.~Faghri, and T.~L. Bergman.
\newblock A comprehensive numerical model for melting with natural convection.
\newblock \emph{Int. J. Heat Mass Transfer.}, 53\penalty0 (9-10):\penalty0
1986--200, 2010.

[Winter and Tartakovsky, 2002]
C.~L. Winter and D.~M. Tartakovsky.
\newblock Groundwater flow in heterogeneous composite aquifers.
\newblock \emph{Water Resources Research}, 38\penalty0 (8), 2002.

[Wu et~al., 2023]
K.~Wu, T.~O'Leary-Roseberry, P.~Chen, and O.~Ghattas.
\newblock {Large-scale Bayesian optimal experimental design with
derivative-informed projected neural network}.
\newblock \emph{Journal of Scientific Computing}, 95\penalty0 (1):\penalty0 30,

[Yarotsky, 2017]
\newblock {Error bounds for approximations with deep ReLU networks}.
\newblock \emph{Neural Networks}, 94:\penalty0 103--114, 2017.

[Yeung et~al., 2019]
E.~Yeung, S.~Kundu, and N.~Hodas.
\newblock {Learning deep neural network representations for Koopman operators
of nonlinear dynamical systems}.
\newblock In \emph{2019 American Control Conference (ACC)}, pages 4832--4839.
IEEE, 2019.

[Yin et~al., 2022]
M.~Yin, E.~Ban, B.~V. Rego, E.~Zhang, C.~Cavinato, J.~D. Humphrey, and
\newblock {Simulating progressive intramural damage leading to aortic
dissection using DeepONet: an operator--regression neural network}.
\newblock \emph{Journal of the Royal Society Interface}, 19\penalty0
(187):\penalty0 20210670, 2022{\natexlab{a}}.

[Yin et~al., 2022]
M.~Yin, E.~Zhang, Y.~Yu, and G.~E. Karniadakis.
\newblock Interfacing finite elements with deep neural operators for fast
multiscale modeling of mechanics problems.
\newblock \emph{Computer methods in applied mechanics and engineering},
402:\penalty0 115027, 2022{\natexlab{b}}.

[Zappala et~al., 2023]
E.~Zappala, A.~H. d.~O. Fonseca, A.~H. Moberly, M.~J. Higley, C.~Abdallah,
J.~A. Cardin, and D.~van Dijk.
\newblock Neural integro-differential equations.
\newblock In \emph{Proceedings of the AAAI Conference on Artificial
Intelligence}, volume~37, pages 11104--11112, 2023.

[Zech, 2018]
\newblock \emph{{{{Sparse-grid approximation of high-dimensional parametric
\newblock PhD thesis, ETH Zurich, 2018.

[Zhang et~al., 2023]
Z.~Zhang, L.~Wing~Tat, and H.~Schaeffer.
\newblock Belnet: basis enhanced learning, a mesh-free neural operator.
\newblock \emph{Proceedings of the Royal Society A}, 479\penalty0
(2276):\penalty0 20230043, 2023.

[Zhu et~al., 2023]
M.~Zhu, H.~Zhang, A.~Jiao, G.~E. Karniadakis, and L.~Lu.
\newblock Reliable extrapolation of deep neural operators informed by physics
or sparse observations.
\newblock \emph{Computer Methods in Applied Mechanics and Engineering},
412:\penalty0 116064, 2023.

[Zou et~al., 2022]
Z.~Zou, X.~Meng, A.~F. Psaros, and G.~E. Karniadakis.
\newblock {NeuralUQ: A comprehensive library for uncertainty quantification in
neural differential equations and operators}.
\newblock \emph{arXiv preprint arXiv:2208.11866}, 2022.

\end{thebibliography}

\newpage
\appendix

\section{Experimental setup}\label{s:exp-setup}

In this section, we describe our experimental setup.

\subsection{Formulation of the learning problems}\label{ss:num-exp-form}

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 \ef{ex:parametric-DE} and involve operators of the form \ef{PDE-operator}, where $_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 \ef{a-affine} and log-transformed parametrizations of the random field $a(x)$. Thus, in general we write
\be{
\label{a-def-general}
a(\bm{x}) = a(\cdot ; \bm{x}) = g \left ( a_0(\cdot) + \sum^{\infty}_{i=1} c_i x_i \varphi_i(\cdot) \right ),
where $g : →$ 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 dimension $d$ and truncate the expansion in \ef{a-def-general} after $d$ terms, giving a map $a_d : [-1,1]^d →$ and measure $μ= a_d ♯ϱ_d$, where $ϱ_d$ is the uniform probability measure on $[-1,1]^d$. We then define the operator $F$ as $F(a_d(x)) = f(x) = u(·; a_d(x))$, where $u(·; a)$ is the solution of the PDE $_a u = 0$.

In alignment with our theorems, we focus on the \textit{in-distribution} performance of the learned approximation $F$. This means we define the encoder only on $supp(μ)$, as $_(X) = x$ when $X = a_d(x)$
with $x ∈[-1,1]^d$. As a result, the encoding-decoding error $E_,q$ in \ef{EX-def} satisfies $E_,q = 0$.

To perform our simulations, we use FEMs to solve the PDE and discretize the output space $$. Let ${ φ_i }^K_i=1 ⊂$ be a FEM basis and set $d_ = K$. Then we define the decoder as
\be{
\label{DY-def}
\cD_{\cY} : \bbR^{K} \rightarrow \cY,\quad \cD_{\cY}(\bm{c}) = \sum^{K}_{i=1} c_i \varphi_i
and set $ = _(^K)$ as the discretization of $$.

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 parameter $a_d(x_i)$. Notice that $Y_i ∈$ 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 approximation $F = _ ∘N ∘_$ over $(μ)$ as $F = _ ∘N(X) = _ ∘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 error with respect to the $L^2_μ(; )$-norm. This effectively means that the encoding-decoding error $E_,q$ in \ef{EX-def} satisfies $E_,q = 0$ as well.

\subsection{Computational setup for the numerical experiments}
\label{app:setup_num}

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.

\begin{itemize}

\item[(i)] \textbf{Implementation}.   We use the
open-source finite element library \texttt{FEniCS}, specifically version 2019.1.0 [Aln{\ae}s et~al., 2015], and Google's \texttt{TensorFlow} version 2.12.0. More information about \texttt{TensorFlow} can be found at \url{https://www.tensorflow.org/}.

\item[(ii)] \textbf{Hardware}.
We train the DNN models in single precision on a large cluster, using Intel Xenon Processor E5-2683 v4 CPUs with either 125GB or 250GB per node. The setup for each of our PDEs is as follows. For each experiment we consider training with 14 sets of points of size $m\in \{10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 200, 300, 400, 500\}$ and for 6 different architectures (4 x 40 and 10 x 100 with ReLU, ELU, and tanh activations) over two dimensions ($d=4$ and $d=8$) and two coefficients ($a_{1,d}$ from \eqref{eq:affine_param_p2} and $a_{2,d}$ from \eqref{eq:affine_param_p5}), giving 336 DNNs to be trained for each trial. For the Poisson PDE and Navier-Stokes-Brinkman PDE we run 12 trials. For the Boussinesq PDE we run 8 trials due to the larger problem size. For the Poisson PDE, we allocate 336 nodes with $1\times32$ core CPUs running 4 threads (totalling 1344 threads, 6.8 GB RAM per node) each running 3 trials, taking approximately 4 hours and 15 minutes to complete. For the Navier-Stokes-Brinkman PDE, we use the same setup allocating 9.88 GB RAM per node and the runs take approximately 9 hours and 13 minutes for each of the two components of the solution to complete. For the Boussinesq PDE, we allocate 336 nodes with $1\times32$ core CPUs running 4 threads per node (totaling 1344 threads, 10 GB RAM per node) each running 2 trials, taking approximately 12 hours and 32 minutes for each of the 3 components to complete. Given this, the total time required to reproduce the results in parallel with the above setup is approximately 60 hours or 2.5 days. Results were stored locally on the cluster and the estimated total space used to store the data for testing and training and results from computation is approximately 50 GB. Trained models were not retained due to space limitations on the cluster.

\item[(iii)] \textbf{Choice of architectures and initialization}. Based on the strategies in [Adcock and Dexter, 2021], we fix the number of nodes per layer $N$ and depth $L$ such that the ratio $\beta := L/N$ is $\beta=0.5$. In addition, we initialize the weights and biases using the {\tt HeUniform} initializer from {\tt keras} setting the seed to the trial number.
We consider the Rectified Linear Unit (ReLU)
\begin{equation*}
\sigma_1(z):= \max \{ 0,z\},
\end{equation*}
hyperbolic tangent (tanh)
\begin{equation*}
\sigma_2(z):= \frac{\E^z-\E^{-z}}{\E^{z}+\E^{-z}},
\end{equation*}
or Exponential Linear Unit (ELU)
\begin{equation*}
\sigma_3(z)=\begin{cases} z & z>0, \\
\E^{z}-1 & z \leq 0
\end{cases}
\end{equation*}
activation functions in our experiments.

\item[(iv)] \textbf{Optimizers for training and parametrization}. To train the DNNs, we use the \texttt{Adam} optimizer  [Kingma and Ba, 2017], incorporating an exponentially-decaying learning rate. We train our models for 60,000 epochs or until converging to a tolerance level of $\epsilon_{\text{tol}}=5\cdot10^{-7}$ in single precision. In light of the nonmonotonic convergence behavior observed during the minimization of the nonconvex loss (see, e.g., [Adcock et~al., 2021, Adcock and Dexter, 2021]), 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.

\item[(v)] \textbf{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 \S\ref{ss:num-exp-form}, we run several trials solving the problem:
\begin{align*}
&\text{Given measurements} \quad  \{(X_i,Y_i)\}_{i=1}^{m}\subset (\cX \times \cY)^m, \;\; X_i \sim_{\mathrm{i.i.d.}} \mu, \;\; Y_i = F(X_i) + E_i \in \cY, \\ &\text{approximate } F \in L^2_{\mu}(\cX;\cY).
\end{align*}
We generate the measurements $Y_i$ using mixed variational formulations of the parametric elliptic, Navier-Stokes-Brinkman and Boussinesq equations discretized using {\tt FEniCS} with input data $X_i$. The noise $E_i\in \cY$ encompasses the discretization errors from numerical solution. Further details of the discretization can be found in \S\ref{app:PDE_details}.
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))\}_{i=1}^{m_k}$, where $m_k$ 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.

\item[(vi)] \textbf{Testing data and error metric}. The testing data is generated similarly to the training data, obtaining solutions at different points $X_i \in \cX$ for  $i=1,\ldots,m_{\text{test}}$. However, the testing data $\{(X_i,Y_i = F(X_i) + E_i)\}_{i=1}^{m_{\mathrm{test}}}$ are generated by using a deterministic high-order sparse grid collocation method [Nobile et~al., 2008].  In particular, we use sparse grid quadrature rules to compute approximations to the Bochner norms
\begin{equation*}
\nm{F}_{L^2_{\mu}(\cX;\widetilde{\cY})}  =
\left( \int_{\cX} \nm{F(X)}_{\widetilde{\cY}}^2 \D \mu(X)) \right)^{1/2} \approx \left( \sum_{i=1}^{m_{\mathrm{test}}} \| F(X_i) \|_{\widetilde{\cY}}^2 w_i \right)^{1/2},
\end{equation*}
where $\mu=a_d \sharp \varrho_d$ is the pushforward measure defined in \eqref{ss:num-exp-form} and $w_i$, $i=1,\ldots,m_{\mathrm{test}}$ are the quadrature weights. We use these approximations to compute the relative $L^2_{\mu}(\cX;\widetilde{\cY})$ error
\begin{equation*}
e_{F}^{\mathrm{test}}=\frac{\left(\sum_{i=1}^{m_{\mathrm{test}}} \|F(X_i) -\widehat{F}(X_i) \|_{\widetilde{\cY}}^{2} w_i\right)^{1/2}}{\left(\sum_{i=1}^{m_{\mathrm{test}}} \|F(X_i)  \|_{\widetilde{\cY}}^{2} w_i\right)^{1/2}}.
\end{equation*}
We use a high order isotropic Clenshaw Curtis sparse grid quadrature rule to evaluate $e^{\mathrm{test}}_F$, as described in [Adcock et~al., 2021]. This method  shows superior convergence over Monte Carlo integration to evaluate the global Bochner error.  The sparse grid rule gives $m_{\mathrm{test}}$ points at a level $\ell$ for $d$ dimensions. We rely on the \texttt{TASMANIAN} sparse grid toolkit [Stoyanov, 2015, Stoyanov and Webster, 2016, Stoyanov et~al., 2013] for the generation of the isotropic rule to study the generalization performance of the DNN.

\item[(vii)] \textbf{Visualization}. The graphs in Figs. 1-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 \cite[Sec. A.1]{adcock2022sparse} for further discussion about this choice.

\end{itemize}

\section{Description of the parametric PDEs used in the numerical experiments and their discretization}
\label{app:PDE_details}

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.

\subsection{Parametric coefficients}
\label{app:parametric_coeffs}

We consider two parametric coefficients of the form \eqref{a-def-general} in our numerical experiments.
Our first coefficient is
\begin{equation}\label{eq:affine_param_p2}
a_1(\bm{z},\x) = 2.62 + \sum_{j =1 }^{\infty} x_j \frac{\sin(\pi z_1j)}{j^{3/2}},
\quad \forall \bm{z}\in \Omega,
\end{equation}
where $x_j ∈[-1,1]$, $∀j$.
Our second example involves a log-transformed coefficient, which is a rescaling of the example from [Nobile et~al., 2008] of a diffusion coefficient with one-dimensional (layered) spatial dependence given by
\be{
\label{eq:affine_param_p5}
\begin{split}
& = \exp\left(1+ x_1 \left(\frac{\sqrt{\pi}\beta}{2}\right)^{1/2} + \sum_{j=2}^{\infty} \zeta_j \vartheta_j(\bm{z})  x_j\right), \quad \forall \bm{z} \in \Omega,  \\
\zeta_j
& := (\sqrt{\pi} \beta)^{1/2} \exp\left( \frac{-\left( \lfloor j / 2\rfloor \pi \beta\right)^2}{8} \right),
\\
\vartheta_j(\bm{z})
& := \begin{cases} \sin\left( \lfloor j / 2\rfloor \pi z_1/\beta_p \right) &  \text{ if $j$ is even} \\  \cos\left(\lfloor j/2\rfloor \pi z_1/\beta_p \right) & \text{ if $j$ is odd} \end{cases},
\end{split}
where $x_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 →$, with $j=1$ corresponding to \eqref{eq:affine_param_p2} and $j=2$ corresponding to \eqref{eq:affine_param_p5}. Our input samples are then ${X_i = a_j,d(_i)}_i=1^m⊂$ with $_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 $$ describing the thermal conductivity of the fluid. See \S \ref{app:parametric_Boussinesq} and \ef{eq:affine_param_p3}.

\subsection{Relevant spaces}

We require several function space definitions for the development of the mixed variation formulations of the Poisson, Navier-Stokes-Brinkman and Boussinesq equations.
Let $Ω⊂^n$, $n∈{ 2,3 }$, be the physical domain of a PDE. We write $^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_μ(; )$.
We denote Sobolev spaces as $W^s,p(Ω)$ for $s∈$ 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., \cite[Sec. 1.2]{boffi2013mixed} for further details).
We also define the following closed subspace of $^1(Ω)$:
$$\rH_0^1(\Omega):= \overline{C^\infty_0(\Omega)}^{\|\cdot\|_{1,\Omega}}.$$
Here $C^∞_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 is given by
|v|_{1,\Omega} := \| \nabla v\|_{\rL^2(\Omega)} \qquad \|v \|_{1,\Omega}:= \left\{|v|_{1,\Omega}^2 + \|v \|_{\rL^2(\Omega)}^2 \right\}^{1/2} \qquad \forall v \in \rH^1(\Omega).
For scalar functions $u$ and vector fields $$̌, we use $\nabla u$ and $\div (\bm{v})$ to denote their gradient and divergence, respectively. For tensor fields $\bm{\sigma}$ and $\bm{\tau}$, represented by $(\sigma_{i,j})^{n}_{i,j =1}$ and $(\tau_{i,j})^{n}_{i,j =1}$, respectively, we define $\Div(\bm{\sigma})$ as the divergence operator $\div$ acting along the rows of $\bm{\sigma}$, and we  define the trace and the tensor inner-product as

\begin{equation*}
\mathrm{tr} (\bm{\sigma}) = \sum^{n}_{i =1}\sigma_{i,i},  \text{  and  } \bm{\tau}:\bm{\sigma} = \sum^{n}_{i,j =1} \tau_{i,j} \sigma_{i,j},
\end{equation*}

Furthermore, we introduce the notation $\mathbf{L}^p(\Omega)$ and $\bbL^p(\Omega)$ to represent the vectorial and tensorial counterparts of $\mathrm L^p(\Omega)$, and $\mathbf{H}^1(\Omega)$ and $\bbH^1(\Omega)$ for the vectorial and tensorial counterparts of $\mathrm H^1(\Omega)$, respectively.
Keeping this in mind, we introduce the Banach spaces

\begin{equation}\label{def-Hdp}
\begin{split}
\Hdftq & := \Big\{ \bm{v} \in \Ldos: \quad \div (\bm{v}) \in \mathrm L^{q}(\Omega) \Big\}
\\
\Htdftq &:= \Big\{ \bm{\tau} \in \mathbb L^2(\Omega): \quad \Div(\bm{\tau}) \in \mathbf L^{q}(\Omega) \Big\}
\end{split}
\end{equation}

provided with the natural norms

\begin{align*}
\|\bm{v}\|_{\Hdftq} \,&:=\, \|\bm{v}\|_{\mathrm{\mathbf{L}}^2(\Omega)} \,+\, \|\div (\bm{v})\|_{\mathrm L^{q}(\Omega)},
\\
\| \bm{\tau}\|_{\Htdftq} \,&:=\, \| \bm{\tau}\|_{\mathrm{\mathbb{L}}^2(\Omega)} \,+\, \|\Div( \bm{\tau})\|_{\mathrm{\mathbf{L}}^{q}(\Omega)}\,.
\end{align*}

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$.

Often, under certain conditions, such as incompressibility conditions \cite[eq.(2.4)]{gatica2021banach}, it is convenient to define variations of these spaces. For example, we define
\begin{equation}
\Ltras \,:=\, \Big\{ \bm{\tau} \in  \bbL^2(\Omega): \mathrm{tr}(\bm{\tau}) \,=\, 0\Big\}\,,
\end{equation}
which represents the space of integrable functions with zero trace over $Ω$. Furthermore, given the decomposition (see, e.g.,  [Gatica, 2014])
\begin{equation}\label{decomp-1}
\Htdft \,=\, \Htdoft \,\oplus\, \mathrm R\,\mathbb{I}\,,
\end{equation}
we may also consider
\begin{equation}\label{decomp-2}
\Htdoft \,:=\, \Big\{\bm{\tau} \in \Htdft: \into \mathrm{tr}(\bm{\tau})=0 \Big\}\,,
\end{equation}
as the space of elements in $$ with zero mean trace. Finally, we define
\begin{equation}
\bbL_{\text{skew}}^2(\Omega) = \{ \bm{\eta} \in \bbL^2(\Omega) : \; \bm{\eta} + \bm{\eta}^t = 0\},
\end{equation}
and the space of $^2(Ω)$ functions with zero integral over $Ω$ as
\begin{equation}
\rL_0^2(\Omega) = \left\{ \nu \in \rL^2(\Omega): \; \int_{\Omega} \nu = 0 \right\}.
\end{equation}

\subsection{The parametric diffusion equation}
\label{app:parametric_diffusion}

We now describe our first example, which is the parametric elliptic diffusion equation.
Let $Ω⊂^2$ be a bounded Lipschitz domain, $ ∂Ω$ be the boundary of $Ω$, $f ∈^2(Ω)$ and  $g ∈^1/2(∂Ω)$.
Given $∈[-1,1]^d$ consider the linear elliptic equation with Dirichlet boundary conditions
\begin{align*}
-\div (a(\x) \nabla u (\x)) &= f,   \quad \text{ in } \Omega   \numberthis\label{eq:Poisson} \\
u(\x) &=g, \quad  \text{ on }  \partial \Omega.
\end{align*}
Here $a(x) = a(·, x) ∈^∞(Ω) =: $ is the parametric diffusion coefficient. The terms $f$ and $g$ are not parametric.

\subsubsection{Mixed variational formulation}\label{ss:elliptic-diffusion-variational}

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 $= ^1_0(Ω)$ (see Remark \ref{app:parametric_diffusion}) and instead consider a \textit{mixed} variational formulation of \eqref{eq:Poisson}. 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 [Gatica, 2014] and references within.

Assume that there exists $r,M>0$ such that
\be{
\label{rM-bounds}
0< r \leq   \mathrm{essinf}_{\bm{z} \in \Omega} a(\bm{z},\x) = : a_{\min}(\x) \text{  and  } a_{\max}(\x) : =\mathrm{esssup}_{\bm{z} \in \Omega} a(\bm{z},\x) \leq M, \quad \forall \x \in [-1,1]^d.
Then,  the problem can be stated as a first-order system: given $∈[-1,1]^d$,  find   $(σ,u)() ∈×̋$ such that
\be{
\label{eq:mixed_poisson}
\begin{split}
d_{a(\x)}(\bm{\sigma}(\x),\bm{\tau})+b(\bm{\tau},u(\x)) &=G(\bm{\tau}), \quad \forall \bm{\tau} \in \H, \\
b(\bm{\sigma},v) &=F(v), \quad  \forall v \in \Ltwo.
\end{split}
Here $d$ and $b$ are the bilinear forms defined by
\begin{align*}
d_{a(\x)}(\bm{\sigma},\bm{\tau})= \into \frac{\bm{\sigma} \cdot \bm{\tau}}{a(\x)}, \quad \forall (\bm{\tau},\bm{\sigma}) \in \H\times \H \\
b(\bm{\tau},v) = \into \div(\bm{\tau})  v , \quad \forall (\bm{\tau}, v) \in \H\times \rL^2(\Omega)
\end{align*}
and the functionals $G ∈()̋'$ and $J ∈()'$ are defined by
\begin{equation}
J(v)= -\into f v, \quad \forall v \in \Ltwo, \text{ and } G(\bm{\tau})= \ip{\gamma(\bm{\tau}) \cdot \bm{n}}{g}_{1/2,\partial \Omega}, \quad \forall \bm{\tau} \in \H.
\end{equation}
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,) = 0.5$ on the bottom of the boundary $(0,1)×{0}$, and zero boundary conditions on the remainder of the boundary.

\subsubsection{Holomorphy assumption}\label{ss:elliptic-holomorphy}

Consider the affine parametrization \eqref{eq:affine_param_p2}. Setting $M = 2.7$ and observing that
\begin{equation}
\left|\sum_{j \in \bbN} x_j \frac{\sin(\pi z_1j)}{j^{3/5}} \right| \leq \sum_{j \in \bbN}   \frac{1}{j^{3/5}}   \approx 2.61238= 2.62-r,
\end{equation}
for some $r<0.00762$, we deduce that \ef{rM-bounds} hold, which makes the mixed variational formulation \eqref{eq:mixed_poisson} well defined, i.e., for each $∈[-1,1]^∞$, there exists a unique solution $(σ,u)() ∈×̋$. Moreover, one can show that the parametric solution map $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_^∞(Ω) = j^-3/2$. See \cite[Prop.\ A.3.2]{moraga2024optimal}.

In view of this property, this example falls within our theory. Note that $b ∈ℓ^p_𝖬()$ for every $p<2/3$. Thus, we expect a theoretical rate of convergence with respect to the amount of samples that is arbitrarily close to $m^1/2-3/2=m^-1$. This holomorphy result applies to the affine diffusion \ef{eq:affine_param_p2}, not the log-transformed diffusion \ef{eq:affine_param_p5}. However, we expect that it is possible to extend \cite[Prop.\ A.3.2]{moraga2024optimal} to the latter case.

\subsubsection{Finite element discretization}

We use so-called \textit{conforming} Finite Element (FE) discretizations \cite[Chp. 3]{ciarlet2002finite}. Given a number of degrees of freedom $K$, this results in finite-dimensional spaces $H_K ⊆$̋ and $Q_K \subseteq \Ltwo$. Specifically, we consider a regular triangulation $\cT_K$  be of $\overline{\Omega}$ made up of triangles $T$ of minimum diameter $h_{\text{min}}=0.0844$ and maximum diameter $h_{\text{max}}=0.1146$. This corresponds to a total number of degrees of freedom $K=2622$. See Fig. <ref> for an illustration of the FE mesh.

The domain $\Omega$ and FE mesh for the parametric diffusion equation.

Fig.  <ref> shows a comparison between a reference solution computed by the FEniCS FEM solver and the approximation obtained by an ELU $4\times 40$ DNN.

The solution $\u(\x)$ of the parametric Poisson problem in (<ref>) for a given parameter $\x=(1,0,0,0)^{\top}$ with affine coefficient $a_{1,d}$ and $d=4$, using a total of $K=2,622$ degrees of freedom (DoF). The left plot shows the solution given by the FEM solver. The right plot show the ELU $4\times 40$ DNN approximation after $60,000$ epochs of training with $m=500$ sample points for training.

§.§ Parametric Navier-Stokes-Brinkman equations

We next consider a parametric model describing the dynamics of a viscous fluid through porous media. Consider a bounded and Lipschitz physical domain $\Omega \subseteq \bbR^2$. Given $\x \in [-1,1]^d$, we consider the incompressible nonlinear stationary Navier-Stokes-Brinkman (NSB) equations: find $\u: [-1,1]^d \times \Omega \rightarrow \bbR^2$ and $p: [-1,1]^d \times \Omega \rightarrow \bbR$ such that

\begin{align*}
\eta \u - \lambda \Div( a(\x) \bm{e}(\u(\x))) + (\u(\x) \cdot \nabla) \u(\x) + \nabla p(\x) &=f, \quad \text{ in } \Omega  \\
\div (\u(\x)) &= 0, \quad \text{ in } \Omega \numberthis \label{eq:NSB}  \\
\u
&= \begin{cases}\u_{D} , \quad &\text{ on } \partial \Omega_{\mathrm{in}} \\
0, \quad &\text{ on } \partial  \Omega_{\mathrm{wall}}    \\
\end{cases}   \\
(a(\x) \nabla  \bm{e}(\u) -p \bbI)\nu &=0, \quad \text{ on } \partial \Omega_{\mathrm{out}}   \\
\int_{\Omega} p &= 0,
\end{align*}

where $\lambda= \text{Re}^{-1}$ and $\text{Re}$ is the Reynolds number, $a(\x) = a(\cdot;\x)\in \cX := \rL^\infty(\Omega)$ with $a: [-1,1]^d \times \Omega \rightarrow \bbR^+$ the random viscosity of the fluid, $\eta \in  \bbR^+$ is the scaled inverse permeability of the porous media, $\u$ is the velocity of the fluid,  $e()̆ = 1/2 (∇+̆ (∇)̆^t)$ is the symmetric part of the gradient and $p$ is the pressure of the fluid, and $f: Ω→$ is an external force  independent of the parameters. Here, the fourth condition imposes a zero normal Cauchy stress
\begin{equation*}
(a(\x) \nabla  \bm{e}(\u) -p \bbI)\nu =0
\end{equation*}
for the output boundary on $∂Ω_out$. Note as well that the incompressibility of the fluid imposes the following compatibility condition on $_̆D$:
\begin{equation*}
\int_{\Gamma} \u_D \cdot \bm{n}=0,\quad \text{ on }\partial \Omega_{\mathrm{in}}.
\end{equation*}
The third condition also imposes a no-slip condition on the walls $Ω_wall$ \cite[eq.(2.3)]{gatica2022non-augmented}.

\subsubsection{Mixed variational formulation}\label{ss:NSB-variational}

The analysis of the detailed mixed formulation used for this problem in the nonparametric case can be found in [Gatica et~al., 2022]. 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 \cite[Sec. 1]{gatica2022non-augmented}.

Based on the analysis in [Gatica et~al., 2022], the mixed variational formulation of the parametric NSB equations in \eqref{eq:NSB} becomes: given  $∈[-1,1]^d$, find $(,̆t,σ,γ)() ∈×^2_tr(Ω) ××_skew^2(Ω)$ such that
\begin{equation} \label{eq:NSB_weak1}
\begin{array}{rcll}
\disp
\lambda \into a_i(\x) \bm{t}(\x) : \bm{s} \,-\,   \into \bm{s}:\bm{\sigma}(\x) \,-\,
\disp \into(\u \otimes \u)(\x):\bm{s} &=& 0, & \\[2ex]
\disp
\into \bm{t}(\x): \bm{\tau} \,+\,   \into \bm{\gamma}(\x):\bm{\tau} \,+\,
\disp \into \u(\x) \cdot \Div(\bm{\tau}) &=& \ip{\bm{\tau} \bm{n}}{\u_D}_{\partial \Omega_{\text{in}}}, & \\[2ex]
\disp
\into   \bm{\delta} : \bm{\sigma}(\x) \,+\,   \into \bm{v}\cdot \Div (\bm{\sigma}(\x)) \,-\,
\disp \into \eta \u(\x) \cdot \v &=& \disp \into \bm{f} \cdot \v, & \\[2ex]
\end{array}
\end{equation}
for all $(v,s,τ,δ) ∈×××_skew^2(Ω)$. Numerically, the skew-symmetry of  $γ$ is imposed by searching for $γ∈L^2(Ω)$ and setting
0 & \gamma \\
-\gamma & 0 \\
\end{bmatrix}.
Moreover, we  impose the Neumann boundary condition via a Nietsche method as in \cite[Sec. 5.2]{gatica2022non-augmented}. Specifically, we add
\begin{equation*}
\kappa\ip{(\bm{\sigma}+ \u \otimes \u)\bm{n})}{\bm{\tau} \bm{n}}_{\partial \Omega_{\text{out}}}=0
\end{equation*}
to the second equation where $κ≫1$ is a large constant (e.g., $κ= 10^4$). As usual in this formulation, the pressure $p ∈^2(Ω)$ can be computed according to the post-processing formula
\begin{equation*}
p = - \dfrac{1}{2} \mathrm{tr} (\bm{\sigma}+(\u \otimes \u)).
\end{equation*}
Note that above  we omitted the term $$ for simplicity.

In our experiments, we consider to approximate solutions to the parametric NSB problem with $λ=0.1$,  a scaled inverse permeability of $η= 10+z_1^2+z_2^2$, an external force $f=(0,-1)^⊤$  and  random viscosity $a_j,d$ as in \eqref{eq:affine_param_p2}--\eqref{eq:affine_param_p5} with $j∈{1,2}$.

We consider the unit square $Ω=(0,1)^2$ as the domain. We consider an inlet boundary defined by $∂Ω_in= (0,1) ×{1}$, an outlet boundary  $∂Ω_out= {1} ×(0,1)$ and walls  defined by $∂Ω_wall  = {0} ×(0,1) ∪(0,1) ×{0}$. For simplicity, we use the same mesh as that of the previous example. See Fig. \ref{Example_affine_fig}.
On the Neumann boundary $∂Ω_out$ we consider a zero normal Cauchy stress, a Dirichlet condition $_̆D=(0.0625)^-1 ((z_2-0.5)(1-z_2),0)$ on  $∂Ω_in$ and a no-slip velocity on $∂Ω_wall$.

Fig.\  \ref{Example_NSB_fig} provides a comparison between a reference solution of the vector field $u$ and pressure $p$ computed by the {\tt FEniCS} FEM solver and the approximation generated by a ELU $4×40$ DNN.

\begin{remark}[Other auxiliary variables]\label{rmk:othervar}
We report the performance of the DNNs approximating $(\bm{u},p)(\x) \in \Lft  \times \Ltwo$. Note that any solver based on the above formulation outputs several other variables, e.g., $(\bm{t},\bm{\sigma},\bm{\gamma})(\x) \in   \bbL^2_{\text{tr}}(\Omega) \times \bbH_0(\Div_{4/3}; \Omega) \times \bbL_{\text{skew}}^2(\Omega)$. 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.
\end{remark}

% ===============================================================================================
\begin{figure}[t]
\centering
\includegraphics[scale=0.15]{NSBuh.pdf}
\includegraphics[scale=0.15]{NSBph.pdf}\\
\includegraphics[scale=0.15]{NSBuhDNN.pdf}
\includegraphics[scale=0.15]{NSBphDNN.pdf}
\caption{The solution $(\u, p)(\x)$ of the parametric NSB problem \eqref{eq:NSB} for a given parameter $\x=(1,0,0,0)^{\top}$ with affine coefficient $a_{1,d}$ and $d=4$, using a total of $1,464$ degrees of freedom (DoF) for $\u$ and $244$ DoF for $p$. The top row shows the solution given by the FEM solver, and the bottom row shows the ELU $4\times 40$ DNN approximation after $60,000$ epochs of training with $m=500$ sample points. The left plots show the vector field $\u$.  The right plots show the points of highest pressure $p$.
\label{Example_NSB_fig}
\end{figure}

To conclude this discussion, in Fig.\ \ref{NSB_res-additional} we plot the numerical results for the approximation of the pressure $p$ in the above problem. This complements Fig.\ \ref{NSB_res}, 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 $m^-1$ and there is no degradation with increasing dimension $d$.

\begin{figure}[t]
\centering
\includegraphics[scale=0.125]{NSB_S3_p_d4_.pdf}
\includegraphics[scale=0.125]{NSB_S3_p_d8_.pdf}
\includegraphics[scale=0.125]{NSB_logKL_p_d4_.pdf}
\includegraphics[scale=0.125]{NSB_logKL_p_d8_.pdf}
\caption{The same as Fig.\ \ref{NSB_res}, except showing results for the pressure $p$.}
\label{NSB_res-additional}
\end{figure}

\subsection{Parametric stationary Boussinesq equation}
\label{app:parametric_Boussinesq}

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 [Cohen and DeVore, 2015, Dexter et~al., 2019, Adcock et~al., 2021] study the more restrictive case of homogeneous Dirichlet boundary conditions, where $u ∈_0^1(Ω)$. 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 is aligns with what we expect. We conjecture that this 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 ($Ω⊂^3$) with two random coefficients affecting different parts of the coupled problem. The nonparametric version of this problem is based on [Colmenares et~al., 2020].

Specifically, we consider the Boussinesq formulation in [Colmenares et~al., 2020] 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 $∈[-1,1]^d$, our goal is to find the velocity $:̆ [-1,1]^d ×Ω→^2$,  pressure $p: [-1,1]^d ×Ω→$ and temperature $φ: [-1,1]^d ×Ω→$ of a fluid such that
\begin{align*}
-  \Div( 2 a(\x) \varpi(\varphi(\x)) \bm{e}(\u(\x))) + (\u(\x) \cdot \nabla) \u(\x) + \nabla p(\x) &=\varphi(\x) \bm{g}, \quad \text{ in } \Omega,  \\
\div (\u(\x)) &= 0, \quad \text{ in } \Omega, \numberthis \label{eq:BSQ}  \\
-  \div(  \mathbb{K}(\x)  \nabla \varphi(\x)) + \u(\x)\cdot \nabla \varphi(\x) &=0, \quad \text{ in } \Omega,  \\
\u   &=  \u_{D} , \text{ on } \partial \Omega,    \\
\varphi   &=  \varphi_{D} , \text{ on } \partial \Omega,\\
\int_{\Omega} p(\x) &= 0.
\end{align*}
Here $g=(0,0,-1)^⊤$ is a gravitational force, $() = (·; )∈^∞(Ω)$  with $: [-1,1]^d ×Ω→^3 ×3$ is a parametric uniformly-positive tensor describing the thermal conductivity of the fluid, and is given by
\begin{equation}\label{eq:affine_param_p3}
\bbK(\bm{z},\x) = \left(1.89 + \sum_{j \in \bbN }  x_j \frac{\sin(\pi z_3 j)}{j^{9/5}} \right) \begin{bmatrix}
\exp(-z_1) & 0 & 0 \\
0 & \exp(-z_2) & 0 \\
0 &  0 & \exp(-z_3)
\end{bmatrix},
\ \forall \bm{z}\in \Omega,
\end{equation}
for $x_j ∈[-1,1]$, $∀j$, $ϖ: →^+$ is the temperature-dependent viscosity given by $ϖ(φ)=0.1+exp(-φ)$ and $a() = a(·; ) ∈^∞(Ω)$with $a: [-1,1]^d ×Ω→$ is a parametric variable affecting the viscosity of the fluid.
As in the previous example $e()̆$ is the symmetric part of $∇$̆. Note that in this case, we have $(a(\bm{x}) , \bbK(\bm{x}) ) \in \cX := \rL^{\infty}(\Omega) \times \bbL^{\infty}(\Omega)$.

§.§.§ Fully mixed variational formulation

The complete derivation of a fully-mixed variational formulation for the non-parametric Boussinesq equation in Banach spaces can be found in <cit.>. To make the presentation simpler  we rewrite it for the parametric case. Given $\x \in [-1,1]^d$, find $(\bm{u},\bm{t},\bm{\sigma},\varphi,\tilde{\bm{t}},\tilde{\bm{\sigma}})(\x)\in \Lft \times \Ltras \times \Htdoft \times \Lf \times \Ldos \times \Hdft$ such that

\begin{equation}\label{weak-Boussinesq}
\begin{array}{rcll}
\disp
-\into \v \cdot \Div (\bm{\sigma}(\x)) \,+\, \dfrac{1}{2} \into \bm{t}(\x) \u(\x) \cdot \v \,-\,
\disp \into \varphi(\x) \bm{g} \cdot \v &=& 0 &\quad \forall\, \v \in \Lft\,,\\[2ex]
\disp \into   2a(\x)\varpi(\varphi(\x)) \bm{t}_{\mathrm{sym}}(\x): \bm{s} \,-\,  \dfrac{1}{2} \disp \into (\u \otimes \u)(\x):\bm{s} &=&
\disp \into \bm{\sigma}(\x) : \bm{s} &\qquad \forall\, \bm{s} \in \Ltras \,,\\[2ex]
\disp
\into \bm{\tau} :\bm{t}(\x) \,+\, \into \u(\x)\cdot \Div (\bm{\tau}) &=&  \ip{\bm{\tau} \bnu}{\u_D}_{\Gamma} &\qquad
\forall \, \bm{\tau} \in \Htdoft \,,\\[2ex]
\disp
-\into \psi \, \div (\tilde{\bm{\sigma}}(\x)) \,+\, \dfrac{1}{2} \into  \psi(\x)   \u(\x) \cdot \tilde{\bm{t}} &=&
0  &\qquad \forall \,\psi \in \Lf \,,\\[2ex]
\disp
\into \mathbb{K}(\x) \tilde{\bm{t}}(\x) \cdot \tilde{\bm{s}}-\dfrac{1}{2} \into \varphi(\x) \u(\x) \cdot \tilde{\bm{s}} &=&  \disp \into \tilde{\bm{\sigma}}(\x) \cdot \tilde{\bm{s}}
&\qquad \forall\, \tilde{\bm{\sigma}}  \in \Ldos\,, \\[2ex]
\disp
\into \tilde{\bm{\tau}}  \cdot \tilde{\bm{t}}(\x)  \,+\, \into \varphi(\x) \, \div (\tilde{\bm{\tau}} ) &=&
\ip{\tilde{\bm{\tau}}  \cdot \bnu}{\varphi_D}_{\Gamma} &\qquad \forall\, \tilde{\bm{\tau}}  \in \Hdft \,, \\[2ex]
\disp
\disp \int_{\Omega} \mathrm{tr}(2\bm{\sigma}+ \u \otimes \u )(\x) &=&
0.  &\\[2ex]
\end{array}
\end{equation}

Here $p \in \rL_0^2(\Omega)$ can be recovered by using

\begin{equation}
p= -\dfrac{1}{6}\mathrm{tr} (2\sigma+2c \mathbb{I}+ \u \otimes \u),  \text{ with } c = -\dfrac{1}{6|\Omega|} \into \mathrm{tr}(\u \otimes \u).
\end{equation}

As in the previous example, we omitted the term $\x$ for simplicity from this equation.
For further details on this formulation we refer to [Colmenares et al., 2020] and references within.

Given $\x \in [-1,1]^d$, we approximate the solution $(\u,p,\varphi)(\x) \in  (\Lft \times \rL_0^2(\Omega) \times \Lf)$ of (<ref>) by using DNNs and study the approximation capabilities as we increase the number of training samples $m$. As in the previous example (see Remark <ref>), we do not aim to approximate the other variables $(\bm{t},\bm{\sigma},\tilde{\bm{t}},\tilde{\bm{\sigma}})(\x)\in \Ltras \times \Htdoft \times \Ldos \times \Hdft$.

In our experiments, we consider the unit cube $\Omega=(0,1)^3$ as the domain in $\bbR^3$. We consider a nonzero boundary condition $u_D=(1,1,0)$ on the bottom face of the cube $\partial \Omega_{\mathrm{bottom}} = (0,1) \times (0,1) \times \{ 0\}$, and zero on the rest of the faces. We set $\varphi_D=\exp(4(-(z_1-0.5)^2-(z_2-0.5)^2))$ on $\partial \Omega_{\mathrm{bottom}}$ and zero otherwise. For simplicity, we consider the same parametric coefficients $a_{1,d}$ and $a_{2,d}$ given by (<ref>) and (<ref>), respectively. See Fig. <ref> for an example of the solution $(\u,p,\varphi)(\x)$ for a given $\x \in [-1,1]^4$.

The solution $(\u, \varphi, p)(\x)$ to the parametric Boussinesq problem in (<ref>) for a given parameter $\x=(1,0,0,0)^{\top}$ with affine coefficient $a_{1,d}$ and $d=4$, using a total of $18,480$ degrees of freedom (DoF) for $\u$ and $528$ DoF for both $φ$ and $p$. The top row shows the solution given by the FEM solver, and on the bottom is the $4×40$ ELU--DNN approximation after $60,000$ epochs of training with $m=500$ sample points. The left plots show streamlines of the vector field $$̆ and their direction indicated with coloured arrows. The middle plots visualize the temperature distribution inside the cube using coloured spheres, with the hottest region at the centre of the cube. The right plots illustrates the points of highest pressure $p$.

To conclude this section, we provide a comparison of the performance of the DNN architectures in approximating the velocity field $\bm{u}$ and pressure $p$ for the Boussinesq PDE. This complements  Fig. <ref>, which showed results for the temperature $\varphi$. As with $\varphi$, the convergence rate for $p$ agrees roughly with the rate $m^{-1}$, and does not appear to deteriorate with the dimension $d$. On the other hand, the convergence rate for the velocity field $\bm{u}$ is somewhat slower.

The same as Fig. <ref>, except showing results for the velocity field $\bm{u}$, where $\cY = \rL^4(\Omega)$, and pressure $p$, where $\cY = \rL^2_0(\Omega)$.

§ OVERVIEW OF THE PROOFS

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.

§.§ Additional notation

§.§.§ Lipschitz constants

Let $(\cX,\nm{\cdot}_{\cX})$ and $(\cY,\nm{\cdot}_{\cY})$ be Banach spaces, $G : \cX \rightarrow \cY$ and $\cB \subseteq \cX$. We define the Lipschitz constant of $G$ as $L = \mathrm{Lip}(G ; \cB , \cY)$ as the smallest constant $L \geq 0$ such that

G(X') - G(X)_ ≤L X'-X_,  ∀X,X' ∈.

§.§.§ Sequence spaces

We require some notation for sequences. Let $(\cZ , \nm{\cdot}_{\cZ})$ be a Banach space, $d = \bbN \cup \{ \infty \}$ and write $\bm{\nu} = (\nu_k)^{d}_{k=1}$ for an arbitrary multi-index in $\bbN^d_0$. If $\Lambda \subseteq \bbN^d_0$ is a finite or countable set of multi-indices and $0 < p \leq \infty$ we define the space $\ell^p(\Lambda ; \cZ)$ as the set of all $\cZ$-valued sequences $\bm{c} = (c_{\bm{\nu}})_{\bm{\nu} \in \Lambda}$ for which $\nm{\bm{c}}_{p;\cZ} < \infty$, where

c_p; =  ( ∑_ν ∈Λ c_ν^p_ )^1/p     0 < p < ∞,

sup_ν ∈Λ c_ν_     p = ∞.

When $(\cZ,\nm{\cdot}_{\cZ}) = (\bbR,\abs{\cdot})$, we just write $\ell^p(\Lambda)$ and $\nm{\cdot}_{p}$. We also define the following:

c_p; = sup_z^* ∈B(^*) z^*(c)_p,  ∀c ∈ℓ^p(Λ; ),

where $z^*(\bm{c}) = (z^*(c_{\bm{\nu}}))_{\bm{\nu} \in \Lambda} \in \bbR^{|\Lambda|}$ for $\bm{c} = (c_{\bm{\nu}})_{\bm{\nu} \in \Lambda}$. Notice that $\tnm{\bm{c}}_{p;\cZ}  \leq \nm{\bm{c}}_{p;\cZ}$.

Given a sequence $\bm{c} = (c_{\bm{\nu}})_{\bm{\nu} \in \Lambda}$, we write

supp(c) = { ν ∈Λ: c_ν ≠0 } ⊆Λ.

For $d \in \bbN \cup \{ \infty \}$, we write $\bm{e}_j$, $j = 1,\ldots,d$, for the canonical basis vectors in $\bbR^d$. We also write $\bm{0}$ for the multi-index of zeros and $\bm{1}$ for the multi-index of ones.
Finally, for multi-indices $\bm{\nu} = (\nu_k)^{d}_{k=1}$ and $\bm{\mu} = (\mu_k)^{d}_{k=1}$, we write $\bm{\nu} \geq \bm{\mu}$ to mean $\nu_k \geq \mu_k$, $\forall k$ and likewise for $\bm{\nu} > \bm{\mu}$.

§.§.§ Weights and weighted sequence spaces

Let $d \in \bbN \cup \{ \infty \}$, $\Lambda \subseteq \bbN^d_0$ and $\bm{w} = (w_{\bm{\nu}})_{\bm{\nu} \in \Lambda} \geq \bm{0}$ be a sequence of nonnegative weights. We define the weighted cardinality of a set $S \subseteq \Lambda$ as

| S |_w  = ∑_ν ∈S w^2_ν .

For a Banach space $(\cZ , \nm{\cdot}_{\cZ})$ and $0 < p \leq 2$, we define the weighted $\ell^p_{\bm{w}}({\Lambda}; \cZ)$ space as the space of $\cZ$-valued sequences $\bm{c} = (c_{\bm{\nu}})_{\bm{\nu} \in {\Lambda}}$ for which $\nm{\bm{c}}_{p,\bm{w};\cZ} < \infty$, where

c_p,u; =  ( ∑_ν ∈Λ w^2-p_ν c_ν^p_ )^1/p .

Notice that $\| \cdot \|_{2,\bm{u} ; \cZ}$ coincides with the unweighted norm $\| \cdot \|_{2;\cZ}$.

Finally, given a sequence $\bm{c} = (c_{\bm{\nu}})_{\bm{\nu} \in {\Lambda}}$ we write

c_0,w = | supp(c) |_w.

§.§.§ Legendre polynomials

We write $\{ P_n \}_{n \in \bbN_0}$ for the classical Legendre polynomials on $[-1,1]$ with normalization $P_n(1) = 1$. Since $\int^{1}_{-1} | P_n(x)|^2 \D x = (n+1/2)^{-1}$, we define the orthonormal (with respect to the uniform probability measure) Legendre polynomials as

ψ_n(x) = √(2n+1) P_n(x),  x ∈[-1,1], n ∈_0.

Let $D = [-1,1]^{\bbN}$ as before,

= { ν = (ν_i )^∞_i=1 ∈^_0 : | supp(ν) | < ∞}

be the set of multi-indices with finitely many nonzero terms and define the multivariate Legendre polynomials as

Ψ_ν(x) = ∏_i ∈ ψ_ν_i(y_i) ≡∏_i ∈(ν) ψ_ν_i(x_i),   ∀x ∈D, ν = (ν_i)^∞_i=1 ∈.

Here, the second equality follows from the fact that $\psi_0 = 1$. Then it is known that the set

{ Ψ_ν }_ν ∈ ⊂L^2_ϱ(D)

constitutes an orthonormal basis for $L^2_{\varrho}(D)$ <cit.>. For later use, we also define the sequences

u = (u_ν)_ν ∈,  where u_ν = Ψ_ν_L^∞_ϱ(D) = ∏_k ∈ √(2 ν_k + 1),  ∀ν ∈,

v = (v_ν)_ν ∈,  where v_ν = u^5+ξ_ν, ∀ν ∈.

Here $\xi \geq 0$ will be chosen later in the proof.

§.§.§ Miscellaneous

Given an optimization problem $\min_t f(t)$, we say that $\hat{t}$ is a $(\sigma,\tau)$-approximate minimizer for some $\sigma \geq 1$ and $\tau \geq 0$ if $f(\hat{t}) \leq \sigma^2 \min_t f(t) + \tau^2$. (In the main paper, we consider only $\sigma = 1$, but for the proofs it is useful to allow $\sigma > 1$).

Finally, for convenience, given $X_1,\ldots,X_m \in \cX$, we now define the semi-norm

G_𝖽𝗂𝗌𝖼,μ = √(1/m ∑^m_i=1 G(X_i)^2_)

of an operator $G : \cX \rightarrow \cY$ and

g_𝖽𝗂𝗌𝖼,ς̃ = g ∘__𝖽𝗂𝗌𝖼,μ =√(1/m ∑^m_i=1 g∘_ (X_i)^2_)

of a function $g : \bbR^{\bbN} \rightarrow \cY$.

§.§ Overview of the proofs

§.§.§ Theorem <ref>

Theorems <ref>-<ref> are based on polynomials, and specifically, procedures for learning efficient Legendre polynomial approximations to holomorphic operators. As observed, these results are based on recent work on learning holomorphic, Banach-valued functions  [Adcock et al., 2021, Adcock et al., 2022, Adcock et al., 2024]. See also Remark <ref> below.

The proof of Theorem <ref> involves three mains steps.

(a) Formulation (<ref>) and analysis (<ref>-<ref>) of a suitable polynomial learning procedure.

(b) Construction of a family of DNNs that approximately emulates the Legendre polynomials (<ref>).

(c) Analysis of the corresponding training problem DNN-training-prob (<ref>–<ref>).

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 poly-set-def 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 poly-min-prob relies heavily on techniques from compressed sensing. In order to bound the encoding-decoding error, we also require several results on Lipschitz continuity (Lemma <ref>) and norm equivalences (Lemma <ref>) for multivariate polynomials.

Step (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 <ref>, then use this to define the DNN family $\cN$ in DNN-arch-def.

We then analyze the DNN training problem DNN-training-prob in step (c). Using this emulation result, we first show that any approximate minimize $\widehat{N}$ of DNN-training-prob yields a polynomial that is also an approximate minimizer of the polynomial training problem poly-min-prob (Lemma <ref>). We may then apply the results shown in Step (a) to prove a generalization bound (Theorem <ref>). 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 (<ref>). Finally, in <ref> we put all these estimates together to complete the proof.

Theorem <ref> is a generalization and improvement of <cit.>, which deals with the case of learning Banach-valued functions rather than operators. Specifically, the setting of [Adcock et al., 2022] can be considered a special case of this paper where $\cX = \ell^{\infty}(\bbN)$, in which case the encoding error $E_{\cX,q} = 0$ as well. Moreover, Theorem <ref> improves the main results of [Adcock et al., 2022] in three key ways. First, the the DNN architectures bounds are much narrower: $\mathrm{width}(\cN) \lesssim (m/L)^{1+\delta}$ versus $\mathrm{width}(\cN) \lesssim m^{3 + \log_2(m)}$ in the latter. Second, Theorem <ref> considers standard training, i.e., $\ell^2$-loss minimization, whereas <cit.> requires regularization. Third, for Banach spaces, the error decay rate with respect to $m$ is roughly doubled: $E_{\mathsf{app},2} = \ord{m^{1-1/p}}$ in Theorem <ref> versus $\ord{m^{1/2(1-1/p)}}$ in <cit.>. Finally, Theorem <ref> also provides bounds in the $L^{\infty}_{\mu}$-norm, whereas the results in [Adcock et al., 2022] only consider the $L^2_{\mu}$-norm.

§.§.§ Theorem <ref>

The proof of Theorem <ref> relies of three key steps.

(a) Using a minimizer of the polynomial training problem poly-min-prob to construct a family of minimizers of the DNN training problem DNN-training-prob.

(b) Analysis of the corresponding minimizers using the previously-derived bound for polynomial minimizers.

(c) Using the Lipschitz continuity of DNNs to show stability of the DNN minimizer and a permutation argument to show the existence of many parameters that lead to equally `good' minimizers.

Given any minimizer of the polynomial training problem poly-min-prob, it is straightforward to define a DNN with the desired generalization bound. Unfortunately, this will generally not be a minimizer of DNN-training-prob. To achieve the aims of step (a) we proceed as follows. First, we note that a DNN will be a minimizer if the corresponding approximation satisfies $\widehat{F}(X_i) = \widetilde{Y}_i$, where $\widetilde{Y}_i$ are the closest points to the $Y_i$ from $\widetilde{\cY} = \cD_{\cY}(\bbR^{d_{\cY}})$. To achieve this, we take the existing DNN then add on a suitable number of additional terms corresponding to the first $r > m$ order-one Legendre polynomials (<ref>). We show that by doing this, we can construct a DNN for which $\widehat{F}(X_i) = \widetilde{Y}_i$ (Lemma <ref>).

In Step (b), we first bound this DNN minimizer in terms of the polynomial minimizer plus the contributions of these additional terms (Lemma <ref>). The latter involves the minimal singular value of a certain $m \times r$ matrix $\bm{B}$, which is the matrix of the linear system that enforces the condition $\widehat{F}(X_i) = \widetilde{Y}_i$. We bound this minimal singular value in <ref>.

We then use this to complete the proof of part (A) of Theorem <ref> in <ref>. In this section, we also complete step (c) to establish parts (B) and (C).

Like Theorem <ref>, Theorem <ref> also relies on ideas from [Adcock et al., 2022]. However, [Adcock et al., 2022] does not address fully-connected DNN architectures. To address this challenge, the proof of Theorem <ref> involves the technical construction described above.

§.§.§ Theorem <ref>

Theorems <ref> is based on <cit.>. The basic idea is to consider a family of affine, holomorphic operators f-def-lower-bd. This allows us to lower bound the quantity $\theta_m(\bm{b})$ by the so-called Gelfand width def:mwidths of a certain weighted unit ball in a finite-dimensional space. Bounds for such Gelfand widths are known, and this allows us to derive the corresponding result.

§.§.§ Theorem <ref>

Theorems <ref> employs similar ideas, but in a more technical manner. We consider a family of linear, holomorphic operators fc-def, which involves a sum over $r$ groups of $m+1$ coefficients. We restrict to coefficients that lie in the null space of the corresponding sampling operator. Then, through a series of inequalities, we can lower bound $\tilde\theta_m(\bm{b})$ by a sum over $r$ terms, each involving the $\ell^1$-norms of certain vectors in the null space of the matrix of the corresponding sampling operator exp-err-split. This matrix is a subgaussian random matrix. We now use a technical estimate from [Nguyen and Vu, 2018] for vectors in the null space of subgaussian random matrices, which shows that they cannot be `spiky'. Applying this and a series of further inequalities yields the result.

§ PROOF OF THEOREM <REF>

§.§ Formulation of an approximate polynomial training problem

Let $\Lambda \subset \cF$ with $\supp(\bm{\nu}) \subseteq \{1,\ldots,d_{\cX} \}$, $\forall \bm{\nu} \in \Lambda$, and write $N = |\Lambda|$. Let $k > 0$ and define the set

= _Λ,k = { S ⊆Λ: | S |_v ≤k }.

Both $\Lambda$ and $k$ will be chosen later in the proof.  For any Banach space $(\cZ,\nm{\cdot}_{\cZ})$ define the space

_;  = { ∑_ν ∈S c_ν Ψ_ν : c_ν ∈, S ∈}.

Then, given the training data DNN-training-data, consider the problem

min_p ∈_;   1/m ∑^m_i=1 Y_i - p ∘_(X_i) ^2_,

where $\widetilde{\cY} = \cD_{\cY}(\bbR^{d_{\cY}})$. Here and throughout the proofs, we slightly abuse notation: the polynomial $p : \bbR^{\bbN} \rightarrow \bbR$, whereas $\cE_{\cX}$ has codomain $\bbR^{d_{\cX}}$. However, by construction, $p$ is independent of all but the first $d_{\cX}$ variable. Hence we may consider $p$ as a function $\bbR^{d_{\cX}} \rightarrow \bbR$.
If $\hat{p}$ is an approximate minimizer of poly-min-prob, then we define the approximation to $F$ as

F ≈F = p̂ ∘_.

Notice that $p \circ \cE_{\cX} = \cD_{\cY} \circ P \circ \cE_{\cX}$, where $P : \bbR^{d_{\cX}} \rightarrow \bbR^{d_{\cY}}$ is a vector-valued polynomial. The idea exploited in this proof is to construct a class of DNNs $\cN$ such that (i) all such polynomials $P$ are approximated by members of $\cN$ and (ii) the polynomial training problem poly-min-prob is approximated by the DNN training problem DNN-training-prob.
The first step in this analysis is therefore to analyze the polynomial training problem poly-min-prob.

§.§ Supporting lemmas

We first require the following lemma, which relates the $L^{\infty}_{\varrho}$-norm of a polynomial to its Pettis $L^2_{\varrho}$-norm.

[Nikolskii inequality for polynomials]

Let $(\cZ , \nm{\cdot}_{\cZ})$ be any Banach space and $p = \sum_{\bm{\nu} \in S} c_{\bm{\nu}} \Psi_{\bm{\nu}}$ for some finite set $S \subset \cF$, where $c_{\bm{\nu}} \in \cZ$. Then

p_L^∞_ϱ(D;) ≤√(|S|_u) p_L^2_ϱ(D;).

By definition

p_L^∞_ϱ(D;) = p_L^∞_ϱ(D;) = sup_z^* ∈B(^*) z^*(p)_L^∞_ϱ(D).

Fix $z^* \in B(\cZ^*)$ and write $z^*(p) = \sum_{\bm{\nu} \in S} z^*(c_{\bm{\nu}}) \Psi_{\bm{\nu}}$. By the triangle inequality and the definition of the weights intrinsic-weights, we have

z^*(p)_L^∞_ϱ(D) ≤∑_ν ∈S |z^*(c_ν)|Ψ_ν_L^∞_ϱ(D) = ∑_ν ∈S u_ν |z^*(c_ν)|.

We now apply the Cauchy–Schwarz inequality, weighted-card 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;) ≤√(|S|_u)  sup_z^* ∈B(^*) z^*(p)_L^2_ϱ(D) = √(|S|_u) p_L^2_ϱ(D;),

as required.

Next, we require the following bound on the Lipschitz constant of a multivariate polynomial.

[Lipschitz continuity for polynomials]

Let $(\cZ , \nm{\cdot}_{\cZ})$ be any Banach space and suppose that $p = \sum_{\bm{\nu} \in S} c_{\bm{\nu}} \Psi_{\bm{\nu}}$ for some finite $S \subset \cF$, where $c_{\bm{\nu}} \in \cZ$. Then the Lipschitz constant of $p$ satisfies

Lip(p; B^∞(), ) ≤1/2 √(|S|_v)·p_L^2_ϱ(D;) ,

where $B^{\infty}(\bbN) = \{ \bm{x} \in \bbR^{\bbN} : \nm{\bm{x}}_{\infty} \leq 1 \}$ is the unit ball of $\ell^{\infty}(\bbN)$.

Fix $z^* \in B(\cZ^*)$ and let $\tilde{p} = z^*(p) = \sum_{\bm{\nu} \in S} z^*(c_{\bm{\nu}}) \Psi_{\bm{\nu}} = :  \sum_{\bm{\nu} \in S} \tilde{c}_{\bm{\nu}} \Psi_{\bm{\nu}} $ be the corresponding scalar-valued polynomial. Let $\bm{x} , \bm{x}' \in D $. Then the mean value theorem gives that

p̃(x') - p̃(x) = ∑^∞_i=1 (x'_i - x_i) ∂p̃/∂x_i (t x + (1-t) x')

for some $0 \leq t \leq 1$. Since $B^{\infty}(\bbN) \equiv D$, this and the fact that $D$ is convex give that

| p̃(x') - p̃(x) | ≤x' - x_∞ ∑^∞_i=1 ∂p̃ / ∂x_i _L^∞_ϱ(D).

We now consider the terms in the sum separately. Using the definition of the Legendre polynomials (see <ref>), we have

∂p̃ / ∂x_i = ∑_ν ∈S c̃_ν u_ν P'_ν_i(x_i) ∏_j ≠i P_ν_j(x_j) .

The unnormalized Legendre polynomials satisy $1 = P_{n}(1) = \nm{P_n}_{L^{\infty}([-1,1])}$ and $n(n+1)/2 = P'_n(1) = \nm{P'_n}_{L^{\infty}([-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_{\bm{\nu}}$ and $v_{\bm{\nu}}$ (see intrinsic-weights and v-weights), we have

(ν_i + 1) ≤∏_j ∈ (2 ν_j+1) = u^2_ν

∑^∞_i=1 ν_i ≤∏_j ∈ (2 ν_j+1) = u^2_ν.

∑^∞_i=1  u_ν ν_i (ν_i+1) ≤u^5_ν ≤v_ν.

Here we also used the fact that $\bm{u} \geq \bm{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 poly-lip-holder now gives

| p̃(x') - p̃(x) | ≤1/2 p̃_L^2_ϱ(D)√(| S|_v) x' - x_∞ .

We now recall that $\tilde{p} = z^*(p)$ and $z^* \in B(\cZ^*)$ was arbitrary to get

p(x') - p(x) _     =  sup_z^* ∈B(^*) | z^*(p(x') - p(x)) |

≤1/2 sup_z^* ∈B(^*) z^*(p)_L^2_ϱ(D) √(| S|_v) x' - x_∞

= 1/2  p_L^2_ϱ(D;)  √(| 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. Note that in this and subsequent results, we consider the measure defined in (A.III) as a measure over the infinite cube $D = [-1,1]^{\bbN}$ that does only involves the first $d_{\cX}$ variables.

[Closeness of $L^2$ norms]

Let $(\cZ , \nm{\cdot}_{\cZ})$ be any Banach space, $\varsigma$ be the measure defined in (A.I) and $\tilde\varsigma$ be the measure defined in (A.III). Suppose that $f \in L^2_{\varsigma}(D;\cZ)$ is Lipschitz continuous with constant $L = \mathrm{Lip}(f ; B^{\infty}(\bbN) , \cZ) < \infty$ and that $f$ depends only on its first $d_{\cX}$ variables.
Then $f \in L^2_{\tilde\varsigma}(D;\cZ)$ and

f_L^2_ς(D;) - δ≤f_L^2_ς̃(D;) ≤f_L^2_ς(D;) + δ,

δ= L ·ι_d_ - ι_d_ ∘_ ∘_ _L^2_μ(; ℓ^∞()).

Fix $z^* \in B(\cZ^*)$, let $g = z^*(f) \in L^2_{\varsigma}(D)$ and set

a = g_L^2_ς̃(D),    b = g_L^2_ς(D).

Notice that

|g(x) - g(x')|≤z^*_^* f(x) - f(x')_ ≤L x - x'_∞.

Now, with slight abuse of notation, $g(\iota(X)) = g(\iota_{d_{\cX}}(X))$ due to the assumption on $f$. Therefore, using this and the Cauchy–Schwarz inequality,

a^2 - b^2     = ∫_D (g(ι_d_ ∘_ ∘_(X)))^2 - (g(ι_d_(X)))^2 μ(X)

≤L ∫_D ι_d_(X) - ι_d_ ∘_ ∘_(X)_∞ ( |g(ι_d_(X)) |  + | g(ι_d_ ∘_ ∘_) | ) μ(X)

≤L ι_d_ - ι_d_ ∘_ ∘__L^2_μ(; ℓ^∞()) ( a + b ).

We deduce that $a - b \leq \delta$. Since $z^*$ was arbitrary, we get

f_L^2_ς̃(D;) = sup_z^* ∈B(^*) z^*(f)_L^2_ς̃(D) ≤sup_z^* ∈B(^*) z^*(f)_L^2_ς(D) + δ= f_L^2_ς(D;) + δ,

which gives the upper bound.
The same argument applied to $b^2 - a^2$ also gives the lower bound.

[Norm equivalences for polynomials]

Let $(\cZ , \nm{\cdot}_{\cZ})$ be any Banach space, $\varsigma$ be the measure defined in (A.I), $\tilde\varsigma$ be the measure defined in (A.III) and $S \subset \cF$ with $\supp(\bm{\nu}) \in \{1,\ldots,d_{\cX}\}$, $\forall \bm{\nu} \in S$. Suppose that

√(|S|_v) ·ι_d_ - ι_d_ ∘_ ∘__L^2_μ(; ℓ^∞())

for some sufficiently small universal constant $c > 0$. Then the norm equivalence

p_L^2_ϱ(D;) ≲p_L^2_ς̃(D;) ≲p_L^2_ϱ(D;)

holds for all $p = \sum_{\bm{\nu} \in S} c_{\bm{\nu}} \Psi_{\bm{\nu}}$, where $c_{\bm{\nu}} \in \cZ$.

Combining Lemmas <ref> and <ref>, we see that

p_L^2_ς(D;) - δ≤p_L^2_ς̃(D;) ≤p_L^2_ς(D;) + δ,

where $\delta = \frac12 \tnm{p}_{L^2_{\varrho}(D;\cZ)} \sqrt{|S|_{\bm{v}}} \nm{\iota_{d_{\cX}} - \iota_{d_{\cX}} \circ \widetilde{\cD}_{\cX} \circ \widetilde{\cE}_{\cX}}_{L^{2}_{\mu}(\cX ; \ell^{\infty}(\bbN))} \leq c \tnm{p}_{L^2_{\varrho}(D;\cZ)} / 2$. By (A.I), there are constants $c_1 \geq c_2 > 0$ such that

c_1 p_L^2_ϱ(D;) ≤p_L^2_ς(D;) ≤c_2 p_L^2_ϱ(D;).

We now take $c = c_1 / 2$.

§.§ Analysis of poly-min-prob

We now analyze poly-min-prob. Our analysis relies on the following result, which shows an error bound subject to a certain discrete metric inequality PS-lower-cond.

[Discrete metric inequality implies error bounds]

Let $\cQ : \cY \rightarrow \widetilde{\cY}  =  \cD_{\cY}(\bbR^{d_{\cY}})$ be a bounded linear operator and $\pi_{\cQ} =  \nm{\cQ}_{\cY \rightarrow \cY} $. Suppose that

p-q_𝖽𝗂𝗌𝖼,ς̃ ≥αmax{ p-q_L^2_ϱ(D;) , p-q_L^2_ς̃(D;) },  ∀p ,q ∈_;

for some $\alpha > 0$.
Then, for any $F \in L^{\infty}_{\mu}(\cX ; \cY)$, $p \in \cP_{\cS;\widetilde{\cY}}$ and $q \in \cP_{\cS;\cY}$, we have

F - p ∘__L^2_μ(;) ≤     F - ∘F_L^2_μ(;) + α^-1 p - ∘q_𝖽𝗂𝗌𝖼,ς̃

+ π_ F - q ∘__L^2_μ(; )

F - p ∘__L^∞_μ(;) ≤     F - ∘F_L^∞_μ(;) + √(2k)/αp - ∘q_𝖽𝗂𝗌𝖼,ς̃

+ π_ F - q ∘__L^∞_μ(; ).

By the triangle inequality and properties of $\cQ$, we have

F - p ∘__L^2_μ(; )

≤F - ∘F_L^2_μ(; )
+ p ∘_ - ∘q ∘__L^2_μ(; )
+∘F - ∘q ∘__L^2_μ(; )

≤F - ∘F_L^2_μ(; )
+ p ∘_ - ∘q ∘__L^2_μ(; )
+π_ F - q ∘__L^2_μ(; ) .

Now, since $p, Q \circ q \in \cP_{S;\widetilde{\cY}}$, the second term can be bounded by

p ∘_ - ∘q ∘__L^2_μ(; )
= p - ∘q_L^2_ς̃(; ) ≤α^-1 p-∘q_𝖽𝗂𝗌𝖼,ς̃.

This yields the first result.

For the second result, we once more write

F - p ∘__L^∞_μ(; )
≤    F - ∘F_L^∞_μ(; )
+ p ∘_ - ∘q ∘__L^∞_μ(; )

+π_ F - q ∘__L^∞_μ(; ) .

For the second term, we use (A.III) to write

p ∘_ - ∘q ∘__L^∞_μ(; ) = p - ∘q_L^∞_ς̃(D;) ≤p - ∘q_L^∞_ϱ(D;).

Notice that $p - \cQ \circ q$ is a polynomial supported in a set $S$ with $|S|_{\bm{u}} \leq 2 k$.
We now apply Lemma <ref> and PS-lower-cond to obtain

p ∘_ - ∘q ∘__L^∞_μ(; ) ≤√(2k) p - ∘q_L^2_ϱ(D;) ≤√(2k)/αp-∘q_𝖽𝗂𝗌𝖼,ς̃.

This gives the result.

[Error bound for polynomial minimizers]

Let $\cQ : \cY \rightarrow \widetilde{\cY} =  \cD_{\cY}(\bbR^{d_{\cY}})$ be a bounded linear operator, $\pi_{\cQ} = \nm{\cQ}_{\cY \rightarrow \cY}$ and suppose that PS-lower-cond holds. Let $\widehat{F}$ be as in Fhat-poly-min for some $(\sigma,\tau)$-approximate minimizer $\hat{p}$ of poly-min-prob. Then

F - F_L^2_μ(;) ≤     F - ∘F_L^2_μ(;) + σ+1/α F - ∘F_𝖽𝗂𝗌𝖼,μ

+ π_ (  F - q ∘__L^2_μ(;) +σ+1/α F - q ∘__𝖽𝗂𝗌𝖼,μ)

+ τ/ α + σ+1/α√(m) E_2;

F - F_L^∞_μ(;) ≤
F - ∘F_L^∞_μ(;) + √(2k)(σ+1)/α F - ∘F_𝖽𝗂𝗌𝖼,μ

+ π_ ( F - q ∘__L^∞_μ(; ) + √(2 k)(σ+1)/α F - q ∘__𝖽𝗂𝗌𝖼,μ )

+ √(2k)τ/α + √(2 k) (σ+1)/α√(m)  E_2;

for all $q \in \cP_{\cS;\cY}$, where $\bm{E} = (E_i)^{m}_{i=1} \in \cY^m$ is the (Banach-valued) vector of noise terms.

We apply the previous lemma with $p = \hat{p}$. This gives

F - F_L^2_μ(;)  ≤      F - ∘F_L^2_μ(;) + α^-1 p̂ - ∘q_𝖽𝗂𝗌𝖼,ς̃

+ π_ F - q ∘__L^2_μ(;)

F - F_L^∞_μ(;) ≤     F - ∘F_L^∞_μ(;) + √(2k)/αp̂ - ∘q_𝖽𝗂𝗌𝖼,ς̃

+ π_ F - q ∘__L^∞_μ(; ).

Consider the second term. We have

p̂ - ∘q_𝖽𝗂𝗌𝖼,ς̃      = p̂ ∘_ - ∘q ∘_ _𝖽𝗂𝗌𝖼,μ
≤F  - ∘q ∘__𝖽𝗂𝗌𝖼,μ + F - p̂ ∘_ _𝖽𝗂𝗌𝖼,μ.

Consider the second term of this expression. By the triangle inequality and the facts that $\hat{p}$ is a $(\sigma,\tau)$-minimizer and $\cQ \circ q$ is feasible, we obtain

F - p̂ ∘_ _𝖽𝗂𝗌𝖼,μ     ≤√(1/m ∑^m_i=1 Y_i - p̂ ∘_^2_ ) + 1/√(m)E_2;

≤σ√(1/m ∑^m_i=1 Y_i - ∘q ∘_^2_ ) + τ+ 1/√(m)E_2;

≤σF - ∘q ∘__𝖽𝗂𝗌𝖼,μ + τ+ σ+1/√(m)E_2;

Therefore, we get

p̂ - ∘q_𝖽𝗂𝗌𝖼,ς̃ ≤(σ+1) F - ∘q ∘__𝖽𝗂𝗌𝖼,μ +τ+  σ+1/√(m)E_2;.

We now estimate the first term in this expression as follows:

F - ∘q ∘__𝖽𝗂𝗌𝖼,μ     ≤F - ∘F_𝖽𝗂𝗌𝖼,μ + ∘F - ∘q ∘__𝖽𝗂𝗌𝖼,μ

≤F - ∘F_𝖽𝗂𝗌𝖼,μ + π_ F -  q ∘__𝖽𝗂𝗌𝖼,μ .

Therefore, we conclude that

p̂ - ∘q_𝖽𝗂𝗌𝖼,ς̃ ≤(σ+1)  F - ∘F_𝖽𝗂𝗌𝖼,μ + (σ+1) π_ F -  q ∘__𝖽𝗂𝗌𝖼,μ + τ+ σ+1/√(m)E_2;.

Combining this with F-hatF-poly-split-1 now gives the result.

§.§ Ensuring PS-lower-cond holds with high probability

For the proof of the next lemma and subsequent steps of the proof, we let $\Lambda = \{\bm{\nu}_1,\ldots,\bm{\nu}_N \}$ and define the matrix

A =(Ψ_ν_j (_(X_i))/√(m) )^m,N_i,j=1 ∈^m ×N.

Let $0 < \epsilon < 1$, $0 < \delta < 1$, $k > 0$ and suppose that

m ≥c_0 ·δ^-2  ·k ·( log(N) ·log^2(k/δ) + log(2/ϵ))

for some universal constant $c_0 > 0 $. Let $T = \{ \bm{c} \in \bbR^N : \nm{\bm{c}}_2 =1,\ \nm{\bm{c}}_{0,\bm{v}} \leq k \}$ and

θ_+ = sup{ Ac^2_2 : c ∈T },  θ_- = inf{ Ac^2_2 : c ∈T }.

A c^2_2 ≥( θ_- - (1+θ_+)c_1 δ)c^2_2 ,  ∀c ∈^N, c_0,u ≤k.

with probability at least $1-\epsilon$, where $c_1 > 0$ is a universal constant.

The proof is based on <cit.>. Since $\nm{\bm{c}}_{1,\bm{v}} \leq \sqrt{\nm{\bm{c}}_{0,\bm{v}} } \nm{\bm{c}}_2 \leq \sqrt{k}$, we see that

T ⊆{ c ∈^N : c_1,v ≤√(k) }.

Define the random vector $\bm{X} \in \bbR^N$ by $\bm{X} = (\Psi_{\bm{\nu}_i}(\cE_{\cX}(X)) )^{N}_{i=1}$ for $X \sim \mu$. Observe that

| Xe_j | = | Ψ_ν_i(_(X)) | ≤u_ν_i ≤v_ν_i

by (A.III) and the definition of the weights. Therefore, by <cit.>, if

m ≥c_0 ·δ^-2 ·k ·log(N) ·log^2(k/c_1),

it holds that