A Priori Estimation of the Approximation, Optimization and Generalization Errors of Random Neural Networks for Solving Partial Differential Equations

1 Introduction

As the development of hardcore and algorithms, deep neural networks have made remarkable progress across various fields, including computer vision [He et al., 2016], natural language processing [Devlin, 2018], reinforcement learning [Silver et al., 2016], and others. These successes have inspired researchers to explore the application of neural networks to scientific challenges, particularly in the modeling of physical systems. In the field of scientific computing, a long-standing problem is solving partial differential equations (PDEs) numerically, which can be hindered by the curse of dimensionality when using classical numerical methods. To tackle PDE-related problems, several neural network-based approaches have been introduced. Among these, Physics-Informed Neural Networks (PINNs) and the Deep Ritz Method (DRM) stand out. PINN incorporates prior information from PDEs into the training in solving forward and inverse problems of PDEs. Specifically, it encodes the PDE constraints into the design of the loss function, which restricts the constructed neural network to follow physical law characterized by the PDE. This framework's flexibility stems from its reliance solely on the PDE's form, making it adaptable for a wide range of PDEs. The DRM, on the other hand, incorporates the variational formulation that is an essential tool in traditional methods, into training the neural networks. This approach typically involves lower-order derivatives, potentially offering greater computational efficiency compared to PINNs. However, the DRM's utility is somewhat limited by the fact that not all PDEs possess a variational formulation.

∗Corresponding author: Ye Li <yeli20@nuaa.edu.cn> and Zhongyi Huang <zhongyih@tsinghua.edu.cn>

The success of neural networks can be partly attributed to their powerful approximation capabilities. The capacity of neural networks with a variety of activation functions to approximate different types of target functions has been extensively studied. This includes continuous functions [Shen et al., 2022], smooth functions [Lu et al., 2021a], as well as Sobolev functions [Belomestny et al., 2023; Yang et al., 2023b; Yang et al., 2023a; Yarotsky, 2017] and Barron functions [Barron, 1993; Lu et al., 2021b; Siegel and Xu, 2020a], among others. In short, neural networks can approximate the aforementioned function classes with arbitrary precision. Despite the widespread adoption and impressive approximation capabilities of neural networks in scientific computing, their practical application can encounter challenges. A fundamental issue arises in solving the optimization problems involving neural networks, which are typically non-convex and nonsmooth. For example, [Krishnapriyan et al., 2021] demonstrates that even for simple problems of convection, reaction, and reaction-diffusion, PINN approach only works for very easy parameter regimes and fails in more challenging physical regimes. By analyzing the loss landscape of PINN, they show that the failure is not due to the limited capacity of the neural network, but rather due to the optimization difficulties caused by the PINN's soft PDE constraint.

Due to the limitations described before, there is a growing interest in the application of random neural networks, whose hidden weights are randomly generated and only the output weights are trainable. Compared to deterministic neural networks, random neural networks can lead to optimization problems that can be efficiently solved. For instance, employing random neural networks in L2 regression problems results in least squares problems, which possess closedform solutions or can be addressed by various optimization algorithms. Because of the favorable properties, random neural networks have been successfully applied not only in traditional machine learning tasks but also in addressing PDErelated problems. In this work, we focus on the utilization of random neural networks in the framework of PINNs. First, we establish the approximation error of neural networks for Barron-type functions under the H2 norm. Subsequently, for applying random neural networks for solving a certain class of second-order elliptic PDEs with Barron-type solutions, we provide the optimization and generalization error analysis, providing a comprehensive understanding of how random neural networks perform when applied to PDEs.

1.1 Contributions

The contribution of this work can be summarized as follows.

• We provides two approximation results for two Barrontype functions via random neural networks under the Sobolev norms. Unlike previous works that only showed the existence of the prior distributions, we give the concise forms of these distributions.
• When applying random neural networks for solving certain second-order elliptic PDEs whose solutions are Barron-type, we design tailored optimization algorithms for addressing these problems. Subsequently, we perform a full error analysis that rigorously bounds the approximation, optimization, and generalization errors. In deriving the generalization bounds, we also show that projected gradient descent can achieve the faster rate O(1/n).
• We have conducted numerical experiments to validate our conclusions.

1.2 RelatedWorks

Random Neural Networks in Machine Learning. Random neural networks, also known as Extreme Learning Machines (ELMs) in the field of machine learning, have drawn significant attention due to their special training methodology and efficiency. The concept of ELMs was first introduced in [Huang et al., 2006] as a two-layer random neural network where the hidden weights are randomly assigned, and only the output weights are trained. As shown in [Huang et al., 2006], this algorithm can provide good generalization performance at extremely fast learning speed and it outperforms deterministic neural networks in regression, classification and real-world complex problems. It is also possible to extend ELM to kernel learning [Huang et al., 2011], which shows that ELM can use a wide type of feature mappings, including random hidden nodes and kernels.

RandomNeural Networks in Scientific Computing. The potential of random neural networks extends beyond traditional machine learning domains, they have also been effectively utilized in scientific computing, particularly for solving problems related to PDEs. For instance, [Dong and Li, 2021] has developed an efficient method based on domain decomposition and random neural networks to solve different types of PDEs, showing the first time when the neural networkbased methods outperform the traditional numerical methods in low dimensions. Later, this method has been extended to high dimensions [Wang and Dong, 2024], producing accurate solutions to high-dimensional PDEs.

[Gonon et al., 2023] has studied approximation based on single-hidden-layer feedforward and recurrent neural networks with randomly generated hidden weights under the L2 norm. Similar to our work, they are also based on an integral representations of the target functions. However, the prior distribution for the hidden weights remains unclear in their study. In contrast, [Neufeld and Schmocker, 2023] derived approximation rates and an explicit algorithm to learn a deterministic function by a random neural network under certain Sobolev norms, but the conditions for the target functions to be approximated are challenging to verify. The work most closely related to ours is [Gonon, 2023], which provided a full analysis of random neural networks for learning sufficiently non-degenerate Black-Scholes type models. However, their approximation results are presented under the L∞ norm, and their generalization analysis is grounded in methods for L2 regression problems. Moreover, their methods cannot be applied to our setting, which involves unsupervised learning as opposed to the supervised problems they consider. To the best of our knowledge, this is the first article that attempts to provide theoretical support for solving PDEs within the framework of PINNs using random neural networks.

1.3 Notations

For x ∈ Rd , ∥x∥p denotes its p-norm (1 ≤ p ≤ ∞). For the activation functions, we write σk(x) for the ReLUk activation function, i.e., σk(x) := [max(0, x)]k . For given probability measure P and a sequence of random variables {Xi}ni=1 distributed according to P, we denote the empirical measure of

2 Preliminaries

In this section, we provide some preliminaries about random neural networks and Barron functions, which are pivotal to our study. Throughout the paper, we only consider the twolayer random neural networks, i.e., a feedforward neural network with one hidden layer and randomly generated hidden weights. For brevity, we still call them random neural networks.

To make it more precise, for m ∈ N, let W1, · · · ,Wm be i.i.d. Rd -valued random vectors and B1, · · · ,Bm be i.i.d. real-valued random variables, where W = (W1, · · · ,Wm) and B = (B1, · · · ,Bm) are independent. Then for any Rmvalued (random) vector A = (A1, · · · ,Am), we have the random neural network

where σ : R → R is a fixed activation function.

In addressing the L2 regression with random neural networks, given training samples {(xi, yi)}ni=1 , where xi is the input and yi is the output for 1 ≤ i ≤ n, we have the following objective function to be optimized.

This results in a least squares problem that has a closed-form solution. In contrast, deterministic neural networks also require optimizing the weights of the hidden layers, which leads to a more complex optimization problem.

In this work, we mainly focus on the Barron-type functions, leveraging their integral representations as highlighted in [Siegel and Xu, 2022]. The Barron space with order s > 0 is defined as

where the infimum is over extensions fe ∈ L1 (Rd ) and ˆ fe is the Fourier transform of fe. Barron introduced this class for s = 1 and showed that two-layer neural networks with sigmoidal activation function can achieve the approximation rate O(1/√n) in the L2 norm. Although the convergence rate does not suffer the curse of dimensionality, the related optimization problems are non-convex and challenging to address. Consequently, we shift our focus toward approximating Barron functions employing random neural networks. Note that we choose 1-norm in the definition (3) just for simplicity. In the following, we assume that Ω is a subset of [0, 1]d .

3 Main Results

3.1 Approximation Results

Our approach to approximating Barron functions by random two-layer neural networks leverages the integral representations of these functions, as presented in [Siegel and Xu, 2020b]. Specifically, for f ∈ Bs (Ω), without loss of generality, assume that the infimum is attained at fe. Then for activation function σ ∈ L1 (R) with ˆσ(a) ̸= 0 and some a ̸= 0, we have

In the paper, we consider the activation function σ(t) = ∑4i=0 (-1)iC i4σ3(t + 2 - i), which is compactly supported on [-2, 2]. The conclusions can be naturally extended to sigmoidal activation functions and tanh activation functions, which are the commonly used smooth activation functions for PINNs. Specific details can be found in Remark 2. The approximation results rely on the polynomial decay condition of the activation function, i.e., |σ(k) (t)| ≤ Cp(1 + |t|)-p for 0 ≤ k ≤ 2 for some p > 1. This condition also appears in [Siegel and Xu, 2020b] for deriving the approximation rate for functions in Bs (Ω) using two-layer deterministic neural networks.

Theorem 1. Let P1 be the uniform distribution of the domain {ω ∈ Rd : ∥ω∥1 ≤ M} and P2 be the uniform distribution of the domain {b ∈ R : |b| ≤ 2M} with a constant M ≥ 2. Let W1, · · · ,Wm ∼ P1 and B1, · · · ,Bm ∼ P2, then there exists a Rm-valued, σ(W,B)-measurable vector A = (A1, · · · ,Am) such that for 3 < s < 6,

and for s ≥ 6,

where

and

≲ indicates that a universal multiplicative constant is omitted.

Remark 1. By appropriately selecting the value of M, we can obtain that

Remark 2. Extending our discussion beyond the specific activation function we have chosen, we explore other common activation functions, such as the Sigmoidal function Sig(x) defined as (1 + e-x )-1 and the Hyperbolic tangent function Tanh(x), expressed as ex-e-x ex+e-x . Then we can construct new activation functions σ(x) = Sig(x + 1) - Sig(x) and σ(x) = T anh(x + 1) - T anh(x), which satisfy that

for 0 ≤ k ≤ 2, where C is a universal constant.

With this exponential decay condition, we can deduce from the proof of Theorem 1 that for 3 < s < 6,

Besides these activation functions, Theorem 3 in [Siegel and Xu, 2020b] implies that, when σ ∈ W m,∞(R) is a nonconstant periodic function, f has a similar integral representation for σ. Moreover, in this case, we do not need to truncate b, which may yields a better approximation rate.

Proof Sketch: In the integral representations of functions within Bs (Ω), the integral can be partitioned into two components, corresponding to bounded and unbounded regions, respectively. Specifically,

where the domain A := {∥ω∥1 ≤ M, |b| ≤ 2M} and we write ˆ f for ˆ fe for simplicity. For the first part f1(x), it can be approximated by a random neural network. The second part, f2(x), can be further decomposed into two terms: one corresponding to the integral over {∥ω∥1 > M, b ∈ R} and another to the integral over {∥ω∥1 ≤ M, |b| > 2M}. We denote these parts as f21(x) and f22(x), respectively. Given that the activation function σ is compactly supported in [-2, 2], we can deduce that σ(ω · x + b) = 0 in {∥ω∥1 ≤ M, |b| > 2M}, since

which implies f22(x) = 0.

Therefore, to ensure that f(x) can be well-approximated well by the random neural network designed to approximate f1(x), it suffices to guarantee that f21(x) has a small enough H2 norm. It is for this reason that we impose the polynomial decay condition on σ. With this condition, we can establish that

which is crucial to achieve our goal. The full proof can be found in the appendix of [Xu et al., 2024].

Remark 3. In [Gonon, 2023], the authors also considered employing a random neural network to approximate the function f belonging to Bs (Ω), with Wi having a strictly positive Lebesgue-density πw on Rd and Bi having a strictly positive Lebesgue-density πb on R. Despite achieving a dimensionindependent approximation rate of O(1/√m), the densities πw and πb are dependent on the unknown decay of the Fourier transform of f, leaving the prior distributions for Wi and Bi unclear. For general s, there may not be densities πw and πb that ensure the constants in the approximation rate are finite. Nonetheless, the approximation results in [Gonon, 2023] remain valid, because the target function they consider possesses a Fourier transform that decays exponentially. Moreover, the method employed in [Gonon, 2023] is specifically tailored for the ReLU activation function and the approximation results are given in the L∞ norm, which renders it inapplicable to the context of PINNs.

For the L2 regression problems, as shown in [Braun et al., 2021], random neural networks can achieve a rate of O(1/√n) for the final prediction error, where n is the sample size. This study assumes that the regression function f satisfies the condition

for all ω ∈ Rd with ∥ω∥2 ≥ 2 and some constant c1 > 0. It is evident that this condition implies that

Thus, the assumption made in [Braun et al., 2021] is stronger than the one in our study and may be more challenging to verify.

Clearly, the approximation rates presented in Theorem 1 are affected by the curse of dimensionality, which is attributed to the unknown decay of the Fourier transform of the target function. When the smoothness parameter s is sufficiently large, we can choose not to truncate ω and b, which may lead to better results. Specifically, assume that s > α+β+5 with α > d, β > 1, then we can derive the following representation for f ∈ Bs (Ω).

where P1(ω) = p1(ω)dω, P2(b) = p2(b)db are probability measures with density functions p1(ω) and p2(b) in Rd and R, respectively, defined as p1(ω) := Cα (1+∥ω∥2)α and p2(b) := Cβ (1+|b|)β .

Building upon our previous discussion, we now consider a scenario where the target function f belongs to a smaller function space. Specifically, we assume that f satisfies the following condition.

Condition 1. Given function f, there is a function fe ∈ L1 (Rd ) such that fe|Ω = f and

with s > α + β + 5 and α > d + 4, 1 < β < 5.

Given the integral representation of the function f that satisfies Condition 1, i.e., (17), we can derive the following approximation results for f.

Theorem2. LetW1, · · · ,Wm ∼ P1 andB1, · · · ,Bm ∼ P2, then there exists a Rm-valued, σ(W,B)-measurable vector A = (A1, · · · ,Am) such that

where

and

for 1 ≤ i ≤ m.

Remark 4. Note that 1 Cα ∼ 1 α-1 Hd-1 (Sd-1 ), i.e., the hypervolume of the (d - 1)-dimensional unit sphere. Therefore, when d is sufficiently large, this term becomes small, and dβ+1 Cα does too.

Remark 5. Here, we aim to discuss the relationship between the class of functions we are considering and the Sobolev class of functions. Lemma 2.1 in [Xu, 2020] shows that for function u and any ϵ > 0,

This demonstrates that Sobolev spaces of sufficiently high order can be embedded into Barron spaces. Furthermore, the form of the definition for the functions considered in Condition 1 can be viewed as an equivalent definition of functions in Hs/2 (Ω). Therefore, overall, these two cases are not very strict assumptions.

3.2 Optimization and Generalization Results

Within the framework of Physics-Informed Neural Networks (PINNs), we begin with the following elliptic partial differential equation, endowed with Dirichlet boundary conditions.

where V, f, g ∈ L∞(Ω).

Let u∗ be the solution of the given equation (21), we consider two scenarios:

(1) u∗ belongs to the Barron space Bs (Ω);

(2) u∗ satisfies Condition 1.

In the first scenario, for brevity and readability, we focus on the situation where s is sufficiently large. Specifically, we take s such that 2(s - 3) ≥ d + 1, i.e., s ≥ d+7 2 . For other values of s, as can be seen from the proof in the appendix of [Xu et al., 2024], the process does not undergo significant changes; the only adjustment required is in the final step concerning the selection of values for m and λ.

According to Theorem 1, we can formulate a random neural network u given by the expression:

where ω1, · · · , ωm ∼ P1, b1, · · · , bm ∼ P2 (as defined in Theorem1) and a = (a1, · · · , am) ∈ Rm is to be determined. For 1 ≤ i ≤ m, we denote the k-th component of ωi by ωi,k.

Given the form of the random neural network u, the loss function of PINNs can be written as

where g1(x) = -f(x), g2(y) = -g(y) and the i-th components of Rm-valued functions f1(x) and f2(y) are defined as

and f i2 (y) = σ(ωi · y + bi), respectively.

To solve for a, we introduce a regularization term to the loss function L(·), resulting in the following objective function:

where λ is a hyperparameter that controls the regularization strength and is to be determined later.

For simplicity, assuming an equal number of samples from the interior and the boundary, the empirical version of F(a) can be expressed as

where

In the first case, according to Theorem 1, there exists a random neural network um = m∑ i=1 ¯aiσ(ωi · x+ bi) such that

where

for 1 ≤ i ≤ m.

Note that ∥¯a∥2 = √∑mi=1 ¯a2i = O(1), thus we can employ the projected gradient descent method for

where C is a constant ensuring that ¯a ∈ A. Specifically, the method of projected gradient descent consists of iteration of the following update rule for t = 1, · · · , T,

where T is the total number of iterations, a1 is an initial approximation and PA is the projection operator onto the convex closed set A. The step size νt is a positive scalar chosen according to a specific schedule or rule.

From the specific form of Fn(·), we can see that it is strongly convex and smooth on A. In convex optimization, it is well-known that for a function f that is α-strongly convex and β-smooth on A, the projected gradient descent with νt = 1/β ensures that for all non-negative integers t, the following inequality holds for the sequence {yt} generated by (31):

where κ = β/α.

The strong convexity of Fn(·) is crucial for the optimization process, ensuring that the algorithm converges to the optimal solution efficiently. Moreover, this property contributes to the generalization performance of the random neural networks, which can achieve the fast rate O(1/n).

Theorem 3. Assume that u∗ ∈ Bs (Ω) with s ≥ d+7 2 , by applying the projected gradient descent in (31) with νt = O(m- d+3 d+1 ) to the empirical loss function Fn(·), we have that

EW,BEX,Y [L(aT ; f1, f2, g1, g2)]

where

and aT is obtained from (31) after T iterations.

Taking λ = 1 √m , m = n d+1 3d+7 and T = O(√n log(n)), we have

The first term represents the approximation error. The second term arises from the addition of the regularization term, which imposes strong convexity to the empirical loss function. This modification allows for the application of the projected gradient descent algorithm, ensuring rapid convergence. The third term corresponds to the optimization error, while the fourth term is the generalization error. The regularization term also contributes to the generalization error, leading to the fast rate O(1/n), which is attributed to the strong convexity of the empirical loss function and the utilization of the localization technique.

Remark 6. In Theorem 3, we have not show the dependency on the dimension, which is an important factor. In fact, from the proof, the final bound provided in Theorem 1 depends quadratically on |Ω| and |∂Ω|. Furthermore, the bound is linear with respect to the constant involved in the Sobolev trace theorem, which is utilized for estimating the approximation error. The generalization error can also be quantified by an appropriate Sobolev norm. Specifically, as demonstrated in [Agmon et al., 1959], we have

where u∗ is the solution of equation (21).

Remark 7. In [Klochkov and Zhivotovskiy, 2021], algorithm stability has been utilized to derive the excess risk for the projected gradient descent method applied to a loss function that is λ-strongly convex and L-Lipschitz continuous. It has been demonstrated that T = O(n2 ) iterations are required to attain the excess risk ˜O( L2 λn ). It is important to note that our approach to deriving the generalization bound relies solely on the strong convexity and Lipschitz continuity of the loss function. Under the same conditions as in [Klochkov and Zhivotovskiy, 2021], by employing the projected gradient descent method with νt = 2 λ(t+1) and aT = 2 T (T +1) ∑Tt=1 tat, our method can achieve the same excess risk within T = O(n) steps. Additionally, because the loss function in our setting is also smooth, T = O(√n log(n)) is sufficient to reach the same excess risk.

In the second case, because the randomly sampled hidden weights are no longer bounded, the loss function can become unbounded. Consequently, projected gradient descent is not suitable for this optimization problem, as the output weights of the random neural networks need to be bounded to achieve the approximation results stated in Theorem 2. Moreover, standard gradient descent is also unsuitable. Although the empirical loss function is convex, it lacks Lipschitz continuity, which means there is no guarantee of convergence for gradient descent. Instead, we can directly apply the closed form of ridge regression. According to the following theorem, we can choose m = n1/4 , which provides an acceptable level of computational complexity for employing closed-form solutions.

Theorem 4. Assume that u∗ satisfies Condition 1, and let an = argmina∈Rm Fn(a), then we can deduce that

where

By taking λ = O(1), m = n1/4 , we have

Moreover, an can be solved by using the closed form of the ridge regression and the complexity is O(m2n) = O(n 32 ).

Remark 8. The error components are distinctly categorized: the first represents the approximation error, the second is attributed to regularization, and the third captures the generalization error. In Theorem 4, an optimization error is notably absent, contrasting with Theorem 3, due to the implementation of the closed-form solution from ridge regression. Furthermore, the unbounded nature of the loss function makes it impossible to to obtain a fast rate for generalization error.

4 Experimental Results

4.1 Verification of Theorem

Theorem 4 shows if m = n 14 , where m is the hidden width and n is the collocation points number, then the testing loss is expected to be O( 1m ). We demonstrate this result from a simple 1D Poisson equation.

The hidden layer's width ranges from 50 to 500, with tanh activations. In Figure 1, we can see that the estimated physics-informed loss decays as the width m increases. The decay trend is consistent with the theoretical decay O( 1m ), and the fluctuations may come from the different random initialization of the first layer.

4.2 Application

We apply the random neural network to solve the ReactionDiffusion (RD) equation, which is significant to nonlinear physics, chemistry, and developmental biology. We consider the spatial-temporal Reaction-Diffusion equation with the following form:

where d1 = d2 = 0.01.

We choose Nb = 300 randomly sampled points on the initial domain and the boundary domain, and Nf = 5, 000 randomly sampled points within the domain Ω = [0, 1]×[-1, 1]. A random neural network with 4 hidden layers, each containing 128 units with tanh activation functions, is used for all computations. The training process runs for 10, 000 epochs with a learning rate of lr = 1e - 3. The relative L2 error is 0.13%. As shown in Figure 2, the random neural network method achieves high accuracy. The training time for conventional PINN method is 1397s in [Li et al., 2024], while the training time for random neural network is 3.2s. This highlights the high efficiency and precision of our approach.

5 Conclusion and Discussion

In this paper, we have established approximation results for Barron-type functions using random neural networks. Our approach diverges from prior studies by providing precise forms for the sampling distributions of hidden weights, rather than relying on empirical selection. Within the PINNs framework, applying random neural networks to solve PDEs translates into formulating regularized least square problems. We address these problems by employing projected gradient descent for one case and leveraging the closed-form solutions for regularized least squares in the other. Combining with the approximation results, we present a comprehensive error analysis. Nevertheless, this work has its limitations. For example, the approximation rate presented in Theorem 1 for functions in Bs (Ω) suffers the curse of dimensionality. Moreover, both approximation results in Theorem 1 and Theorem 2 are only valid for relatively small function spaces. Expanding the approximation capabilities of random neural networks to more general functions presents an interesting direction for future research. Furthermore, the methods we employed for deriving the optimization and generalization errors could be extended to other methods involving random neural networks for solving PDEs, such as the Deep Ritz Method.

Acknowledgments

This work was partially supported by the National Natural Science Foundation of China (No. 12025104, No. 62106103, ), and the basic research project (ILF240021A24). This work is also partially supported by High Performance Computing Platform of Nanjing University of Aeronautics and Astronautics.