A Theoretical Analysis of Deep Neural Networks and Parametric PDEs

We derive upper bounds on the complexity of ReLU neural networks approximating the solution maps of parametric partial differential equations. In particular, without any knowledge of its concrete shape, we use the inherent low dimensionality of the solution manifold to obtain approximation rates which are significantly superior to those provided by classical neural network approximation results. Concretely, we use the existence of a small reduced basis to construct, for a large variety of parametric partial differential equations, neural networks that yield approximations of the parametric solution maps in such a way that the sizes of these networks essentially only depend on the size of the reduced basis.


Introduction
In this work, we analyze the suitability of deep neural networks (DNNs) for the numerical solution of parametric problems. Such problems connect a parameter space with a solution state space via a so-called parametric map [51]. One special case of such a parametric problem arises when the parametric map results from solving a partial differential equation (PDE) and the parameters describe physical or geometrical constraints of the PDE such as, for example, the shape of the physical domain, boundary conditions, or a source term. Applications that lead to these problems include modeling unsteady and steady heat and mass transfer, acoustics, fluid mechanics, or electromagnetics [33].
Solving a parametric PDE for every point in the parameter space of interest individually typically leads to two types of problems. First, if the number of parameters of interest is excessive-a scenario coined many-query application-then the associated computational complexity could be unreasonably high. Second, if the computation time is severely limited, such as in real-time applications, then solving even a single PDE might be too costly.
A core assumption to overcome the two issues outlined above is that the solution manifold, i.e., the set of all admissible solutions associated with the parameter space, is inherently low-dimensional. This assumption forms the foundation for the so-called reduced basis method (RBM). A reduced basis discretization is then a (Galerkin) projection on a low-dimensional approximation space that is built from snapshots of the parametrically induced manifold [60].
Constructing the low-dimensional approximation spaces is typically computationally expensive because it involves solving the PDEs for multiple instances of parameters. These computations take place in a so-called offline phase-a step of pre-computation, where one assumes to have access to sufficiently powerful computational resources. Once a suitable low-dimensional space is found, the cost of solving the associated PDEs for a new parameter value is significantly reduced and can be performed quickly and online, i.e., with limited resources [5,56]. We will give a more thorough introduction to RBMs in Sect. 2. An extensive survey of works on RBMs, which can be traced back to the seventies and eighties of the last century (see, for instance, [22,49,50]), is beyond the scope of this paper. We refer, for example, to [33,Chapter 1.1], [16,29,57] and [12,Chapter 1.9] for (historical) studies of this topic.
In this work, we show that the low-dimensionality of the solution manifold also enables an efficient approximation of the parametric map by DNNs. In this context, the RBM will be, first and foremost, a tool to model this low-dimensionality by acting as a blueprint for the construction of the DNNs.

Statistical Learning Problems
The motivation to study the approximability of parametric maps by DNNs stems from the following similarities between parametric problems and statistical learning problems: Assume that we are given a domain set X ⊂ R n , n ∈ N and a label set Y ⊂ R k , k ∈ N. Further assume that there exists an unknown probability distribution ρ on X × Y .
Given a loss function L : Y × Y → R + , the goal of a statistical learning problem is to find a function f , which we will call prediction rule, from a hypothesis class H ⊂ {h : X → Y } such that the expected loss E (x,y)∼ρ L( f (x), y) is minimized [14]. Since the probability measure ρ is unknown, we have no direct access to the expected loss. Instead, we assume that we are given a set of training data, i.e., pairs (x i , y i ) N i=1 , N ∈ N, which were drawn independently with respect to ρ. Then one finds f by minimizing the so-called empirical loss over H . We will call optimizing the empirical loss the learning procedure.
In view of PDEs, the approach proposed above can be rephrased in the following way. We are aiming to produce a function from a parameter set to a state space based on a few snapshots only. This function should satisfy the involved PDEs as precisely as possible, and the evaluation of this function should be very efficient even though the construction of it can potentially be computationally expensive.
In the above-described sense, a parametric PDE problem almost perfectly matches the definition of a statistical learning problem. Indeed, the PDEs and the metric on the state space correspond to a (deterministic) distribution ρ and a loss function. Moreover, the snapshots are construed as the training data, and the offline phase mirrors the learning procedure. Finally, the parametric map is the prediction rule.
One of the most efficient learning methods nowadays is deep learning. This method describes a range of learning procedures to solve statistical learning problems where the hypothesis class H is taken to be a set of DNNs [24,40]. These methods outperform virtually all classical machine learning techniques in sufficiently complicated tasks from speech recognition to image classification. Strikingly, training DNNs is a computationally very demanding task that is usually performed on highly parallelized machines. Once a DNN is fully trained, however, its application to a given input is orders of magnitudes faster than the training process. This observation again reflects the offline-online phase distinction that is common in RBM approaches. Based on the overwhelming success of these techniques and the apparent similarities of learning problems and parametric problems, it appears natural to apply methods from deep learning to statistical learning problems in the sense of (partly) replacing the parameter-dependent map by a DNN. Very successful advances in this direction have been reported in [17,34,39,42,58,68].

Our Contribution
In the applications [17,34,39,42,58,68] mentioned above, the combination of DNNs and parametric problems seems to be remarkably efficient. In this paper, we present a theoretical justification of this approach. We address the question to what extent the hypothesis class of DNNs is sufficiently broad to approximately and efficiently represent the associated parametric maps. Concretely, we aim at understanding the necessary number of parameters of DNNs required to allow a sufficiently accurate approximation. We will demonstrate that depending on the target accuracy the required number of parameters of DNNs essentially only scales with the intrinsic dimension of the solution manifold, in particular, according to its Kolmogorov N -widths. We outline our results in Sect. 1.2.1. Then, we present a simplified exposition of our argument leading to the main results in Sect. 1.2.2.

Approximation Theoretical Results
The main contributions of this work are given by an approximation result with DNNs based on ReLU activation functions. Here, we aim to learn a variation of the parametric map where Y is the parameter space and H is a Hilbert space. In our case, the parameter space will be a compact subset of R p for some fixed, but possibly large p ∈ N, i.e., we consider the case of finitely supported parameter vectors.
We assume that there exists a basis of a high-fidelity discretization of H which may potentially be quite large. Let u y be the coefficient vector of u y with respect to the high-fidelity discretization. Moreover, we assume that there exists a RB approximating u y sufficiently accurately for every y ∈ Y. Theorem 4.3 then states that, under some technical assumptions, there exists a DNN that approximates the discretized solution map Y y → u y up to a uniform error of > 0, while having a size that is polylogarithmical in , cubic in the size of the reduced basis, and at most linear in the size of the high-fidelity basis.
This result highlights the common observation that if a low-dimensional structure is present in a problem, then DNNs are able to identify it and use it advantageously. Concretely, our results show that a DNN is sufficiently flexible to benefit from the existence of a reduced basis in the sense that its size in the complex task of solving a parametric PDE does not or only weakly depend on the high-fidelity discretization and mainly on the size of the reduced basis.
The main result is based on four pillars that are described in detail in Sect. 1.2.2: First, we show that DNNs can efficiently solve linear systems, in the sense that, if supplied with a matrix and a right-hand side, a moderately sized network outputs the solution of the inverse problem. Second, the reduced-basis approach allows reformulating the parametric problem, as a relatively small and parametrized linear system. Third, in many cases, the map that takes the parameters to the stiffness matrices with respect to the reduced basis and right-hand side can be very efficiently represented by DNNs. Finally, the fact that neural networks are naturally compositional allows combining the efficient representation of linear problems with the NN implementing operator inversion.
In practice, the approximating DNNs that we show to exist need to be found using a learning algorithm. In this work, we will not analyze the feasibility of learning these DNNs. The typical approach here is to apply methods based on stochastic gradient descent. Empirical studies of this procedure in the context of learning deep neural networks were carried out in [34,39,42,58,68]. In particular, we mention the recent study in [23], which analyzes precisely the setup described in this work and finds a strong impact of the approximation-theoretical behavior of DNNs on their practical performance.

Simplified Presentation of the Argument
In this section, we present a simplified outline of the arguments leading to the approximation result described in Sect. 1.2.1. In this simplified setup, we think of a ReLU neural network (ReLU NN) as a function where L ∈ N, T 1 , . . . , T L are affine maps, and : R → R, (x) := max{0, x} is the ReLU activation function which is applied coordinate-wise in (1.2). We call L the number of layers of the NN. Since T are affine linear maps, we have for all x ∈ dom T that T (x) = A (x) + b for a matrix A and a vector b . We define the size of the NN as the number of nonzero entries of all A and b for ∈ {1, . . . , L}. This definition will later be sharpened and extended in Definition 3.1.
1. As a first step, we recall the construction of a scalar multiplication operator by ReLU NNs due to [69]. This construction is based on two observations. First, defining g : [0, 1] → [0, 1], g(x) := min{2x, 2−2x}, we see that g is a hat function. Moreover, multiple compositions of g with itself produce saw-tooth functions. We set, for s ∈ N, g 1 := g and g s+1 := g • g s . It was demonstrated in [69] that The second observation for establishing an approximation of a scalar multiplication by NNs is that we can write g(x) = 2 (x) − 4 (x − 1/2) + 2 (x − 2) and therefore g s can be exactly represented by a ReLU NN. Given that g s is bounded by 1, it is not hard to see that f n converges to the square function exponentially fast for n → ∞. Moreover, f n can be implemented exactly as a ReLU NN by previous arguments. Finally, the parallelogram identity, where the additional log 2 term in m inside the brackets appears since each of the approximations of the sum needs to be performed with accuracy /m. It is well known that the Neumann series m s=0 A s converges exponentially fast to (Id R d − A) −1 for m → ∞. Therefore, under suitable conditions on the matrix A, we can construct a NN inv that approximates the inversion operator, i.e., the map A → A −1 up to accuracy > 0. This NN has size O(d 3 log q 2 (1/ )) for a constant q > 0. This is made precise in Theorem 3.8. 4. The existence of inv and the emulation of approximate matrix-vector multiplications yield that there exists a NN that for a given matrix and right-hand side approximately solves the associated linear system. Next, we make two assumptions that are satisfied in many applications as we demonstrate in Sect. 4.2: • The map from the parameters to the associated stiffness matrices of the Galerkin discretization of the parametric PDE with respect to a reduced basis can be well approximated by NNs. • The map from the parameters to the right-hand side of the parametric PDEs discretized according to the reduced basis can be well approximated by NNs.
From these assumptions and the existence of inv and a ReLU NN emulating a matrix-vector multiplication, it is not hard to see that there is a NN that approximately implements the map from a parameter to the associated discretized solution with respect to the reduced basis. If the reduced basis has size d and the implementations of the map yielding the stiffness matrix and the right-hand side are sufficiently efficient, then, by the construction of inv , the resulting NN has size O(d 3 log q 2 (1/ )). We call this NN rb . 5. Finally, we build on the construction of rb to establish the result of Sect. 1.2.1.
First of all, let D be the size of the high-fidelity basis. If D is sufficiently large, then every element from the reduced basis can be approximately represented in the highfidelity basis. Therefore, one can perform an approximation to a change of bases by applying a linear map V ∈ R D×d to a vector with respect to the reduced basis. The first statement of Sect. 1.2.1 now follows directly by considering the NN V • rb . Through this procedure, the size of the NN is increased to O(d 3 log q 2 (1/ )) + d D). The full argument is presented in the proof of Theorem 4.3.

Potential Impact and Extensions
We believe that the results of this article have the potential to significantly impact the research on NNs and parametric problems in the following ways: • Theoretical foundation: We offer a theoretical underpinning for the empirical success of NNs for parametric problems which was observed in, e.g., [34,39,42,58,68]. Indeed, our result, Theorem 4.3, indicates that properly trained NNs are as efficient in solving parametric PDEs as RBMs if the complexity of NNs is measured in terms of free parameters. On a broader level, linking deep learning techniques for parametric PDE problems with approximation theory opens the field up to a new direction of thorough mathematical analysis. • Understanding the role of the ambient dimension: It has been repeatedly observed that NNs seem to offer approximation rates of high-dimensional functions that do not deteriorate exponentially with increasing dimension [24,45]. In this context, it is interesting to identify the key quantity determining the achievable approximation rates of DNNs. Possible explanation for approximation rates that are essentially independent from the ambient dimension have been identified if the functions to be approximated have special structures such as compositionality [48,55], or invariances [45,53]. In this article, we identify the highly problem-specific notion of the dimension of the solution manifold as a key quantity determining the achievable approximation rates by NNs for parametric problems. We discuss the connection between the approximation rates that NNs achieve and the ambient dimension in detail in Sect. 5. • Identifying suitable architectures: One question in applications is how to choose the right NN architectures for the associated problem. Our results show that NNs of sufficient depth and size are able to produce very efficient approximations. Nonetheless, it needs to be mentioned that our results do not yield a lower bound on the number of layers and thus it is not clear whether deep NNs are indeed necessary.
This work is a step toward establishing a theory of deep learning-based solutions of parametric problems. However, given the complexity of this field, it is clear that many more steps need to follow. We outline a couple of natural further questions of interest below: • General parametric problems: Below we restrict ourselves to coercive, symmetric, and linear parametric problems with finitely many parameters. There exist many extensions to, e.g., non-coercive, non-symmetric, or nonlinear problems [10,11,25,38,67,70], or to infinite parameter spaces, see, e.g., [2,4]. It would be interesting to see if the methods proposed in this work can be generalized to these more challenging situations. • Bounding the number of snapshots: The interpretation of the parametric problem as a statistical learning problem has the convenient side-effect that various techniques have been established to bound the number of necessary samples N , such that the empirical loss (1.1) is very close to the expected loss. In other words, the generalization error of the minimizer of the learning procedure is small, meaning that the prediction rule performs well on unseen data. (Here, the error is measured in a norm induced by the loss function and the underlying probability distribution.) Using these techniques, it is possible to bound the number of snapshots required for the offline phase to achieve a certain fidelity in the online phase. Estimates of the generalization error in the context of high-dimensional PDEs have been deduced in, e.g., [7,19,20,26,59]. • Special NN architectures: This article studies the feasibility of standard feedforward NNs. In practice, one often uses special architectures that have proved efficient in applications. First and foremost, almost all NNs used in applications are convolutional neural networks (CNNs) [41]. Hence, a relevant question is to what extent the results of this work also hold for such architectures. It was demonstrated in [54] that there is a direct correspondence between the approximation rates of CNNs and that of standard NNs. Thus, we expect that the results of this work translate to CNNs. Another successful architecture is that of residual neural networks (ResNets) [32]. These neural networks also admit skip-connections, i.e., do not only connect neurons in adjacent layers. This architecture is by design at least as powerful as a standard NN and hence inherits all approximation properties of standard NNs. • Necessary properties of neural networks: In this work, we demonstrate the attainability of certain approximation rates by NNs. It is not clear if the presented results are optimal or if there are specific necessary assumptions on the architectures, such as a minimal depth, a minimal number of parameters, or a minimal number of neurons per layer. For approximation results of classical function spaces, such lower bounds on specifications of NNs have been established, for example, in [9,28,53,69]. It is conceivable that the techniques in these works can be transferred to the approximation tasks described in this work. • General matrix polynomials: As outlined in Sect. 1.2.2, our results are based on the approximate implementation of matrix polynomials. Naturally, this construction can be used to define and construct a ReLU NN based functional calculus. In other words, for any d ∈ N and every continuous function f that can be well approximated by polynomials, we can construct a ReLU NN which approximates the map A → f (A) for any appropriately bounded matrix A. A special instance of such a function of interest is given by f (A) := e tA , t > 0, which is analytic and plays an important role in the treatment of initial value problems. • Numerical studies: In a practical learning problem, the approximation-theoretical aspect only describes one part of the problem. Two further central factors are the data generation and the optimization process. It is conceivable that in comparison with these issues, approximation theoretical considerations only play a negligible role. To understand the extent to which the result of this paper is relevant for applications, comprehensive studies of the theoretical setup of this work should be carried out. A first one was published recently in [23].

Related Work
In this section, we give an extensive overview of works related to this paper. In particular, for completeness, we start by giving a review of approximation theory of NNs without an explicit connection to PDEs. Afterward, we will see how NNs have been employed for the solution of PDEs.

Review of Approximation Theory of Neural Networks
The first and most fundamental results on the approximation capabilities of NNs were universality results. These results claim that NNs with at least one hidden layer can approximate any continuous function on a bounded domain to arbitrary accuracy if they have sufficiently many neurons [15,35]. However, these results do not quantify the required sizes of NNs to achieve these rates. One of the first results in this direction was given in [6]. There, a bound on the sufficient size of NNs with sigmoidal activation functions approximating a function with finite Fourier moments is presented. Further results describe approximation rates for various smoothness classes by sigmoidal or even more general activation functions [43,44,46,47].
For the non-differentiable activation function ReLU, first rates of approximation were identified in [69] for classes of smooth functions, in [53] for piecewise smooth functions, and in [27] for oscillatory functions. Moreover, NNs mirror the approximation rates of various dictionaries such as wavelets [62], general affine systems [9], linear finite elements [31], and higher-order finite elements [52].

Neural Networks and PDEs
A well-established line of research is that of solving high-dimensional PDEs by NNs assuming that the NN is the solution of the underlying PDE, e.g., [7,19,20,30,36,37,59,63]. In this regime, it is often possible to bound the size of the involved NNs in a way that does not scale exponentially with the underlying dimension. In that way, these results are quite related to our approaches. Our results do not seek to represent the solution of a PDE as a NN, but a parametric map. Moreover, we analyze the complexity of the solution manifold in terms of Kolmogorov N -widths. Finally, the underlying spatial dimension of the involved PDEs in our case would usually be moderate. However, the dimension of the parameter space could be immense.
One of the first approaches analyzing NN approximation rates for solutions of parametric PDEs was carried out in [61]. In that work, the analyticity of the solution map y → u y and polynomial chaos expansions with respect to the parametric variable are used to approximate the map y → u y by ReLU NNs of moderate size. Moreover, we mention the works [34,39,42,58,68] which apply NNs in one way or another to parametric problems. These approaches study the topic of learning a parametric problem but do not offer a theoretical analysis of the required sizes of the involved NNs. These results form our motivation to study the constructions of this paper.
Finally, we mention that the setup of the recent numerical study [23] is closely related to this work.

Outline
In Sect. 2, we describe the type of parametric PDEs that we consider in this paper, and we recall the theory of RBs. Section 3 introduces a NN calculus which is the basis for all constructions in this work. There we will also construct the NN that maps a matrix to its approximate inverse in Theorem 3.8. In Sect. 4, we construct NNs approximating parametric maps. First, in Theorem 4.1, we approximate the parametric maps after a high-fidelity discretization. Afterward, in Sect. 4.2, we list two broad examples where all assumptions which we imposed are satisfied.
We conclude this paper in Sect. 5 with a discussion of our results in light of the dependence of the underlying NN complexities in terms of the governing quantities.
To not interrupt the flow of reading, we have deferred all auxiliary results and proofs to the appendices.

Notation
We denote by N = {1, 2, . . .} the set of all natural numbers and define N 0 := N ∪ {0}. Moreover, for a ∈ R we set a := max{b ∈ Z : b ≤ a} and a := min{b ∈ Z : b ≥ a}. Let n, l ∈ N. Let Id R n be the identity and 0 R n be the zero vector on R n . Moreover, for A ∈ R n×l , we denote by A T its transpose, by σ (A) the spectrum of A, by A 2 its spectral norm and by A 0 := #{(i, j) : A i, j = 0}, where #V denotes the cardinality of a set V , the number of nonzero entries of A. Moreover, for v ∈ R n we denote by |v| its Euclidean norm. Let V be a vector space. Then we say that , the set of all scalar-valued, linear, continuous functions equipped with the operator norm. For a compact set ⊂ R n , we denote by C r ( ), r ∈ N 0 ∪ {∞}, the spaces of r times continuously differentiable functions, by L p ( , R n ), p ∈ [1, ∞] the R n -valued Lebesgue spaces, where we set L p ( ) := L p ( , R) and by H 1 ( ) := W 1,2 ( ) the first-order Sobolev space.

Parametric PDEs and Reduced Basis Methods
In this section, we introduce the type of parametric problems that we study in this paper. A parametric problem in its most general form is based on a map P : Y → Z, where Y is the parameter space and Z is called solution state space Z. In the case of parametric PDEs, Y describes certain parameters of a partial differential equation, Z is a function space or a discretization thereof, and P(y) ∈ Z is found by solving a PDE with parameter y.
We will place several assumptions on the PDEs underlying P and the parameter spaces Y in Sect. 2.1. Afterward, we give an abstract overview of Galerkin methods in Sect. 2.2 before recapitulating some basic facts about RBs in Sect. 2.3.

Parametric Partial Differential Equations
In the following, we will consider parameter-dependent equations given in the variational form which fulfills certain wellposedness conditions specified in Assumption 2.1, (iv) f y ∈ H * is the parameter-dependent right-hand side of (2.1), (v) u y ∈ H is the solution of (2.1).

Assumption 2.1
Throughout this paper, we impose the following assumptions on Eq. (2.1).
• The parameter set Y: We assume that Y is a compact subset of R p , where p ∈ N is fixed and potentially large.
Remark In [12, Section 1.2], it has been demonstrated that if Y is a compact subset of some Banach space V , then one can describe every element in Y by a sequence of real numbers in an affine way. To be more precise, there exist (v i ) ∞ i=0 ⊂ V such that for every y ∈ Y and some coefficient sequence c y whose elements can be bounded in absolute value by 1 there holds implying that Y can be completely described by the collection of sequences c y . In this paper, we assume these sequences c y to be finite with a fixed, but possibly large support size.
• Symmetry, uniform continuity, and coercivity of the bilinear forms: We assume that for all y ∈ Y the bilinear forms b y are symmetric, i.e., Moreover, we assume that the bilinear forms b y are uniformly continuous in the sense that there exists a constant C cont > 0 with Finally, we assume that the involved bilinear forms are uniformly coercive in the sense that there exists a constant C coer > 0 such that Hence, by the Lax-Milgram lemma (see [57, Lemma 2.1]), Eq. (2.1) is well-posed, i.e., for every y ∈ Y and every f y ∈ H * there exists exactly one u y ∈ H such that (2.1) is satisfied and u y depends continuously on f y .
• Uniform boundedness of the right-hand side: We assume that there exists a constant C rhs > 0 such that • Compactness of the solution manifold: We assume that the solution manifold

Remark
The assumption that S(Y) is compact follows immediately if the solution map y → u y is continuous. This condition is true (see [57, Proposition 5.1, Corollary In fact, there exists a multitude of parametric PDEs, for which the maps . Moreover, we refer to [57, Remark 5.2] and the references therein for a discussion under which circumstances it is possible to turn a discontinuous parameter dependency into a continuous one ensuring the compactness of S(Y).

High-Fidelity Approximations
In practice, one cannot hope to solve (2.1) exactly for every y ∈ Y. Instead, if we assume for the moment that y is fixed, a common approach toward the calculation of an approximate solution of (2.1) is given by the Galerkin method, which we will describe shortly following [33, Appendix A] and [57,Chapter 2.4]. In this framework, instead of solving (2.1), one solves a discrete scheme of the form Moreover, up to a constant, we have that u disc y is a best approximation of the solution u y of (2.1) by elements in U disc . To be more precise, by Cea's Lemma, [ Let us now assume that U disc is given. Moreover, if N := dim U disc , let (ϕ i ) N i=1 be a basis for U disc . Then, the matrix is non-singular and positive definite. The solution u disc y of (2.2) satisfies Typically, one starts with a high-fidelity discretization of the space H, i.e., one chooses a finite but potentially high-dimensional subspace for which the computed discretized solutions are sufficiently accurate for any y ∈ Y.
To be more precise, we postulate the following: We assume that, for every y ∈ Y, sup y∈Y u y − u h y H ≤ˆ for an arbitrarily small, but fixedˆ > 0. In the following, similarly as in [16], we will not distinguish between H and U h , unless such a distinction matters.
In practice, following this approach, one often needs to calculate u h y ≈ u y for a variety of parameters y ∈ Y. This, in general, is a very expensive procedure due to the high-dimensionality of the space U h . In particular, given , one needs to solve high-dimensional systems of linear equations to determine the coefficient vector u h y . A well-established remedy to overcome these difficulties is given by methods based on the theory of reduced bases, which we will recapitulate in the upcoming subsection.
Before we proceed, let us fix some notation. We denote by G : (2.4)

Theory of Reduced Bases
In this subsection and unless stated otherwise, we follow [57,Chapter 5] The starting point of this theory lies in the concept of the Kolmogorov N -width which is defined as follows. [16] For N ∈ N, the Kolmogorov N -width of a bounded subset X of a normed space V is defined by

Definition 2.3
This quantity describes the best possible uniform approximation error of X by an at most N -dimensional linear subspace of V . We discuss concrete upper bounds on W N (S(Y)) in more detail in Sect. 5. The aim of RBMs is to construct the spaces U rb in such a way that the quantity sup y∈Y dist u y , U rb is close to W d(˜ ) (S(Y)).

The identification of the basis vectors
of U rb usually happens in an offline phase in which one has considerable computational resources available and which is usually based on the determination of high-fidelity discretizations of samples of the parameter set Y. The most common methods are based on (weak) greedy procedures (see for instance [57,Chapter 7] and the references therein) or proper orthogonal decompositions (see, for instance, [57,Chapter 6] and the references therein). In the last step, an orthogonalization procedure (such as a Gram-Schmidt process) is performed to obtain an orthonormal set of basis vectors . Afterward, in the online phase, one assembles for a given input y the corresponding low-dimensional stiffness matrices and vectors and determines the Galerkin solution by solving a low-dimensional system of linear equations. To ensure an efficient implementation of the online phase, a common assumption which we do not require in this paper is the affine decomposition of (2.1), which means that there exist as well as As has been pointed out in [57,Chapter 5.7], in principal three types of reduced bases generated by RBMs have been established in the literature-the Lagrange reduced basis, the Hermite reduced basis and the Taylor reduced basis. While the most common type, the Lagrange RB, consists of orthonormalized versions of high-fidelity given expansion point y ∈ Y. In this paper, we will later assume that there exist small Note that all types of RBs discussed above satisfy this assumption. The next statement gives a (generally sharp) upper bound which relates the possibility of constructing small snapshot RBs directly to the Kolmogorov N -width.
Translated into our framework, Theorem 2.4 states that for every Remark 2. 5 We note that this bound is sharp for general X , V . However, it is not necessarily optimal for special instances of S(Y). If, for instance, W N (S(Y)) decays polynomially, then W N (S(Y)) decays with the same rate (see [8,Theorem 3 Taking the discussion above as a justification, we assume from now on that for every˜ ≥ˆ there exists a RB space D is chosen to be as small as possible, at least fulfilling ). In addition, we assume that the vectors (ψ i ) form an orthonormal system in H, which is equivalent to the fact that the columns of G 1/2 V˜ are orthonormal (see [57,Remark 4.1]). This in turn implies as well as For the underlying discretization matrix, one can demonstrate (see, for instance, [57, Moreover, due to the symmetry and the coercivity of the underlying bilinear forms combined with the orthonormality of the basis vectors (ψ i ) implying that the condition number of the stiffness matrix with respect to the RB remains bounded independently of y and the dimension d(˜ ). Additionally, the discretized right-hand side with respect to the RB is given by and, by the Bessel inequality, we have that f rb y,˜ ≤ f y H * ≤ C rhs . Moreover, let be the coefficient vector of the Galerkin solution with respect to the RB space. Then, the Galerkin solution u rb y,˜ can be written as In the following sections, we will emulate RBMs with NNs by showing that we are able to construct NNs which approximate the maps u rb ·,˜ ,ũ h ·,˜ such that their complexity depends only on the size of the reduced basis and at most linearly on D.
The key ingredient will be the construction of small NNs implementing an approximate matrix inversion based on Richardson iterations in Sect. 3. In Sect. 4, we then proceed with building the NNs the realizations of which approximate the maps u rb ·,˜ ,ũ h ·,˜ , respectively.

Neural Network Calculus
The goal of this chapter is to emulate the matrix inversion by NNs. In Sect. 3.1, we introduce some basic notions connected to NNs as well as some basic operations one can perform with these. In Sect. 3.2, we state a result which shows the existence of NNs the ReLU-realizations of which take a matrix A ∈ R d×d , A 2 < 1 as their input and calculate an approximation of Id R d − A −1 based on its Neumann series expansion. The associated proofs can be found in "Appendix A."

Basic Definitions and Operations
We start by introducing a formal definition of a NN. Afterward, we introduce several operations, such as parallelization and concatenation that can be used to assemble complex NNs out of simpler ones. Unless stated otherwise, we follow the notion of [53] where most of this formal framework was introduced. First, we introduce a terminology for NNs that allows us to differentiate between a NN as a family of weights and the function implemented by the NN. This implemented function will be called the realization of the NN.  1 , b 1 ), (A 2 , b 2 ), . . . , (A L , b L ) , If is a NN as above, K ⊂ R n , and if : R → R is arbitrary, then we define the associated realization of with activation function over K (in short, therealization of over K ) as the map R K ( ) : K → R N L such that where x L results from the following scheme: and where acts componentwise, that is, We First of all, we note that it is possible to concatenate two NNs in the following way.
be two NNs such that the input layer of 1 has the same dimension as the output layer of 2 . Then, 1 2 denotes the following L 1 + L 2 − 1 layer NN: We call 1 2 the concatenation of 1 , 2 .
In general, there is no bound on M( 1 2 ) that is linear in M( 1 ) and M( 2 ). For the remainder of the paper, let be given by the ReLU activation function, i.e., (x) = max{x, 0} for x ∈ R. We will see in the following, that we are able to introduce an alternative concatenation which helps us to control the number of nonzero weights. Toward this goal, we give the following result which shows that we can construct NNs the ReLU-realization of which is the identity function on R n . We now introduce the sparse concatenation of two NNs. We will see later in Lemma 3.6 the properties of the sparse concatenation of NNs. We proceed with the second operation that we can perform with NNs. This operation is called parallelization.
Now, let be a NN and L ∈ N such that L( ) ≤ L. Then, define the NN Finally, let˜ 1 , . . . ,˜ k be NNs which have the same input dimension and let Then, we define P ˜ 1 , . . . ,˜ k := P EL ˜ 1 , . . . , EL ˜ k .
The following lemma was established in [21,Lemma 5.4] and examines the properties of the sparse concatenation as well as of the parallelization of NNs. Lemma 3.6 [21] Let 1 , . . . , k be NNs.
(a) If the input dimension of 1 , which shall be denoted by n 1 , equals the output dimension of 2 , and n 2 is the input dimension of 2 , then (b) If the input dimension of i , denoted by n, equals the input dimension of j , for all i, j, then for the NN P 1 , 2 , . . . , k and all x 1 , . . . x k ∈ R n we have as well as

A Neural Network-Based Approach Toward Matrix Inversion
The goal of this subsection is to emulate the inversion of square matrices by realizations of NNs which are comparatively small in size. In particular, Theorem 3.8 shows that, for d ∈ N, ∈ (0, 1/4), and δ ∈ (0, 1), we are able to efficiently construct NNs 1−δ,d inv; the ReLU-realization of which approximates the map To stay in the classical NN setting, we employ vectorized matrices in the remainder of this paper. Let A ∈ R d×l . We write vec(A) := (A 1,1 , . . . , A d,1 , . . . , A 1,l , . . . , A d,l ) Moreover, for a vector v = (v 1,1 , . . . , v d,1 , . . . , v 1,d , . . . , v d,l ) T ∈ R dl we set In addition, for d, n, l ∈ N and Z > 0 we set The basic ingredient for the construction of NNs emulating a matrix inversion is the following result about NNs emulating the multiplication of two matrices. with n · (d + l)-dimensional input, dl-dimensional output such that, for an absolute constant C mult > 0, the following properties are fulfilled: Based on Proposition 3.7, we construct in "Appendix A.2" NNs emulating the map A → A k for square matrices A and k ∈ N. This construction is then used to prove the following result.

Remark 3.9
In the proof of Theorem 3.8, we approximate the function mapping a matrix to its inverse via the Neumann series and then emulate this construction by NNs. There certainly exist alternative approaches to approximating this inversion function, such as, for example, via Chebyshev matrix polynomials (for an introduction of Chebyshev polynomials, see, for instance, [65,Chapter 8.2]). In fact, approximation by Chebyshev matrix polynomials is more efficient in terms of the degree of the polynomials required to reach a certain approximation accuracy. However, emulation of Chebyshev matrix polynomials by NNs either requires larger networks than that of monomials or, if they are represented in a monomial basis, coefficients that grow exponentially with the polynomial degree. In the end, the advantage of a smaller degree in the approximation through Chebyshev matrix polynomials does not seem to set off the drawbacks described before.

Neural Networks and Solutions of PDEs Using Reduced Bases
In this section, we invoke the estimates for the approximate matrix inversion from Sect. 3.2 to approximate the parameter-dependent solution of parametric PDEs by NNs. In other words, for˜ ≥ˆ , we construct NNs approximating the maps Here, the sizes of the NNs essentially only depend on the approximation fidelity˜ and the size d(˜ ) of an appropriate RB, but are independent or at most linear in the dimension of the high-fidelity discretization D.
We start in Sect. 4.1 by constructing, under some general assumptions on the parametric problem, a NN emulating the mapsũ h ·,˜ and u rb ·,˜ . In Sect. 4.2, we verify these assumptions on two examples.

Determining the Coefficients of the Solution
Next, we present constructions of NNs the ReLU-realizations of which approximate the mapsũ h ·,˜ and u rb ·,˜ , respectively. In our main result of this subsection, the approximation error of the NN approximationũ h ·,˜ will be measured with respect to the | · | G -norm since we can relate this norm directly to the norm on H via Eq. (2.4). In contrast, the approximation error of the NN approximating u rb ·,˜ will be measured with respect to the | · |-norm due to Eq. 2.8.
As already indicated earlier, the main ingredient of the following arguments is an application of the NN of Theorem 3.8 to the matrix B rb y,˜ . As a preparation, we show in Proposition B.1 in Appendix, that we can rescale the matrix B rb y,˜ with a constant factor α := (C coer + C cont ) −1 (in particular, independent of y and d(˜ )) so that with C coer δ := αC coer We will fix these values of α and δ for the remainder of the manuscript. Next, we state two abstract assumptions on the approximability of the map B rb ·,˜ which we will later on specify when we consider concrete examples in Sect. 4.2. In addition to Assumption 4.1, we state the following assumption on the approximability of the map f rb ·,˜ . Assumption 4. 2 We assume that for every˜ ≥ˆ , > 0, and a corresponding RB Now we are in a position to construct NNs the ReLU-realizations of which approximate the coefficient mapsũ h ·,˜ , u rb ·,˜ . Theorem 4.3 Let˜ ≥ˆ and ∈ (0, α/4 · min{1, C coer }) . Moreover, define := / max{6, C rhs }, := /3 · C coer , := 3/8 · αC 2 coer and κ := 2 max {1, C rhs , 1/C coer }. Additionally, assume that Assumption 4.1 and Assumption 4.2 hold. Then, there exist NNs u,rb , and u,h , such that the following properties hold: that the constants C u L , C u M depend on the constants C coer , C cont , C rhs in the following way (recall that C coer 2C cont ≤ δ = C coer C coer +C cont ≤ 1 2 ): • C u L depends affine linearly on log 2 2 log 2 (δ/2) log 2 (1 − δ) , log 2 1 C coer + C cont , log 2 (max {1, C rhs , 1/C coer }) .
• C u M depends affine linearly on suggests that (y, x) → u y (x) can be approximated with essentially the cost of approximating the respective function in (4.1). Many basis elements that are commonly used for the high-fidelity representation can indeed be approximated very efficiently by realizations of NNs, such as, e.g., polynomials, finite elements, or wavelets [31,52,62,69].

Examples of Neural Network Approximation of Parametric Maps
In this subsection, we apply Theorem 4.3 to a variety of concrete examples in which the approximation of the coefficient maps u rb ·,˜ ,ũ h ·,˜ can be approximated by comparatively small NNs. We show that the sizes of these NNs depend only on the size of associated reduced bases by verifying Assumption 4.1 and Assumption 4.2, respectively. We will discuss to what extent our results depend on the respective ambient dimensions D, p in Sect. 5.
We will state the following examples already in their variational formulation and note that they fulfill the requirements of Assumption 2.1. We also remark that the presented examples represent only a small portion of problems to which our theory is applicable.

Example I: Diffusion Equation
We consider a special case of [57, Chapter 2.3.1] which can be interpreted as a generalized version of the heavily used example −div(a∇u) = f , where a is a scalar field (see for instance [13,61] and the references therein). Let n ∈ N, ⊂ R n , be a Lipschitz domain and H := H 1 0 ( ) = u ∈ H 1 ( ) : u| ∂ = 0 . We assume that the parameter set is given by a compact set T ⊂ L ∞ ( , R n×n ) such that for all T ∈ T and almost all x ∈ the matrix T(x) is symmetric, positive definite with matrix norm that can be bounded from above and below independently of T and x. As we have noted in Assumption 2.1, we can assume that there exist some We restrict ourselves to the case of finitely supported sequences ( To be more precise, let p ∈ N be potentially very high but fixed, let Y := [−1, 1] p and consider for y ∈ Y and some fixed f ∈ H * the parametric PDE Then, the parameter dependency of the bilinear forms is linear, hence analytic whereas the parameter dependency of the right-hand side is constant, hence also analytic, implying that W (S(Y)) decays exponentially fast. This in turn implies existence of small RBs (ψ i ) In conclusion, in this example, we have, for˜ , > 0, Theorem 4.3 hence implies the existence of a NN approximating u rb ·,˜ up to error with a size that is linear in p, polylogarithmic in 1/ , and, up to a log factor, cubic in d(˜ ). Moreover, we have shown the existence of a NN approximatingũ h ·,˜ with a size that is linear in p, polylogarithmic in 1/ , linear in D and, up to a log factor, cubic in d(˜ ).
The parameter dependency of the right-hand side is linear (hence analytic), whereas the parameter dependency of the bilinear forms is rational, hence (due to the choice ofỹ 1,2 ,ỹ 2,2 ) also analytic and W N (S(Y)) decays exponentially fast implying that we can choose d(˜ ) to depend polylogarithmically on˜ . It is now easy to see that Assump- Thus, Theorem 4.3 implies the existence of NNs approximating u rb ·,˜ up to error with a size that is polylogarithmic in 1/ , and, up to a log factor, cubic in d(˜ ). Moreover, there exist NNs approximatingũ h ·,˜ up to error with a size that is linear in D, polylogarithmic in 1/ , and, up to a log factor, cubic in d(˜ ).
For a more thorough discussion of this example (a special case of the linear elasticity equation which describes the displacement of some elastic structure under physical stress on its boundaries), we refer to [

Discussion: Dependence of Approximation Rates on Involved Dimensions
In this section, we will discuss our results in terms of the dependence on the involved dimensions. We would like to stress that the resulting approximation rates (which can be derived from Theorem 4.3) differ significantly from and are often substantially better than alternative approaches. As described in Sect. 2, there are three central dimensions that describe the hardness of the problem. These are the dimension D of the high-fidelity discretization space U h , the dimension d(˜ ) of the reduced basis space, and the dimension p of the parameter space Y. Dependence on D: Examples I and II above establish approximation rates that depend at most linearly on D; in particular, the dependence on D is not coupled to the depen-dence on . Another approach to solve these problems would be to directly solve the linear systems from the high-fidelity discretization. Without further assumptions on sparsity properties of the matrices, the resulting complexity would be O(D 3 ) plus the cost of assembling the high-fidelity stiffness matrices. Since D d(˜ ), this is significantly worse than the approximation rate provided by Theorem 4.3. Dependence on d(˜ ) : If one assembles and solves the Galerkin scheme for a previously found reduced basis, one typically needs O d(˜ ) 3 operations. By building NNs emulating this method, we achieve essentially the same approximation rate of O d(˜ ) 3 log 2 (d(˜ )) · C( ) where C( ) depends polylogarithmically on the approximation accuracy .
Note that, while having comparable complexity, the NN-based approach is more versatile than using a Galerkin scheme and can be applied even when the underlying PPDE is fully unknown as long as sufficiently many snapshots are available.

Dependence on p:
We start by comparing our result to naive NN approximation results which are simply based on the smoothness properties of the map y →ũ h ·,˜ without using its specific structure. For example, if these maps are analytic, then classical approximation rates with NNs (such as those provided by [ N and a constant c( p, n). In this case, the dependence on D is again linear, but coupled with the potentially quickly growing term − p/n . Similarly, when approximating the map y →ũ rb ·,˜ , one would obtain an approximation rate of − p/n . In addition, our approach is more flexible than the naive approach in the sense that Assumptions 4.1 and 4.2 could even be satisfied if the map y → B rb y,˜ is non-smooth. Now we analyze the dependence of our result to p in more detail. We recall from For problems of the type (2.6), where the involved maps θ q are sufficiently smooth and the right-hand side is parameter-independent, one can show (see, for instance, [1,Equation 3.17] or [51] that W N (S(Y)) (and hence also W N (S(Y))) scales like e −cN Q b for some c > 0. This implies for the commonly studied case Q b = p (such as in Example I of Sect Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

A.1 Proof of Proposition 3.7
In this subsection, we will prove Proposition 3.7. As a preparation, we first prove the following special instance under which M( 1 2 ) can be estimated by max M( 1 ), M( 2 ) .  (((a, 0)) ) ≤ M ( ).

Lemma A.1 Let be a NN with m-dimensional output and d-dimensional input. If
In particular, it holds that M ((a, 0) ) ≤ M( ). Moreover, if D ∈ R d×n such that, for every k ≤ d there is at most one l k ≤ n such that D k,l k = 0, then, for all = 1, . . . , L( ), In particular, it holds that M ( ((D, 0 R d ))) ≤ M( ).
Proof Let = (A 1 , b 1 ), . . . , (A L , b L ) , and a, D as in the statement of the lemma. Then, the result follows if and It is clear that aA L 0 is less than the number of nonzero columns of A L which is certainly bounded by A L 0 . The same argument shows that ab L 0 ≤ b L 0 . This yields (A.1). We have that for two vectors p, q ∈ R k , k ∈ N and for all μ, ν ∈ R where I (γ ) = 0 if γ = 0 and I (γ ) = 1 otherwise. Also, where, for a matrix G, G l,− denotes the l-th row of G. Moreover, we have that for all l ≤ n As a consequence, we obtain Now we are ready to prove Proposition 3.7.
As a final step, we define Z ,d,n,l mult;˜ := P completing the proof of (iii). By construction and using the fact that for any N ∈ R d×l there holds This demonstrates that (v) holds and thereby finishes the proof.

A.2 Proof of Theorem 3.8
The objective of this subsection is to prove of Theorem 3.8. Toward this goal, we construct NNs which emulate the map A → A k for k ∈ N and square matrices A. This is done by heavily using Proposition 3.7. First of all, as a direct consequence of Proposition 3.7 we can estimate the sizes of the emulation of the multiplication of two squared matrices. Indeed, there exists a universal constant C 1 > 0 such that for all d ∈ N, Z > 0, ∈ (0, 1) (vec(A), vec(B)) One consequence of the ability to emulate the multiplication of matrices is that we can also emulate the squaring of matrices. We make this precise in the following definition.
which has d 2 -dimensional input and d 2 -dimensional output. By Lemma 3.6 we have that there exists a constant C sq > C 1 such that for all d ∈ N, Z > 0, ∈ (0, 1) Our next goal is to approximate the map A → A k for an arbitrary k ∈ N 0 . We start with the case that k is a power of 2 and for the moment we only consider the set of all matrices the norm of which is bounded by 1/2.
as well as M .
By the triangle inequality, we obtain for any vec( (vec(A)) 2 + A 2 j matr R Therefore, using the triangle inequality and the fact that · 2 is a submultiplicative operator norm, we derive that (vec(A)) where the penultimate estimate follows by the induction hypothesis (A.10) and < 1/4. Hence, since · 2 is a submultiplicative operator norm, we obtain where we used A 2 j 2 ≤ 1/4 and the induction hypothesis (A.10). Applying (A. 13) and (A.12) to (A.11) yields A direct consequence of (A.14) is that The estimates (A.14) and (A.15) complete the proof of the assertions (iv) and (v) of the proposition statement. Now we estimate the size of 1/2,d 2 j+1 ; . By the induction hypothesis and Lemma 3.6(a)(i), we obtain ≤ C sq · log 2 (1/ ) + log 2 (d) + log 2 (4) + j log 2 (1/ ) which implies (i). Moreover, by the induction hypothesis and Lemma 3.6(a)(ii), we conclude that implying (ii). Finally, it follows from Lemma 3.6(a)(iii) in combination with the induction hypothesis as well Lemma 3.6(a)(iv) that which finishes the proof.
We proceed by demonstrating, how to build a NN that emulates the map A → A k for an arbitrary k ∈ N 0 . Again, for the moment we only consider the set of all matrices the norms of which are bounded by 1/2. For the case of the set of all matrices the norms of which are bounded by an arbitrary Z > 0, we refer to Corollary A.5. with d 2 -dimensional input and d 2 -dimensional output satisfying the following properties: Proof We prove the result per induction over k ∈ N 0 . The cases k = 0 and k = 1 hold trivially by defining the NNs 1/2,d 0; 1/2,d 1; For the induction hypothesis, we claim that the result holds true for all k ≤ k ∈ N.
If k is a power of two, then the result holds per Proposition A.3; thus, we can assume without loss of generality, that k is not a power of two. We define j := log 2 (k) such that, for t := k − 2 j , we have that 0 < t < 2 j . This implies that A k = A 2 j A t . Hence, by Proposition A. 3 By construction and Lemma 3.6(a)(iv), we first observe that M L Moreover, we obtain by the induction hypothesis as well as Lemma 3.6(a)(iii) in combination with Lemma 3.6(b)(iv) that This shows (iii). To show (iv), we perform a similar estimate as the one following (A.11). By the triangle inequality, (A. 16) By the construction of . We have by Lemma 3.6(a)(i) in combination with Lemma 3.6(b)(i) and by the induction hypothesis that which implies (i). Finally, we address the number of nonzero weights of the resulting NN. We first observe that, by Lemma 3.6(a)(ii),  .
Without loss of generality we assume that L := L where we have used the definition of the parallelization for the first two equalities, Lemma 3.6(b)(iii) for the third equality, Lemma 3.6(a)(ii) for the fourth inequality as well as the properties of Id d 2 ,L in combination with Proposition A.3(iii) for the last inequality. Moreover, by the definition of the parallelization of two NNs with different numbers of layers, we conclude that Combining the estimates on I , II (a), and II (b), we obtain by using the induction hypothesis that Proposition A.4 only provides a construction of a NN the ReLU-realization of which emulates a power of a matrix A, under the assumption that A 2 ≤ 1/2. We remove this restriction in the following corollary by presenting a construction of a NN Z ,d k; the ReLU-realization of which approximates the map A → A k , on the set of all matrices A the norms of which are bounded by an arbitrary Z > 0.
Corollary A.5 There exists a universal constant C pow > C sq such that for all Z > 0, d ∈ N and k ∈ N 0 , there exists some NN Z ,d k; with the following properties: as well as M Proof Let ((A 1 , b 1 ), . . . , (A L , b L )) := fulfills all of the desired properties.
We have seen how to construct a NN that takes a matrix as an input and computes a power of this matrix. With this tool at hand, we are now ready to prove Theorem 3.8.
This completes the proof.
We proceed with estimating II. It is not hard to see from Assumption  (y) ≤ + C rhs ≤ C coer + C rhs ≤ 1 + C rhs ≤ κ.