{'fraction_non_alphanumeric': 0.04225526131543182, 'fraction_numerical': 0.012725999752893208, 'mean_word_length': 4.408017817371937, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Momentum Diminishes the Effect of Spectral Bias in Physics-Informed Neural Networks

Physics-informed neural network (PINN) algorithms have shown promising results in solving a wide range of problems involving partial differential equations (PDEs). However, they often fail to converge to desirable solutions when the target function contains high-frequency features, due to a phenomenon known as spectral bias. In the present work, we exploit neural tangent kernels (NTKs) to investigate the training dynamics of PINNs evolving under stochastic gradient descent with momentum (SGDM). This demonstrates SGDM significantly reduces the effect of spectral bias. We have also examined why training a model via the Adam optimizer can accelerate the convergence while reducing the spectral bias. Moreover, our numerical experiments have confirmed that wide-enough networks using SGDM still converge to desirable solutions, even in the presence of high-frequency features. In fact, we show that the width of a network plays a critical role in convergence.


Introduction
Physics-informed neural networks (PINNs) have been proposed as alternatives to traditional numerical partial differential equations (PDEs) solvers [1,2,3,4]. In PINNs, a PDE which describes the physical domain knowledge of a problem is added as a regularization term to an empirical loss function. Although PINNs have shown remarkable performance in solving a wide range of problems in science and engineering [5,6,7,8,9], they often fail to accurately construct solutions when the target function contains "high-frequency" features [10,11]. This phenomenon exists in even the simplest linear PDEs. Since PINNs are relatively new, only a few studies have explored the potential causes of the convergence problems in these algorithms. Here, we briefly explain their approaches.
It has been shown that PINNs normally fail to converge to the solutions when the target functions have "highfrequency" features [10,11]. The issue in the training of the "high-frequency" features is not limited to PINNs. Rahman et. al [12] empirically and Cao et. al [13] theoretically showed all fully-connected feed-forward neural networks (NNs) exhibit spectral bias, meaning NNs are biased to learn "low-frequency" features that correspond to less complex functions [12]. To better understand the spectral bias of PINNs, by assuming an infinitely wide network, Wang et. al [14] employed neural tangent kernel (NTK) theory [15] to analyze the training convergence of the networks. By examining the gradient flow of PINNs during the training, they showed that under the stochastic gradient descent (SGD) optimization process a very strong spectral bias was exhibited. Moreover, Wang et. al [11] concluded that in the presence of high-frequency features the convergence to the solution under SGD may completely fail. They proposed to assign a weight to each term of the loss function and dynamically update it. Although the results showed some improvements, this approach could not completely minimize the effect of spectral bias. For target functions with higher frequency terms, their PINNs still failed to converge to the solution. Moreover, as assigning weights may result in indefinite kernels, the training process could become extremely unstable.
Krishnapriyan et. al [10] took a different approach. Instead of modifying the loss function terms, they implemented a curriculum learning method. In their approach, the model was gradually trained on target functions with lower frequencies, and, at each step, the optimized weights were used as the warm initialization for higher-frequency target functions. Their method showed good performance on some PDEs. However, as the frequency terms become larger, more steps are required and the step sizes between one frequency to the next one are not known. Thus, several trials and errors are required to find the optimized initialization terms.
In the present work, we prove that an infinitely-wide PINN under the vanilla SGD optimization process will converge to the solution. However, for high-frequency terms the learning rate needs to become very small. This makes the convergence extremely slow, hence in practice not possible. Moreover, we prove that for infinitely-wide networks, using the SGD with momentum (SGDM) optimizer can reduce the effect of spectral bias in the networks, while significantly outperforming vanilla SGD. We also discuss why the Adam optimizer can accelerate the optimization process relative to SGD. Although the concept of infinitely-wide PINNs is not new, to our knowledge this is the first time that the gradient flow of the output of PINNs under the SGDM is being analyzed, and its relation to solving spectral bias is provided. Moreover, until now, infinitely-wide NNs were only used theoretically to facilitate the analysis of the dynamics of the learning process. In this work, by using the numerical experiments, we show that when dealing with high-frequency features, for fixed-depth networks, the training error decreases as the networks become wider. Hence, depending on the problem in hand, a sufficient width may be reached such that, for practical use, the training error becomes sufficiently small.

PINNs General Form
The general form of a well-posed PDE on a bounded domain (˙⊂ R d ) is defined as: where D is a differential operator and u(x) is the solution (of the PDE), in which x = (x 1 , x 2 , . . . , x d ). Note that for time-dependent equations, t = (t 1 , t 2 , . . . , t d ) are viewed as additional coordinates within x. Hence, the initial condition is viewed as a special type of Dirichlet boundary condition and included in the second term.
Using PINNs, the solution of Eq. 1 can be approximated as u(x, w). That is the approximation resulting from minimizing the following loss function: i =1 are boundary and collocation points respectively, and w describes the neural network weights. The first term of Eq. 2, which we will denote as L b (w), corresponds to the mean squared error of the boundary (and initial) condition data points. The second term, denoted L r (w), encapsulates the physics of the problem using the randomly selected collocation points. Similar to all other NNs, minimizing the loss function L (w) results in finding the optimized solution for u(x, w).

Infinitely-Wide Neural Networks
The output of a fully-connected infinitely-wide NN with L hidden layers can be written as [15]: whereˆh and b h are respectively the weight matrices and the bias vectors in the layers h = 1, . . . , L, N is the width of the layer, x is the input vector, and σ(·) is a β-smooth activation function (such as tanh(·)). The final output of the NN is written as: At each time step t we can determine the change of the output with respect to the input, which defines a neural tangent kernel (NTK): One major consequence of dealing with infinitely-wide NNs is that if the last layer of the network is linear the NTK becomes static (does not change in time) [16]. Another major consequence is that the norm of the Hessian of the output becomes smaller as the width become larger. In fact, as N → ∞ the norm of the Hessian becomes 0.

Convergence of the Gradient Descent Algorithm in the Presence of High-Frequency Features
Generally, the optimization problems corresponding to over-parametrized systems, even on a local scale, are non-convex. A loss function L (w, x) of a µ-uniformly conditioned NN satisfies the µ-PL * condition on a set S ⊂ R m if: where µ is the lower bound of the tangent kernel K (w) of the NN. It has been shown that infinitely-wide networks satisfy the µ-PL * condition, a variant of the Polyak-Lojasiewicz condition, and as a result the (stochastic) gradient descent (SGD/GD) optimization algorithms will converge to the optimal solution [17]. The following proposition makes use of the µ-PL * condition to provide a convergence analysis for an infinitely-wide PINN with a loss function as in Eq. 2 optimized with SGD (Appendix A).

Proposition 1.
Let λ b max and λ r max , respectively, be the largest eigenvalues of the Hessians ∇ 2 L b (w t ) and ∇ 2 L r (w t ).
Consider an infinitely-wide PINN optimized with the following update rule: where η is a constant learning rate. Then, provided η = O (1/(λ b max + λ r max )), the µ-PL * condition is satisfied and the PINN will converge.
Unfortunately, for PDEs exhibiting stiff dynamics (those with large parameters) the eigenvalues (of the Hessian) dictating convergence are often very large. For example, take the one-dimensional Poisson equation: . Thus as C grows so does the bound on the eigenvalues of H (w t ), indicating that at least one of λ b max or λ r max will be large. This makes the learning rate necessary in order to guarantee the convergence of the PINN prohibitively small, and therefore SGD is unusable in practice.

Wide PINN Optimization using SGD with Momentum
In order to accelerate learning, we go beyond the vanilla gradient descent and investigate the training dynamics of PINNs under SGD with momentum (SGDM) [18]: where α and η are fixed learning rates.
Specifically, we analyze the gradient flow of the update rule in Eq. 5 by leveraging the notion of neural tangent kernels (NTK). The following theorem (proof in Appendix C) reveals that SGDM has a different convergence behavior from SGD:

Theorem 1. For an infinitely-wide PINN, the gradient flow of SGDM is:
where µ is related to momentum by: α = m m+µ´t and K is: and: are three NTKs associated to the boundary and residual terms. Moreover, let γ = µ/2m, κ i be the eigenvalues of K, and Then, the solutions to the gradient flow are of the form: where A 1 and A 2 are constants.
Examining the analogous gradient flow for the vanilla SGD, the training error decays at the rate of e −κ i t [11]. That is, the eigenvalues of the kernel indicate how fast the training error of SGD decreases. Spectral bias arises as eigenvectors of larger eigenvalues of the kernel matrix represent lower frequencies [12,19]. Thus, it is clear that under vanilla SGD the high-frequency components are extremely slow to converge.
Once we add momentum, the decay rate analysis becomes more involved as Eq. 6 yields three different cases of solutions. Each of the three cases are analogous to one of the solutions of a damped harmonic oscillator in physics [20,21] (Fig. 1). • Under-damped: Imaginary roots (γ 2 < κ i ) • Critically-damped: Real and equal roots (γ 2 = κ i ) • Over-damped: Real roots (γ 2 > κ i ).
Under-damped case As γ 2 − κ i has imaginary roots, Eq. 6 can be rewritten as: where ω 1 = κ i − γ 2 , and A and φ are two constants corresponding to the amplitude and the phase of the damped oscillation. In physics, this solution corresponds to an oscillatory motion in which the amplitude is decaying exponentially.

Critically-damped case
The general solution of the critically-damped case can be written as: where B 1 and B 2 are constants. As the oscillation motion is not present, the decaying rate for this case is much faster (Fig. 1).
Over-damped case Lastly, in the over-damped case the general solution is simplified as: where ω 2 = γ 2 − κ i and C 1 and C 2 are constants. Similar to the critically-damped case, the above equation states a fast decay.
Thus, depending on the absolute value of the eigenvalues of the kernel matrix (|λ i |), the dynamics of the training error can differ. For larger eigenvalues the training error corresponds to an under-damped solution, in which the amplitude of an oscillatory motion is decaying exponentially, whereas for smaller eigenvalues the correspondence is with an over-damped or critically-damped oscillation, with much faster decay. Thus, when using SGDM instead of vanilla SGD, the training process for high-frequency components of the target (corresponding to the smaller eigenvalues) will decay faster (undergoing the over-damped or critically-damped motions) while the low-frequency components (corresponding to the larger eigenvalues) will decay at a slower rate. As a consequence, the effect of spectral bias will be less prominent (compared to vanilla SGD). Of note, the fastest decay is observed for the critically-damped case, in practice many small eigenvalues are near the value of γ and the dynamics of the training process becomes very close to the critically-damped case.

Analysis of Adam for Band-Limited Functions
Adam is a commonly-used optimizer which can be interpreted as SGDM adapted for variance [22,23]. Optimizing with any gradient-descent-based method, such as Adam, we sample a batch X of N data points, and compute the stochastic gradient g (w; X ) = x∈X ∇ (w; x)/N (where in our case is the respective loss for whether x is a boundary or collocation point). Define the pointwise variance σ = var g (w; X ), and denote by σ i and ∇L i their i -th entries (henceforth we stop writing w for clarity). At each epoch, Adam then updates the i -th entry of w with a magnitude inversely proportional to an estimator of σ 2 i /∇L 2 i [23]. The end-goal of the training process is to find a w which minimizes the expected loss across all collocation and boundary points L, approximated by L in Eq. 2.
Our numerical experiments confirm that for wide NNs both SGDM and Adam can converge to desirable solutions, but that Adam is even faster than SGDM (see Section 4). This may be due to variance adaptation alleviating spectral bias even further. The connection is readily seen when the solution to our PDE u(x) is (at least well-approximated by) a band-limited function, meaning u(x) = k∈K α k exp(2πi kx) for some finite set of frequencies K .
For such solutions we eventually satisfy the bound: with high probability [24,25]. This indicates the loss is only guaranteed to be small if the k ∈ K are themselves small, precisely when u(x) consists of low frequencies. In other words, L can be very large for high-frequency solutions. Since our NN is wide and thus satisfies the µ-PL * condition ∇L 2 ≥ µL of Eq. 3, such solutions imply ∇L and thus the ∇L i also have loose bounds.
Convergence happens when ∇L i ≈ 0 for all i , so if any ∇L i 0 the NN is far away from the solution. However, this is precisely the situation in which Adam takes large steps and accelerates toward the solution fastest, as σ 2 i /∇L 2 i is close to zero. This phenomenon certainly gives Adam an edge over SGD, and may contribute to it outperforming SGDM too.

Numerical Experiments
In this section, using Poisson's equation, transport function, and the reaction-diffusion problem we provide results from our numerical experiments.

Poisson's Equation
Poisson's equation is a well-known elliptic PDE in physics. For example, in electromagnetism, the solution to Poisson's equation is the potential field of a given electric charge. In Eq. 4 the general form of the one-dimensional Poisson's equation was presented. Here, we write the Poisson's equation for a specific source function and a specific boundary condition: with u(x) = sin(C x) as the solution. The corresponding loss function is written as: whereû is the output of the network. We chose N r = 1000 and N b = 1000 and built a simple NN with 1 hidden layer of width 500. As expected, when C = 5π, after only 45000 epochs, trained models via all three algorithms could accurately estimate the solution. The relative error (u −û)/û of the output of NN using the vanilla SGD was on order of 10 −2 , indicating that the model could successfully predict the solution for the PDE. The relative error using SGDM and Adam were respectively on order of 10 −3 and 10 As the parameters of a PDE grow larger, converging to the solution becomes more challenging. For Poisson's equation when C = 15π, after 55000 epochs, the obtained solution from training the network with the vanilla SGD exhibited a relative error on order of 10 −1 . Meanwhile the trained model with SGDM, after being trained for 55000 epochs, provided an acceptable solution with a relative error on order of 10 −2 . Moreover, Adam converged to the solution after training the model for only 25000 epochs. The relative error was also considerably smaller and on order of 10 −3 . The comparison between the results of training the model using the three algorithms is shown in Fig. 2. Of note, using vanilla SGD and training the model for 180000 epochs resulted in a small relative error on the order of 10 −2 (Fig. C.2). This confirmed our discussion presented in Section 3.1 that implementing vanilla SGD will result in convergence under large parameters, though due to the presence of high-frequency features, the convergence becomes extremely slow.

Transport Equation
The transport equation is a hyperbolic PDE that models the concentration of a substance flowing in a fluid. Here, we focus on a one-dimensional transport equation: is the initial condition, and β is a constant parameter. Using the methods of characteristics, the transport function has a well-defined analytical solution: u(x, t ) = g (x − βt ). The corresponding loss function is written as: 2 whereû is the output of the network. We used N b = 450 and N r = 4500. We also chose the boundary and initial conditions to be: We

Reaction-Diffusion Equation
A reaction-diffusion equation contains a reaction and a diffusion term. Its general form is written as: where u(x, t) is the solution describing the density/concentration of a substance, Δ is the Laplace operator, ν is a diffusion coefficient, and f (u) is a smooth function describing processes that change the present state of u (for example, birth, death or a chemical reaction). Here, we assume a one-dimensional equation, where f (u) = ρu (1 − u). For a constant value of ρ and with the defined f (u), Equ. 4.3, can be solved analytically [10]. To estimate the solution, using PINNs, the loss function for the 1D reaction-diffusion PDE is written as: We further assigned N r = 1000, and N n = 100. We also chose the initial and boundary conditions as u(x, 0) = exp{− (x−π) 2

Width of Neural Networks
An interesting observation was that for PDEs with large parameters the width of the neural networks played an important role. For example, in Poisson's equation with C = 20π, we used an 8-layer NN with a width of 50 neurons (per layer), and trained the model via Adam optimizer. The final output had a large relative error of order of 10 0 , while a shallow network with a width of 500 neurons could converge to the solution and had small error on order of 10 −3 (Fig. 5(a) and (b)). Moreover, for a fixed depth of 2, we trained several networks at different widths. As the networks became wider the relative error decreased (Fig. 5(c)). Also, in the transport function, for β = 100, a 4-layer network with a width of 50 had a large relative error on order of 10 −1 . However, another 4-layer network with a width of 500 neurons had a considerably smaller relative error on order of 10 −2 . Similarly, in the reaction-diffusion equation for ν = 23 and ρ = 25, a 4-layer network with a width of 50 had a large relative error on order of 10 −1 , while training another 4-layer network with a width of 500 resulted in a much smaller error on order of 10 −3 (Fig. 6).

Conclusion
In the present study, through the lens of NTKs, we examined the dynamics of training PINNs via SGDM. We also showed that under SGDM the conversion rate of low-frequency features becomes slower (analogous to the underdamped oscillation) and many low-frequency features undergo much faster conversion dynamics (analogous  to over-damped and critically-damped oscillations). Thus, the effect of spectral bias becomes less prominent. Moreover, we discussed how training a PINN via Adam can also accelerate the convergence, and showed it can be much faster than SGDM. Although we analyzed the dynamics of convergence by assuming an extreme case of infinitely-wide networks, our experiments confirmed that even when the width of the networks are not infinite, the estimated solutions obtained from the trained models via SGDM and Adam had high accuracy. Moreover, we observed that when dealing with target functions with high-frequency features, wider networks (compared to their narrower counterparts) could estimate solutions with much smaller errors. The fact that the width of the networks significantly enhanced the performance of the models indicates that the neural architecture may play an important role in PINNs. Thus, we are planning to explore more neural architecture choices. Moreover, one major future step is to examine if PINNs can solve more complicated PDEs. For example, researching if they are capable of solving highly non-linear equations, and if spectral bias becomes a more serious problem for those PDEs.
where H (w t ) = ∇ 2 L b (w t ) + ∇ 2 L r (w t ). Following the approach of [11], let Q be an orthogonal matrix diagonalizing H (w t ) and v = L (w t )/ L (w t ) . With y = Qv, we have the following: where the λ b i and λ r i are the respective eigenvalues of the Hessians of L b and L r ordered non-decreasingly, and the summation is taken over the components of y.
From here, fixing any bound B such that λ b max + λ r max ≤ B, we obtain the following inequality: Witness that when η = 1/B, we can further simplify this to: Since our NN satisfies the µ-PL * condition from Eq. 3 due to its width, we therefore have:

B Hessian for Poisson's Equation
The Hessian of the PDE part of the loss function is calculated following the methods introduced in [11], where the gradient of loss function is: Here, u(x) is the target solution admitting some parameter C , and u w (x) is the NN approximation of the output. Assuming that the approximation is a good representation of the actual solution, it can be written as u w (x) = u(x) w (x), where w (x) is a smooth function taking values in [0, 1], such that | w (x) − 1|< δ for some δ > 0 and ∂x k | ≤ δ. The Hessian of the loss function will therefore be: where I = ∂L r ∂w , and can be calculated to be: The above equation contains 3 terms, which we call I 1 , I 2 , and I 3 respectively, and define ∂I ∂w = ∂I 1 ∂w + ∂I 2 ∂w + ∂I 3 ∂w . We calculate these terms as follows:

F Comparison
Here, we compare our numerical solution of the transport function, the reaction-diffusion problem and the Poisson's equation based on training the model using the Adam optimizer for a wide network, with the solutions presented in [10] and [14]. In the transport function, β = 40 was the largest parameter that was analysed in [10], and in the reaction-diffusion problem, ν = 5 and ρ = 6 were the largest shown parameters. Similarly, the solution of the Poisson's equation in [14] was shown when C = 4π. We used the same parameters to show the comparison between our methodology and the mentions studies. In Table 1, the relative and absolute errors of the estimated solutions (based on our approach and method presented in [10] and [14]) are presented. Our solution for the transport function when β = 40 is as good as the presented solution in [10] (the estimated solution from both methods have errors on order of 10 −2 ). However, in the case of reaction-diffusion and the Poisson's equations our estimated solutions are more accurate and have smaller relative and absolute errors. Of note, we have already shown that for much larger choices of parameters (see Section 4.4) our trained models are also capable of estimating accurate solutions.
To build our networks, we used the same number of hidden layers as used in [10] and [14]. However, we made our networks much wider. Moreover, for Poisson's equation we trained our model, using the same number of epochs as proposed in [14]. The comparison between the choices of basic hyperparemters such as number of hidden layers, width of each layer, and number of epochs are shown in Table 2 [10] and [14] for 1D transport function, the reaction-diffusion equation, and Poisson's equation are shown in the third column. The relative and absolute errors of the solutions presented using our methodology are also presented in the fourth column.
Neural Network Parameters Network of [10] for Transport Function  The neural network parameters such as number of hidden layers, width of each layer, and number of epochs are compared between our model and the presented models in [10] and [14].