Parametric Complexity Bounds for Approximating PDEs with Neural Networks
Year: 2022

Abstract:
Recent experiments have shown that deep networks can approximate solutions to high-dimensional PDEs, seemingly escaping the curse of dimensionality. However, questions regarding the theoretical basis for such approximations, including the required network size remain open. In this paper, we investigate the representational power of neural networks for approximating solutions to linear elliptic PDEs with Dirichlet boundary conditions. We prove that when a PDE’s coefficients are representable by small neural networks, the parameters required to approximate its solution scale polynomially with the input dimension d and proportionally to the parameter counts of the coefficient networks. To this end, we develop a proof technique that simulates gradient descent (in an appropriate Hilbert space) by growing a neural network architecture whose iterates each participate as subnetworks in their (slightly larger) successors, and converge to the solution of the PDE. We bound the size of the solution showing a polynomial dependence on d and no dependence on the volume of the domain.

heading: 1 Introduction
text: A partial differential equation (PDE) relates a multivariate function defined over some domain to its partial derivatives. Typically, one’s goal is to solve for the (unknown) function, often subject to additional constraints, such as the function’s value on the boundary of the domain. PDEs are ubiquitous in both the natural and social sciences, where they model such diverse processes as heat diffusion [7, 32], fluid dynamics [2, 39], and financial markets [6, 10]. Because most PDEs of interest lack closed-form solutions, computational approximation methods remain a vital and an active field of research [1]. For low-dimensional functions, dominant approaches include the finite differences and finite element methods [26], which discretize the domain. After partitioning the domain into a mesh, these methods solve for the function value at the vertices of the mesh. However, these techniques scale exponentially with the input dimension, rendering them unsuitable for high-dimensional problems.
Following breakthroughs in deep learning for approximating high-dimensional functions in such diverse domains as computer vision [21, 33] and natural language processing [3, 9, 40], a burgeoning line of research leverages neural networks to approximate solutions to PDEs. This line of work has produced promising empirical results for common PDEs such as the Hamilton-Jacobi-Bellman and Black-Scholes equations [17, 16, 37]. Because they do not explicitly discretize the domain, and given their empirical success on high-dimensional problems, these methods appear not to suffer the curse of dimensionality. However, these methods are not well understood theoretically, leaving open questions about when they are applicable, what their performance depends on, and just how many parameters are required to approximate the solution to a given PDE.
Over the past three years, several theoretical works have investigated questions of representational power under various assumptions. Exploring a variety of settings, [22], [16], and [19], proved that the number of parameters required to approximate a solution to a PDE exhibits a less than exponential
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
dependence on the input dimension for some special parabolic PDEs that admit straightforward analysis. [15] consider elliptic PDEs with Dirichlet boundary conditions. However, their rate depends on the volume of the domain, and thus can have an implicit exponential dependence on dimension (e.g., consider a hypercube with side length greater than one).
In this paper, we focus on linear elliptic PDEs with Dirichlet boundary conditions, which are prevalent in science and engineering (e.g., the Laplace and Poisson equations). Notably, linear elliptic PDEs define the steady state of processes like heat diffusion and fluid dynamics. Our work asks: Question. How many parameters suffice to approximate solutions to a linear elliptic PDE up to a specified level of precision using a neural network?
We show that when the coefficients of the PDE are expressible as small neural networks (note that PDE coefficients are functions), the number of parameters required to approximate the PDE’s solution is proportional to the number of parameters required to express the coefficients. Furthermore, we show that the number of parameters depends polynomially on the dimension and does not depend upon the volume of the domain.

heading: 2 Overview of Results
text: To begin, we formally define linear elliptic PDEs. Definition 1 (Linear Elliptic PDE [11]). Linear elliptic PDEs with Dirichlet boundary condition can be expressed in the following form:{
(Lu) (x) ≡ (−div (A∇u) + cu) (x) = f(x),∀x ∈ Ω, u(x) = 0,∀x ∈ ∂Ω,
where Ω ⊂ Rd is a bounded open set with a boundary ∂Ω. Further, for all x ∈ Ω, A : Ω→ Rd×d is a matrix-valued function, s.t. A(x) 0, and c : Ω→ R, s.t. c(x) > 0. 1
We refer to A and c as the coefficients of the PDE. The divergence form in Definition 1 is one of two canonical ways to define a linear elliptic PDE [11] and is convenient for several technical reasons (see Section 4). The Dirichlet boundary condition states that the solution takes a constant value (here 0) on the boundary ∂Ω.
Our goal is to express the number of parameters required to approximate the solution of a PDE in terms of those required to approximate its coefficients A and c. Our key result shows: Theorem (Informal). If the coefficients A, c and the function f are approximable by neural networks with at most N parameters, the solution u? to the PDE in Definition 1 is approximable by a neural network with O ( dlog( 1 )N ) parameters.
This result, formally expressed in Section 5, may help to explain the practical efficacy of neural networks in approximating solutions to high-dimensional PDEs with boundary conditions [37, 27]. To establish this result, we develop a constructive proof technique that simulates gradient descent (in an appropriate Hilbert space) through the very architecture of a neural network. Each iterate, given by a neural network, is subsumed into the (slightly larger) network representing the subsequent iterate. The key to our analysis is to bound both (i) the growth in network size across consecutive iterates; and (ii) the total number of iterates required.
Organization of the paper We introduce the required notation along with some mathematical preliminaries on PDEs in Section 4. The problem setting and formal statement of the main result are provided in Section 5. Finally, we provide the proof of the main result in Section 6.

heading: 3 Prior Work
text: Among the first papers to leverage neural networks to approximate solutions to PDEs with boundary conditions are [23], [24], and [29]. However, these methods discretize the input space and thus are not
1Here, div denotes the divergence operator. Given a vector field F : Rd → Rd, div(F ) = ∇·F = ∑d
i=1 ∂Fi ∂xi
suitable for high-dimensional input spaces. More recently, mesh-free neural network approaches have been proposed for high-dimensional PDEs [17, 34, 35], achieving impressive empirical results in various applications. [37] design a loss function that penalizes failure to satisfy the PDE, training their network on minibatches sampled uniformly from the input domain. They also provide a universal approximation result, showing that for sufficiently regularized PDEs, there exists a multilayer network that approximates its solution. However, they do not comment on the complexity of the neural network or how it scales with the input dimension. [20] also prove universal approximation power, albeit with networks of size exponential in the input dimension. Recently, [16, 19] provided a better-thanexponential dependence on the input dimension for some special parabolic PDEs, for which the simulating a PDE solver by a neural network is straightforward.
Several recent works [5, 22, 28, 27] show (experimentally) that a single neural network can solve for an entire family of PDEs. They approximate the map from a PDE’s parameters to its solution, potentially avoiding the trouble of retraining for every set of coefficients. Among these, only [22] provides theoretical grounding. However, they assume the existence of a finite low-dimensional space with basis functions that can approximate this parametric map—and it is unclear when this would obtain. Our work proves the existence of such maps, under the assumption that the family of PDEs has coefficients described by neural networks with a fixed architecture (Section 7).
In the work most closely related to ours, [15] provides approximation rates polynomial in the input dimension d for the Poisson equation (a special kind of linear elliptic PDE) with Dirichlet boundary conditions. They introduce a walk-on-the-sphere algorithm, which simulates a stochastic differential equation that can be used to solve a Poisson equation with Dirichlet boundary conditions (see, e.g. [31]’s Theorem 9.13). The rates provided in [15] depend on the volume of the domain, and thus depend, implicitly, exponentially on the input dimension d. Our result considers the boundary condition for the PDE and is independent of the volume of the domain. Further, we note that our results are defined for a more general linear elliptic PDE, of which the Poisson equation is a special case.

heading: 4 Notation and Definitions
text: We now introduce several key concepts from PDEs and some notation. For any open set Ω ⊂ Rd, we denote its boundary by ∂Ω and denote its closure by Ω̄ := Ω ∪ ∂Ω. By C0(Ω), we denote the space of real-valued continuous functions defined over the domain Ω. Furthermore, for k ∈ N, a function g belongs to Ck(Ω) if all partial derivatives ∂αg exist and are continuous for any multi-index α, such that |α| ≤ k. Finally, a function g ∈ C∞(Ω) if g ∈ Ck(Ω) for all k ∈ N. Next, we define several relevant function spaces:
Definition 2. For any k ∈ N ∪ {∞}, Ck0 (Ω) := {g : g ∈ Ck(Ω), supp(g) ⊂ Ω}. Definition 3. For a domain Ω, the function space L2(Ω) consists of all functions g : Ω → R, s.t. ‖g‖L2(Ω) <∞ where ‖g‖L2(Ω) = (∫ Ω |g(x)|2dx ) 1 2 . This function space is equipped with the inner product
〈g, h〉L2(Ω) = ∫
Ω
g(x)h(x)dx.
Definition 4. For a domain Ω and a function g : Ω → R, the function space L∞(Ω) is defined analogously, where ‖g‖L∞(Ω) = inf{c ≥ 0 : |g(x)| ≤ c for almost all x ∈ Ω}. Definition 5. For a domain Ω and m ∈ N, we define the Hilbert space Hm(Ω) as
Hm(Ω) := {g : Ω→ R : ∂αg ∈ L2(Ω), ∀α s.t. |α| ≤ m} Furthermore, Hm(Ω) is equipped with the inner product, 〈g, h〉Hm(Ω) = ∑ |α|≤m ∫ Ω (∂αg)(∂αh)dx and the corresponding norm
‖g‖Hm(Ω) =  ∑ |α|≤m ‖∂αg‖2L2(Ω)  12 . Definition 6. The closure of C∞0 (Ω) in Hm(Ω) is denoted by Hm0 (Ω).
Informally, Hm0 (Ω) is the set of functions belonging to H m(Ω) that can be approximated by a sequence of functions φn ∈ C∞0 (Ω). This also implies that if a function g ∈ Hm0 (Ω), then g(x) = 0 for all x ∈ ∂Ω. This space (particularly with m = 1) is often useful when analyzing elliptic PDEs with Dirichlet boundary conditions. Definition 7 (Weak Solution). Given the PDE in Definition 1, if f ∈ L2(Ω), then a function u : Ω→ R solves the PDE in a weak sense if u ∈ H10 (Ω) and for all v ∈ H10 (Ω), we have∫
Ω
(A∇u · ∇v + cuv) dx = ∫
Ω
fvdx (1)
The left hand side of (1) is also equal to 〈Lu, v〉L2(Ω) for all u, v ∈ H10 (Ω) (see Lemma 6), whereas, following the definition of the L2(Ω) norm, the right side is simply 〈f, v〉L2(Ω). Having introduced these preliminaries, we now introduce some important facts about linear PDEs that feature prominently in our analysis. Proposition 1. For the PDE in Definition 1, if f ∈ L2(Ω) the following hold:
1. The solution to Equation (1) exists and is unique.
2. The weak solution is also the unique solution of the following minimization problem:
u? = argmin v∈H10 (Ω) J(v) := argmin v∈H10 (Ω)
{ 1
2 〈Lv, v〉L2(Ω) − 〈f, v〉L2(Ω)
} . (2)
This proposition is standard (we include a proof in the Appendix, Section A.1 for completeness) and states that there exists a unique solution to the PDE (referred to as u?), which is also the solution we get from the variational formulation in (2). In this work, we introduce a sequence of functions that minimizes the loss in the variational formulation. Definition 8 (Eigenvalues and Eigenfunctions, [11]). Given an operator L, the tuples (λ, ϕ)∞i=1, where λi ∈ R and ϕi ∈ H10 (Ω) are (eigenvalue, eigenfunction) pairs that satisfy Lϕ = λϕ, for all x ∈ Ω. Since ϕ ∈ H10 (Ω), we know that ϕ|∂Ω = 0. The eigenvalue can be written as
λi = inf u∈Xi 〈Lu, u〉L2(Ω) ‖u‖2L2(Ω) , (3)
where Xi := span{ϕ1, . . . , ϕi}⊥ = {u ∈ H10 (Ω) : 〈u, ϕj〉L2(Ω) = 0 ∀j ∈ {1, · · · , i}} and 0 < λ1 ≤ λ2 ≤ · · · . Furthermore, we define by Φk the span of the first k eigenfunctions of L, i.e., Φk := span{ϕ1, · · · , ϕk}.
We note that since the operator L is self-adjoint and elliptic (in particular, L−1 is compact), the eigenvalues are real and countable. Moreover, the eigenfunctions form an orthonormal basis of H10 (Ω) (see [11], Section 6.5).

heading: 5 Main Result
text: Before stating our results, we provide the formal assumptions on the PDEs of interest:
Assumptions:
(i) Smoothness: We assume that the coefficient A ∈ Ω → Rd×d is a symmetric matrix valued function, i.e., A = (aij(x)) and aij(x) ∈ L∞(Ω) is three times differentiable for all i, j ∈ [d]. Furthermore, the function c ∈ L∞(Ω) is twice differentiable and c(x) ≥ ζ > 0 for all x ∈ Ω . We define a constant C := (2d2 + 1) max { maxα:|α|≤3 maxi,j ‖∂αaij‖L∞(Ω),maxα:|α|≤2 ‖∂αc‖L∞(Ω) } . Further, the func-
tion f ∈ L2(Ω) is infinitely differentiable and we assume there exists a function fspan ∈ Φk, such that for any multi-index α, ‖∂αf − ∂αfspan‖L2(Ω) ≤ span.
(ii) Ellipticity: There exist constants M ≥ m > 0 such that, for all x ∈ Ω and ξ ∈ Rd,
m‖ξ‖2 ≤ d∑
i,j=1
aij(x)ξiξj ≤M‖ξ‖2.
(iii) Neural network approximability: There exist neural networks Ã and c̃ with NA, Nc ∈ N parameters, respectively, that approximate the functions A and c, i.e., ‖A− Ã‖L∞(Ω) ≤ A and ‖c − c̃‖L∞(Ω) ≤ c, for small A, c ≥ 0. We assume that for all u ∈ H10 (Ω) the operator L̃ defined as,
L̃u = −div(Ã∇u) + c̃u. (4) is elliptic with (λ̃i, ϕ̃i)∞i=1 (eigenvalue, eigenfunction) pairs. We also assume that there exists a neural network fnn ∈ C∞ with Nf ∈ N parameters such that for any multi-index α, ‖∂αf − ∂αfnn‖L2(Ω) ≤ nn. By Σ, we denote the set of all (infinitely differentiable) activation functions used by networks Ã, c̃, and fnn. By Σ′, we denote the set that contains all the n-th order derivatives of the activation functions in Σ, ∀n ∈ N0
Intuitively, ellipticity of L in a linear PDE Lu = f is analogous to positive definiteness of a matrix Q ∈ Rd in a linear equation Qx = k, where x, k ∈ Rd. In (iii), we assume that the coefficients A and c, and the function f can be approximated by neural networks. While this is true for any smooth functions given sufficiently large NA, Nc, Nf , our results are most interesting when these quantities are small (e.g. subexponential in the input dimension d). For many PDEs used in practice, approximating the coefficients using small neural networks is straightforward. For example, in heat diffusion (whose equilibrium is defined by a linear elliptic PDE) A(x) defines the conductivity of the material at point x. If the conductivity is constant, then the coefficients can be written as neural networks with O(1) parameters.
Intuitively, our assumption (iii) that there exists an fspan in H10 (Ω) that lies in the span of the lowlying k eigenfunctions of the operator L can be thought of as a smoothness condition on f . For instance, if L = −∆ (the Laplacian operator), the Dirichlet form satisfies 〈Lu,u〉L2(Ω)‖u‖2
L2(Ω)
= ‖∇u‖L2(Ω) ‖u‖L2(Ω) ,
so eigenfunctions corresponding to higher eigenvalues tend to exhibit a higher degree of spikiness. The reader can also think of the eigenfunctions corresponding to larger k as Fourier basis functions corresponding to higher frequencies.
Finally, in (i) and (iii), while the requirement that the function pairs (f , fnn) and (f , fspan) are close not only in their values, but their derivatives as well is a matter of analytical convenience, our key results do not necessarily depend on this precise assumption. Alternatively, we could replace this assumption with similar (but incomparable) conditions: e.g., we can also assume closeness of the values and a rapid decay of the L2 norms of the derivatives. We require control over the derivatives because our method’s gradient descent iterations involve repeatedly applying the operator L to f—which results in progressively higher derivatives.
We can now formally state our main result: Theorem 1 (Main Theorem). Consider a linear elliptic PDE satisfying Assumptions (i)-(iii), and let u? ∈ H10 (Ω) denote its unique solution. If there exists a neural network u0 ∈ H10 (Ω) with N0 parameters, such that ‖u? − u0‖L2(Ω) ≤ R, for some R <∞, then for every > 0, there exists a neural network u with size
O ( d2T (N0 +NA) + T (Nf +Nc) ) such that ‖u? − u ‖L2(Ω) ≤ + ̃ where,
̃ := span λ1
+ δ
λ1
23/2‖f‖L2(Ω) γ − δ
+δ‖u?‖L2(Ω)+(max{1, T 2Cη})T ( span + nn +
23/2δ‖f‖L2(Ω) γ − δ
)
and κ = λ̃k+λ̃1 λ̃k−λ̃1 , η = 2 λ̃1+λ̃k
, T = O (
log(R/ ) log κ ) , and δ = 1min{m/ A,ζ/ c} and α is a multi-index.
Furthermore, the activation functions used in u belong to the set Σ ∪ Σ′ ∪ {ρ} where ρ(y) = y2 for all y ∈ R is the square activation function.
This theorem shows that given an initial neural network u0 ∈ H10 (Ω) containing N0 parameters, we can recover a neural network that is close to the unique solution u?. The number of parameters in u depend on how close the initial estimate u0 is to the solution u?, and N0. This results in a trade-off, where better approximations may require more parameters, compared to a poorer approximation with fewer parameters.
Note that can be taken arbitrarily close to 0, while ̃ is a “bias” error term that does not go to 0. The first three terms in the expression for ̃ result from bounding the difference between the solutions to the equations Lu = f and L̃u = fspan, whereas the third term is due to difference between f and fnn and the fact that our proof involves simulating the gradient descent updates with neural networks.
The term T = O (
log(R/ ) log(κ)
) comes from the fact that we are simulating T steps of a gradient
descent-like procedure on a strongly convex loss with condition number κ to reach an -approximate optimum. The parameters λk and λ1 can be thought of as the effective Lipschitz and strong-convexity constants of the loss. Finally, to give a sense of what R looks like, we show in Corollary 2 (see Section B in the Appendix) that if u0 is initialized to be identically zero then R = ‖f‖L2(Ω) λ1 .
We make few remarks about the theorem statement: Remark 1. While we state our convergence results in L2(Ω) norm, our proof works for the H10 (Ω) norm as well. This is because in the space defined by the top-k eigenfunctions of the operator L, L2(Ω) and H10 (Ω) norm are equivalent (shown in Proposition 2). Further, note that even though we have assumed that u? ∈ H10 (Ω) is the unique solution of (1) from the boundary regularity condition, we have that u? ∈ H2(Ω) (see [11], Theorem 4 in Chapter 6). This ensures that the solution u? is twice differentiable as well. Remark 2. To get a sense of the scale of λ1 and λk, when L = −∆ (the Laplacian operator), the eigenvalue λ1 = infu∈H10 (Ω) ‖∇u‖L2(Ω) ‖u‖L2(Ω)
= 1Cp , where Cp is the Poincaré constant (see Theorem 2 in Appendix). For geometrically well-behaved sets Ω (e.g. convex sets with a strongly convex boundary, like a sphere), Cp is even dimension-independent. Further from the Weyl’s law operator ([11], Section 6.5) we have
lim k→∞
λ d/2 k
k =
(2π)d
vol(Ω)α(d)
where α(d) is the volume of a unit ball in d dimensions. So, if vol(Ω) ≥ 1/α(d), λk grows as O(k2/d), which is a constant so long as log k d. Remark 3. The choice of activation functions (in particular, the requirement that Σ only contains differentiable functions) is for mathematical convenience. Namely, by standard results from approximation theory [41, 4] one can approximate a neural network with one choice of nonlinearity via a (comparably sized) neural network with another choice of nonlinearity (under mild assumptions) by incurring a dimension-independent increase in size. Thus, if the coefficients are approximable by a neural network with a non-differentiable activation function, they can also be approximated by a slightly larger network with a differentiable activation over any compact domain. Similarly, the activations on the network resulting from Theorem 1 can be replaced by a different activation at the expense of a slight blowup in size. Details are included in Appendix C.4.

heading: 6 Proof of Main Result
text: First, we provide some intuition behind the proof, via an analogy between a uniformly elliptic operator and a positive definite matrix in linear algebra. We can think of finding the solution to the equation Lu = f for an elliptic L as analogous to finding the solution to the linear system of equations Qx = k, where Q is a d×d positive definite matrix, and x and k are d-dimensional vectors. One way to solve such a linear system is by minimizing the strongly convex function ‖Qx − b‖2 using gradient descent. Since the objective is strongly convex, after O(log(1/ )) gradient steps, we reach an -optimal point in an l2 sense.
Our proof uses a similar strategy. First, we show that for the operator L, we can define a sequence of functions that converge to an -optimal function approximation (in this case in the L2(Ω) norm) after O(log(1/ ) steps—similar to the rate of convergence for strongly convex functions. Next, we inductively show that each iterate in the sequence can be approximated by a small neural network. More precisely, we show that given a bound on the size of the t-th iterate ut, we can, in turn, upper bound the size of the (t + 1)-th iterate ut+1 because the update transforming ut to ut+1 can be simulated by a small neural network (Lemma 4). These iterations look roughly like ut+1 ← ut − η(Lut − f), and we use a “backpropagation” lemma (Lemma 7) which bounds the size of the derivative of a neural network.

heading: 6.1 Defining a Convergent Sequence
text: The rough idea is to perform gradient descent in L2(Ω) [30, 12, 13], and define a convergent sequence whose iterates converge to u? in L2(Ω) norm (and following Remark 1, in H10 (Ω) as well). However, there are two obstacles to defining the iterates as simply ut+1 ← ut − η(Lut − f), (1) L is unbounded—so the standard way of choosing a step size for gradient descent (roughly the ratio of the minimum and maximum eigenvalues of L) would imply choosing a step size η = 0, and (2) L does not necessarily preserve the boundary conditions, so if we start with ut ∈ H10 (Ω), it may be that Lut − f does not even lie in H10 (Ω). We resolve both issues by restricting the updates to the span of the first k eigenfunctions of L. More concretely, as shown in Lemma 1, if a function u in Φk, then the function Lu will also lie in Φk. We also show that within the span of the first k eigenfunctions, L is bounded (with maximum eigenvalue λk), and can therefore be viewed as an operator from Φk to Φk. Further, we use fspan instead of f in our updates, which now have the form ut+1 ← ut − η(Lut − fspan). Since fspan belongs to Φk, for a ut in Φk the next iterate ut+1 will now remain in Φk. Continuing the matrix analogy, we can choose the usual step size of η = 2λ1+λk . Precisely, we show:
Lemma 1. Let L be an elliptic operator. Then, for all v ∈ Φk it holds:
1. Lv ∈ Φk.
2. λ1‖v‖L2(Ω) ≤ 〈Lv, v〉L2(Ω) ≤ λk‖v‖L2(Ω)
3. ∥∥∥(I − 2λk+λkL)u∥∥∥L2(Ω) ≤ λk−λ1λk+λ1 ‖u‖L2(Ω)
The proof of this lemma is provided in Section C.1.
Further, note that we will use a slight variant of the updates and instead set ut+1 ← ut−η(L̃u−f̃span) as the iterates of the convergent sequence. Here, the operator L̃ (defined in (4)) has the neural network approximations of A and c as its coefficients, and f̃span is a function that lies in span of the first k eigenfunctions of L̃ (denoted by Φ̃k), such that ‖fspan− f̃span‖L2(Ω) is small (for an exact statement, see Lemma 11 in the Appendix). In Section 6.2, we will see that updates defined thusly will be more convenient to simulate via a neural network.
The sequence defined so far satisfies two important properties. First, the convergence point of the sequence and u∗, the solution to the original PDE, are not too far from each other. Concretely:
Lemma 2. Assume that ũ?span is the solution to the PDE L̃u = f̃span, where f̃span : H10 (Ω) → R and f̃span ∈ Φ̃k. Given Assumptions (i)-(iii), we have ‖u? − ũ?span‖L2(Ω) ≤ , such that = span λ1 + δλ1 23/2‖f‖L2(Ω) γ−δ + δ‖ũ ? span‖L2(Ω), where γ = 1λk − 1 λk+1 and δ = 1min{m/ A,ζ/ c} .
The proof for Lemma 2 is provided in the Appendix (Section D.1). Each of the three terms in the final error captures different sources of perturbation: the first term comes from approximating f by fspan; the second term comes from applying Davis-Kahan [8] to bound the “misalignment” between the eigenspaces Φk and Φ̃k (hence, the appearance of the eigengap between the k and (k + 1)-st eigenvalue of L−1); the third term is a type of “relative” error bounding the difference between the solutions to the PDEs Lu = f̃span and L̃u = f̃span.
The second property of the sequence of functions is that they converge exponentially fast. The rate is characterized in the following lemma:
Lemma 3 (Convergence of gradient descent in L2). Let ũ?span denote the unique solution to the PDE L̃u = f̃span, where f̃span ∈ Φ̃k, and the operator L̃ satisfies the conditions in Lemma 1. Then for any u0 ∈ H10 (Ω) such that u0 ∈ Φ̃k, we define the sequence
ut+1 ← ut − 2
λ̃1 + λ̃k (L̃ut − f̃span) (t ∈ N) (5)
where for all t ∈ N, ut ∈ H10 (Ω). Then for any ≥ 0 we have ‖uT − ũ?span‖L2(Ω) ≤ after T iterations where, T ≥ log (‖u0−ũ?span‖L2(Ω) ) log ( λ̃k+λ̃1 λ̃k−λ̃1
) The proof for Lemma 3 is similar to the analysis of the convergence time of gradient descent for strongly convex losses and can be found in Section C.2 of the Appendix. Finally, combining the results from Lemma 2 and Lemma 3 via triangle inequality, we have:
‖u? − uT ‖L2(Ω) ≤ ‖u? − ũ?span‖L2(Ω) + ‖ũ?span − uT ‖L2(Ω) and the first term on the RHS subsumes the first three summands of ̃ defined in Theorem 1.

heading: 6.2 Approximating iterates by neural networks
text: In Lemma 3, we show that there exists a sequence of functions (5) which converge fast to a function close to u?. The next step in the proof is to approximate the iterates by neural networks.
The main idea is as follows. Suppose first the iterates ut+1 = ut − η(L̃ut − f̃span) are such that f̃span is exactly representable as a neural network. Then, the iterate ut+1 can be written in terms of three operations performed on ut, a and f : taking derivatives, multiplication and addition. Moreover, if g is representable as a neural network with N parameters, the coordinates of the vector∇g can be represented by a neural network withO(N) parameters. This is a classic result (Lemma 7), essentially following from the backpropagation algorithm. Finally, addition or multiplication of two functions representable as neural networks with sizes N1, N2 can be represented as neural networks with size O(N1 +N2) (see Lemma 8).
Using these facts, we can write down a recurrence upper bounding the size of neural network approximation ut+1, denoted by ût+1, in terms of the number of parameters in ût (which is the neural network approximation to ut). Formally, we have: Lemma 4 (Recursion Lemma). Given the Assumptions (i)-(iii), consider the update equation
ût+1 ← ût − 2
λ̃1 + λ̃k
( L̃ût − fnn ) (6)
If at step t, ût : Rd → R is a neural network with Nt parameters, then the function ût+1 is a neural network with O(d2(NA +Nt) +Nt +Nf̃ +Nc) parameters. Proof. Expand the update ût+1 ← ût − η ( L̃ût − fnn ) as follows:
ût+1 ← ût − η  d∑ i,j=1 ãij∂ij ût + d∑ j=1 ( d∑ i=1 ∂iãij ) ∂j ût + c̃ût − fnn  . Using Lemma 7, ∂ij ût, ∂j ût and ∂iãij can be represented by a neural network with O(Nt), O(Nt) and O(NA) parameters, respectively. Further, ∂iãij∂ju and ãij∂ij û can be represented by a neural network with O(NA +Nt) parameters, and c̃ût can be represented by a network with O(Nt +Nc) parameters, from Lemma 8. Hence ût+1 can be represented in O(d2(NA +Nt) +Nf +Nc +Nt) parameters. Note that, throughout the entire proofs O hides independent constants.
Combining the results of Lemma 3 and Lemma 4, we can get a recurrence for the number of parameters required to represent the neural network ût:
Nt+1 ≤ d2Nt + d2NA +Nt +Nf̃ +Nc
Unfolding this recurrence, we get NT ≤ d2TN0 + d 2(dT−1) d2−1 NA + T (Nf ) +Nc).
Hence, the total number of parameters required for a neural network to approximate a solution to a PDE of the form in Definition 1
O ( d2 log  ‖u0−ũ?span‖L2(Ω)  log κ (N0 +NA) + log (‖u0−ũ?span‖L2(Ω) ) log κ (Nf +Nc) )
Finally, we have to deal with the fact that f̃span is not exactly a neural network, but only approximately so. The error due to this discrepancy can be characterized through the following lemma:
Lemma 5 (Error using fnn). Consider the update equation in (6), where fnn is a neural network with Nf . Then the neural network ût approximates the function ut such that ‖ut − ût‖L2(Ω) ≤ (t) nn where (t)nn is
O ( (max{1, t2ηeC})t ( span + nn +
23/2δ‖f‖L2(Ω) γ − δ )) where δ = 1min{m/ A,ζ/ c} , γ = 1 λk − 1λk+1 , and α is a multi-index.
The main strategy to prove this lemma involves tracking the “residual” non-neural-network part of the iterates. Precisely, we can write the update ut+1 = ut−η(L̃ut− (fnn +r)), for a “residual” function r = f̃span − fnn. If ut was exactly a neural network, the first part of the update, ut − η(L̃ut − fnn) can be written as a neural network as in Lemma 4, and ηr can be treated as an error. Thus, in order to bound the error for using fnn instead of f̃span, we maintain a decomposition of ut into a neural network part, and a residual part, and inductively bound the total residual part. Given the recurrent structure of our updates, at each iteration t, there will be an accumulation of the earlier t− 1 applications of the operator L̃, that results in an increasing number of higher order derivatives of r to be bounded at each step. This is why we require that fnn is close to f not only in the L2 sense but also in terms of their higher order derivatives.

heading: 7 Applications to Learning Operators
text: A number of recent works attempt to simultaneously approximate the solutions for an entire family of PDEs by learning a parametric map that takes as inputs (some representation of) the coefficients of a PDE and returns its solution [5, 28, 27]. For example, given a set of observations that {aj , uj}Nj=1, where each aj denotes a coefficient of a PDE with corresponding solution uj , they learn a neural network G such that for all j, uj = G(aj). Our parametric results provide useful insights for why simultaneously solving an entire family of PDEs with a single neural network G is possible in the case of linear elliptic PDEs.
Consider the case where the coefficients aj in the family of PDEs are given by neural networks with a fixed architecture, but where each instance of a PDE is characterized by a different setting of the weights in the models representing the coefficients. Lemma 4 shows that each iteration of our sequence (5) constructs a new network containing both the current solution and the coefficient networks as subnetworks. We can view our approximation as not merely approximating the solution to a single PDE but to every PDE in the family, by treating the coefficient networks as placeholder architectures whose weights are provided as inputs. Thus, our construction provides a parametric map between the coefficients of an elliptic PDE in this family and its solution.

heading: 8 Conclusion and Future Work
text: We derive parametric complexity bounds for neural network approximations for solving linear elliptic PDEs with Dirichlet boundary conditions, whenever the coefficients can be approximated by are neural networks with finite parameter counts. By simulating gradient descent in function spaces using neural networks, we construct a neural network that approximates the solution of a PDE. We show that the number of parameters in the neural network depends on the parameters required to represent the coeffcients and has a poly(d) dependence on the dimension of the input space, therefore avoiding the curse of dimensionality.
An immediate open question is related to the tightening our results: our current error bound is sensitive to the neural network approximation lying close to Φk which could be alleviated by relaxing (by adding some kind of “regularity” assumptions) the dependence of our analysis on the first k eigenfunctions. Further, the dependencies in the exponent of d on R and κ in parametric bound may also be improvable. Finally, the idea of simulating an iterative algorithm by a neural network to derive a representation-theoretic result is broadly applicable, and may be a fertile ground for further work, both theoretically and empirically, as it suggest a particular kind of weight tying.

heading: 9 Acknowledgement
text: This paper is based upon work funded and supported by the Department of Defense under contract FA8702-15-D-0002

Review summaries:
The authors attempt to bound a neural network's complexity that approximates a solution to a self-adjoint elliptic PDE. The paper lacks mathematical rigor on the PDE theory side: the most important one is that their gradient descent step doesn't make sense from a functional analysis point-of-view (see main review). The authors need to compare the complexity of their neural network to the bound obtained by combining standard PDE regularity theory (see Evan's Section 6.3) with the complexity of neural networks for approximating smooth functions. Since the main theorem doesn't take into account the smoothness of the boundary of \Omega, I suspect that the latter gives much stronger bounds for many practical domains. 

The authors show that the solution of a linear elliptic PDE (with boundary condition) can be represented with a neural network, whose number of parameters is at most polynomial in the input dimension d, under the assumption that the coefficients of the PDE are given by NNs with poly(d) parameters.  This is achieved by first showing that gradient descent with a particular choice of stepsize, provides a sequence in $L_2$, that converges to the solution of the PDE. Then the iterates of the sequence are approximated by NNs in a recursive manner, and counting the number of parameters that are added at each iteration.  In essence, the result aims to help the understanding of how learning a solution of a high-dimensional PDE from samples might be computationally feasible using parametrized families of functions, like NNs.

The paper proved an $\epsilon$ approximation theory showing that the neural network approximation for linear elliptical PDEs is bounded by the number of NN parameters required for approximating the PDE coefficients as well as the dimension of the domain. The proof is mainly based on an operator theoretic point of view and the bounds are obtained by approximating gradient descent using neural networks. The results also provide insights into the recent success of using NN to approximate high dimensional as well as parametric PDEs.