Physics informed neural networks for elliptic equations with oscillatory differential operators
Authors: Arnav Gangal
Luis Kim
Sean Patrick Carney
Venue: Accepted by TMLR ()

Abstract: Physics informed neural network (PINN) based solution methods for differential equations have recently shown success in a variety of scientific computing applications. Several authors have reported difficulties, however, when using PINNs to solve equations with multiscale features. The objective of the present work is to illustrate and explain the difficulty of using standard PINNs for the particular case of divergence-form elliptic partial differential equations (PDEs) with oscillatory coefficients present in the differential operator. We show that if the coefficient in the elliptic operator $a^{\epsilon}(x)$ is of the form $a(x/\epsilon)$ for a 1-periodic coercive function $a(\cdot)$, then the Frobenius norm of the neural tangent kernel (NTK) matrix associated to the loss function grows as $1/\epsilon^2$. This implies that as the separation of scales in the problem increases, training the neural network with gradient descent based methods to achieve an accurate approximation of the solution to the PDE becomes increasingly difficult. Numerical examples illustrate the stiffness of the optimization problem.

Published in Transactions on Machine Learning Research (10/2023)
Physics informed neural networks for elliptic equations with oscillatory differential operators
Arnav Gangal
agengal00@@mail.com com
Deportment of Mathematics University of California, Los Angeles
Luis Kim
lnixyookin022@@mmill.coo
Department Mathematics
The University of Texas at Austin
Sean P. Carney
spcorney@@aath ucla edu
Department of Mathematics University of California, Los Angeles
Reviewed on OpenReview: https: // openreview. net/ forum? id- Q.y
Abstract
Physics informed neural network (PINN) based solution methods for differential equations have recently shown success in variety of scientific computing applications. Several authors have reported difficulties, however, when using PINNs to solve equations with multiscale features. The objective of the present work is to illustrate and explain the difficulty of using standard PINNs for the particular case of divergence-form elliptic partial differential equations (PDEs) with oscillatory coefficients present in the differential operator. We show that if the coefficient the elliptic operator a" (x:) is of the form a(x/c) for a 1-periodic coercive function a(-), then the Frobenius norm of the neural tangent kernel (NTK) matrix associated to the loss function grows 1/c². This implies that the separation of scales in the problem increases, training the neural network with gradient descent based methods to achieve an accurate approximation of the solution to the PDE becomes increasingly difficult. Numerical examples illustrate the stiffness of the optimization problem.
1 Introduction
Recent developments in deep learning have shown great promise for advancing computational and applied mathematics 2021). Physics informed neural networks (PINNs) have recently emerged method for scientific computation. Building off carlier work Lagaris et al. (9998) Psichogios & Ungar (1992). PINNs were introduced in Raissi et al. (2019) and seek to approximate the true solution of a differential equation by a neural network u(x; 8) parameterized by weights and biases 0. Consider the following general partial differential equation (PDE) defined on some Ω C R by differential operator /
popular
[1) f(x).
I Ω
u(x) g(x),
E aΩ.
(1)
where x and Q represent some general space-time coordinate and domain, respectively. Here and throughout this paper, the boundary conditions are taken to be of Dirichlet type. By the uniform approximation theorem (Horiik et al. 1989 Leshno et. al. 1993), given sufficient. data, a neural network can uniformly approximate
classical smooth solutions to (1) whenever they exist.
In practice, the network parameters are determined by minimizing a loss function that enforces (
hold for set of N. collocation and N₂, boundary points, N. E Q and (st) 82, respectively. The

---

Published in Transactions on Machine Learning Research (10/2023)
loss is
1 N. No
Nu u(s;; 0) g(si) -g((s)
[()) = Σ Mu](<;;;) f(x;))
(2)
Nb. ==1
where À is a tunable parameter that weighs the relative importance of the boundary conditions. We refer to Karniadakis et al. (2021) for a review of this methodology applied in a wide range of contexts in scientific computing.
Despite the success of PINNs a wide variety of applications, numerous authors have reported difficulties applying the technique to problems with multiscale features. Some examples include a scalar, nonlincar hyperbolic equation from model of (wo-phase immiscible fluid flow in porous media. (Fuks & Tchelepi 2020), which kinetics (Ji ct al. 2021), which exhibit stiff dynamics that evolve over wide range of time scales. Difficulties have also been reported for the Helmholtz equation (Wang el al. 2021a), , as well as linear hyperbolic problems, for example for example the one-dimensional advection equation (Krishnapriyan et al., 2021 and the wave equation (Wang et al. 2021b).
support shock waves, and systems of ordinary differential equations governing chemical
The focus of the current work is to illustrate and explain the difficulty of using standard PINN solution methods for linear elliptic boundary value problems (BVPs) of the form
(z) = f(x), -V
x EEΩ
g(x),
x E an,
where the coefficient tensor a' R" - R" is uniformly bounded and coercive in e and assumed to consist of entries that contain frequencies on the order of c-1 for < c < 1. The measure of the domain Ω C R is assumed to be (1), and hence, (3) a multiscale problem that models, for example, steady-state heat conduction a composite material or porous media flow governed by Darcy's 's law.
Theoretical analysis and numerical experiments presented below illustrate that, whenever standard PINN architectures are used in conjunction with gradient descent based training for the multiscale problem 3 the resulting optimization problem becomes increasingly difficult. as the scale separation in the BVP increases, i.e. as vanishes. After motivating the present study in Section 2 we show in Section 3 that the neural tangent kernel matrix associated with the PINN approximation to 3 has a Frobenius norm that becomes unbounded as 4 O. Numerical examples Section | illustrrate that during training, the ordinary differential equation that governs the evolution of the BVP residuals indeod becomes increasingly stiff as c 4 O, translating poor PINN performance for problems with a large separation of scales.
2 Motivation
The motivation for attempting to use physics informed neural network solutions to the oscillatory problem 3 is to investigate whether connection can be established between asymptotic homogenization theory and the so-called "frequency principle" in deep learning.
Recall from mathematical homogenization theory (Bensoussan et al. 2011 that, under suitable conditions, the solution to is well approximated as by the solution to a homogenized equation
V Vu(x)) = -(x), ((x) = g(x),
x E Ω
x an.
Both the homogenized coefficients 7 and the solution TI do not contain e-scale oscillations; the latter approx- imates the large-scale, low-frequency features of the oscillatory function w'.
Additionally, when neural networks learn a target function, they are known to learn the low frequency components more rapidly than the large frequencics (Rahaman al.] 2019 Xu al. 2020). This "frequency principle" was shown to hold for gradient descent training Luo al. (2021); Markidis (2021): empirically the result be observed for PINNs applied to relatively simple problems (Wang et al. 2021b even when
2

---

Published in Transactions on Machine Learning Research (10/2023)
more widely used optimizers are used for training, e.g. Adam (inggma & Ba 2014 representative example the PINN solution to the one-dimensional Poisson BVP
Consider as brief
d
de u = sin(x) + sin + sin(15r) 15x) + sin((55) =: f(x)
(4)
for x E (=), with homogeneous Dirichlet boundary conditions u(- u(TT) = 0. For collocation points
{ar Jimi E (-T,, T) the loss function becomes
1 N ((0) No
Σ dez
u 0) s(x))) (1 +
(5)
de
Figure 1 shows the evolution (as function of training iteration) of the complex modulus of
6 (ūrre) (ØNN)
(6)
where Utrue is the true solution to the BVP, uNN is the neural network approximation, denotes the discrete Fourier transform, and the index & denotes the frequency of the coefficient. A rolling average is used to make the trajectories more legible. The PINN solution is computed with fully connected neural network with four hidden layers of sixty nodes each, N₂ = 512 equispaced collocation points, and À = =0; see Appendix for a complete description the training process that generates the final neural network output. One can clearly observe the frequency principle this simple example; the low frequencics of the target BVP solution more rapidly learned than the high frequencies. See also Markidis (2021) for another example in two-dimensional Poisson problem, well as Wang al. (2021b) for mathematical analysis of this phenomenon for Poisson equations.
250
200
150
k=1
V
k=5
==
100
k= 55
10
10
10a
10a
10
Training iteration
Figure 1: A simple illustration of the frequency principle: evolution a function of training iteration of the magnitude of the discrete Fourier transform of the error (6) between the true solution to the BVP and the neural network approxination.
Given some evidence that PINNs can learn low frequency features PDE solutions, it reasonable to ask if the homogenized solution 'naturally' ariscs when trying to learn solutions to the multiscale problem ( and, if not, for what reason? The homogenized solution itself is of course useful in a variety of applications; it additionally could be used to construct coarse grid solutions for multigrid solver of the full, oscillatory

---

No Published in Transactions on Machine Learning Research (10/2023)
problem as in Engquist & Luo (1997) The natural smoothing properties of standard iterative solvers, e-g. the method of Gauss-Seidel successive over-relaxation (Briggs et al. 2000), would then rapidly reduce the error of high-frequency solution components that are not easily approximated by the neural networks; such a 'hybrid' multigrid strategy was explored in Markidis
2021 for Poisson equations.
Recent theoretical results Shin et al. (2020) guarantee that as the number collocation N₂ and boundary N₂ points tend to infinity, the sequence of minimizers to a "Hölder regularized" version of the loss function 2 will converge to classical solutions of elliptic and parabolic PDEs (whenever they exist) in the limit of infinite learning data. Nevertheless, the authors in both Han & Lee (0223) and Leung et al. (2022) reported difficulties training neural networks to achieve such minimizers of the loss function associated to the multiscale problem (3). The purpose of the current work is to characterize the reason why optimization is difficult with both theory and numerical examples.
The theoretical results presented below describe the neural tangent kernel matrix of the PINN solution the BVP 3) this matrix appears in the ordinary differential equation that, under gradient descent. training dynamics, governs the evolution of the PDE and boundary residuals that define the PINN loss function. Under few technical assumptions, This result, along with mumerical evidence presented in Section ! suggests that as that as the scale scparation in the differential operator in (3) increases, the optimization problem that determines the neural network approximation to the PDE solution becomes increasingly stiff.
show that the matrix Frobenius norm becomes unbounded as / 0.
3
Neural tangent kernel matrix theory
3.1 Theoretical results
We now describe the neural tangent kernel (NTK) matrix for the physics informed neural network approx- imation to the multiscale elliptic equation before showing that it can grow arbitrarily large in norm + O. See Jacot al. (2018) for the original development of the NTK theory for least-squares regression problems, as well as Wang al. (2221b) and Wang et al. (2022) for extension of the theory to PINNs.
For the presentation below it will be useful to define the linear multiscale differential operator
== -V- (V)
where ip is some suitably regular function. The multiscale elliptic boundary value problem
(3) then becomes
u = g an. L'u =
Ω
(7)
For simplicity, we assumee throughout that the entries of the coefficient tensor a are at least once continuously differentiable, and hence bounded on any compact domain Ω.
Given some neural network u(x; 9) paramcterized by Np weights and biases 0, as well as collocation { E
No Ω and boundary points { E an, the loss function associated to
7 becomes
No.
No
(0)
Σ
Tfu(x;;0) f(xz)
u(s;; 0) g(8;)
(8)
N:
2
No
2
If the residual values
T Tpdo 0) = --u(zx;)-f((am
i I...,,No No
(9)
and
(b(s;; = u(s;; g(ss),
(10)
are grouped together into single vector y(#), and if the parameters A evolve according to the gradient flow
d8
dt = VoL.
(11)
4

---

Published in Transactions on Machine Learning Research (10/2023)
then y(f(t)) will evolve according to initial value problem of the form
d
dt v() K (x((t( (t))
(12)
where K is called the neural tangent kernel matrix. Note here that both 0 and implicitly depend on e; however, the dependence is not explicitly marked.
We next derive 12 noting again that similar results can be found in Wang el al. (2021b) and Wang et al. (2022).
Lemma 3.1. For neural network parameters RR,,, let y(0) Ro be the vector of residual values
v(0) = (Fpde(21; .)..... ( 0), rb(si: ... (s)
where the entries Vdde and Tb are defined by 19 and (10) Suppose that the parameters 0 evolve from some initial value du according to the gradient flow 11 Then y(O(t)) evolves from the initial condition y(fb) according to
d K (y(()) (
di v
(13)
where the explicit dependence on 0 in 13 is dropped for convenience. The neural tangent kernel matrix is given by
(w)(t) (t) (t)
(+)
(14)
Kos(t)
where the subblocks have entries
N₂
(Kin) =
L
L (t))
L' ((t);
i S No.
N
00,
00,
1=1
N₂
X
(t
uss 0(t)
u(s;; 0(t))
15 $ No.
No
a0,
00,
and
(K))) X
a u(s;: 0(t)).
L' xx; f(t))
1 No. 1 1 jisa No.
No
00,
00,
while Kin is simply the scaled transpose of Kūb:
N₂ 1 (Kba)i X No (KW))s
1 No. 15j5 No.
The proof is computation that follows from the chain rule from differential calculus; it can be found in Appendix D We note that while the K and Kí subblocks of the NTK matrix are symmetric, overall, Ke is not symmetric unless N₂- = AN..
Next we show that the NTK matrix 14) associated to the multiscale elliptic problem in particular the Kin.. KG and Kí subblocks can become arbitrarily large for vanishing e. In contrast, because the boundary condition g in (7) is indepedent of c, there is explicit c-dependence in the Kys subblock.
For simplicity we assume the problem is one-dimensional (n 1) and that the multiscale coefficient a*(x) a(z/e) is e-periodic. showing the desired result for 1, we briefly sketch below an additional proof for the higher dimension case n > under additional technical assumptions.
In dimension n =1, 1, becomes
d (a(x/c) d dr dx u (x) = f(x),
E (a,
wf (a) vaa
(b) ub.
(15)
5

---

Published in Transactions on Machine Learning Research (10/2023)
The loss function associated to 15 is equation
with No 2 and elliptic operator L given by
L [= 0) = - (a(x/e) dx²
u(2; 0) + a' (x/c) =z((z))) 0)):
(16)
dr
where it is assumed the oscillatory function differentiable. The 1/c factor the second term will be the source of the divergent behavior of the NTK as + 0.
We first show the result for neural network with one hidden layer before extending it to more general fully connected networks.
Theorem 3.2. Let : R + R be neural network with one hidden layer of width d and smooth activation function a, so that
d u(x; 0) Σ wh a (w x b/o) + +6"
(17)
1=1
where 600) E R, (()) E Ixd and (() 6(0) E Rdx1, and assume that Vo E
0 <') < oc.
(18)
Let aly) be a one-periodic, non-constant C function that is bounded and coercive, so that Vy E R, there exist some X, A E IR such that
0 o<\Saly) SA,
and let (x) = a(x//c). For a fixed number of collocation points No, let K be the neural langent. kernel
matric [4] associated to the loss function (8) of the boundary value problem 115 and suppose the network
paramelers 0 evolve according to the gradient flow (11) for O S T for some T 0.
Suppose also that the porameters 0 evolving according to 11 are bounded uniformly in both t and e:
sup o<cect 0S1ST
sup
|0(t)l <C
(19)
for some co > Finally, let be the of E [0 for which there exists least one entry of (()) that is asymptotically larger than e; that is, te Y and only if
lim e/,, (0) (f) 40
for some 1 SIS d. Then VE r
lim [1)ll <0
where I| the Frobenius norm.
Before proving the result, we brielly remark on the theorem assumptions. First note that the number of collocation points No is assumed to be fixed. Since we will show that the entries of the NTK matrix are proportional to
1
N₂ ?
the singular nature of the evolution the residuals (12) could in principle be prevented by taking No ~ - However, this can make the computational cost of gradient descent based optimization prohibitively expensive whenever the problem scale separation sufficiently large.
Next, note the assumption (18 do do/dr /dr holds for sigmoidal-type activation functions commonly used for PINNs (Lu et al. 2021b), such as the hyberbolic tangent and logistic functions. Although it does not hold for the ReLU function (x) = inax(0,22) these are of course not suitable for PINNs and other NN based solution methods that are based on the strong formulation of PDEs with second order (or higher) differential operators.
Regarding the set T C j0, T), note that if t & Y. then each entry of ()) will converge to in the limit J. O: in this scenario the neural network 17 would limit to constant function and hence not be suitable

---

Published in Transactions on Machine Learning Research (10/2023)
the network parameters A did not hold, then an unstable, divergent training process could result. (Wang et al. 2022;; ; the bound helps prevent the so-called "vanishing gradient" problem that can occur when components of VeL vanish, preventing descent in that direction.
solution 15 in general. Finally, if the uniform bound 19
Proof of Theorem 3 The proof, like that of Lemma 3.1 is a consequence of direct computation; note that it is sufficient for single entry of matrix to diverge to ensure the Frobenius norm diverges. particular we show that entries in the K subblock of the full NTK matrix K scale as 1/(N Indeed, recall from Lemma 3.1 that
1 N Ne Σ L' (x f(t)))
Ô a0,
(K (na
(20)
for any 1 S a S Ne. Since 20 is a sum of squares, if just one entry in the sum blows up as + 0, then the entire sum of course will as well. Below we drop the explicit dependence of the parameters 0 on for notational convenience.
First, let (2; 0) be the Rd valued function whose k-th component equals
(x(x;) σ (w + b²").
so that
u(x; 0) = Σ
( (z; + +0())
k=1
By linearity,
0 L u(:; 0) = 1, g, (x;
aw
for any S S d. So, to prove the desired result, it suffices to show that
1
L' gy 0]] 00
(21)
No
as &. Using
16
the left-hand side of equals
& (a(aa/0) dz29- (Fa:
2
2 (a(taal()
d² (r29)(ra:00) (xa/()) d dx (r(:
Ne
c
+ (d(aaa))) de 9 ())())))))))
d
(22)
Since is periodic C function, both a(ra/c) and d' d'(a/(t are bounded independent of safe E RR. The derivatives of gy are
du 22 0) 0 o'(w:0) + (100 W
d
(23)
and
&
d (x; o" (w:(0> b(0) (W?
(24)
Since the activation function is smooth and the network parameters are uniformly bounded in and 4, both
(23) and (24) are also bounded for all xa E [a, b). Consequently, (22) is dominated by the 1/e² term for vanishingly small. The positivity assumption (18) on do do/dx and the uniform boundedness of the network
parameters imply also that the (absolute value of the) first derivative 23 is bounded below by constant
independent of e; this implies that for any t E T
lim g1 (z; = 00 dx

---

Published in Transactions on Machine Learning Research (10/2023)
holds for at least one E {l....d}. a' (xa/t) o 0 for any (such sequence exists since a is periodic and non-constant). Then (22 indeed limits to positive infinity, giving the desired result.
d.. Take y and let vanish monotonically to zero in such way that
Without giving a formal proof, we now briefly describe sufficient conditions under which Theorem 3.2 could extended the more general case n (assuming one relains the positivity assumption (88 on do / da and the uniform boundedness assumption 19 on the network parameters). A neural network with one hidden layer of width d would be
n
u(x; 0) = W
w kI xi + b() + b();
here 27 denotes the I-th component of the function input x E R". Following the same argument just given, define
gx(2;)
W 00) 21 + 1
and assume that there exists some 1 μ S and 13 V S n such that W is asymptotically larger than (as before, no such entries exist, then W limits the zero matrix as c 4 O, and hence u(x; 0) limits constant function). The Frobenius norm of the PINN's neural tangent kernel matrix will blow-up if [L ga (2;0)² 0) blows-up at some collocation point This can occur if the entries for which V are nonzero coincide with the entries where (7()) is larger than O (e); more precisely, if there is sequence I 0 such that
Σ dai o'
dz, )
Σ 1 at Iı + b00) w nu #0
for any collocation point x, then the main result will generalize.
Returning to the case of one-dimension (n = 1), we now extend Theorem 3.2 to more general fully connected neural networks under an additional assumption on the network parameter's 's asymptotic bchavior.
Theorem 3.3. Let R + E. be fully connected neural network with A hidden layers of widths d, d2,.... .2....,d,, so that
(00) = x, a (x)
(x) = o (W + b0),
1 15ISA,
u(r;0) 0) W(A+1) (M) (x) + ((++1).
(25)
If all of the assumptions from Theorem 3.2 are retained, but T is now defined to be the set of E [0, T] such that the magnitude of every entry of each matria W 1 1 < A, are asymptotically larger than e1/A then
Nt E T
lim (+]l 0o, <10
where K the NTK mairia associated the loss function of the boundary value problem 15)
The proof proceeds nearly identically to that of Theorem 3.2 and hence it is Appendix C
An immediate corollary of Theorem 3.3 is that the spectral radius of K must also blow up 0; in contrast, because there is no explicit e-dependence in the K& subblock, one expects its spectral radius to be bounded independently of which is indeed numerically observed in section 4.1
Corollary 3.4. Let p( Ke (t)) denote the spectral radius of (t)(()- Under the same assumptions as Theorem
3.3 we have for any E T that
lim lp( K (t)) oc.
+.0
See Appendix D for the simple proof.
8

---

Published in Transactions on Machine Learning Research (10/2023)
3.2 Discussion
Fundamentally, using PINNs to numerically solve PDEs involves optimizing multiobjective loss functional, for example of the form
( = LPDE( + LBC(O)
for boundary values problems. The present work builds on the previous studies Wang et al. (2021a). Wang
et al. (2021b) and Wang et al. (2022). where certain PINN failure modes were identified for multiscale problems that, broadly speaking, involve inbalances between CPDE and LBC.
These provious analyses focused unsatisfactory PINN training in the particular case of Poisson-type PDEs, i.e. where the differential operator is the Laplacian (or, in the one-dimensional case, the second derivative). The multiscale character of the problems originated from oscillatory forcing functions; an illustrative example is the boundary value problem
de
dx u (x) (x) sin(x) + (1/c) sin(s/nn(////)
(26)
on [0, I] along with appropriate boundary conditions (see e.g. Eq. (2.6) in Wang al (0021b)). The present work focuses on elliptic problems that feature a more general class of divergence-form differential operators with oscillatory coefficients a" (z), as in Eq. (3] The differences between the two cases can be considcrable.
For example, from Wang et al. 2021a it is known that PINNs can fail to train when the magnitudes of the gradients of the two different loss components (CPDE and LBC) with respect to the network parameters are imbalanced, i.c. when either |VoLDDE > VooBcl or VVCPDE < |VoLBel. For multiscale Poisson equations such (26), the 1/E term in the right-hand side forcing results in VoLrDE/ scaling like 1/c. In contrast, for Darcy-type problems such as (3 (or, in one-dimension,
(15)). VaLpDEl is more singular; it scales like 1/c². In both cases, /aoncl is indepedent of c, so the discrepancy is greater for problems considered in the present work.
An alternative PINN failure mode analyzed in Wang et al. 2021bb and Wang et al. 2022 is when the the eigenvalues of the different NTK matrix subblocks differ greally in modulus. The authors reported a large imbalance between the eigenvalues of the Kuu and Kis subblocks for Poisson-style problems. However, for equations such (26), it is clear from the definitions in (14) that the Kau subblock is independent of c, since the forcing function f (which is responsible for the equation's 's multiscale character) is independent of 0.
In contrast, the results from Section 3.1 show that while the Kas subblock is independent of c for Darcy-type problems, the spectral radius of the Kau subblock increases as 1/e2. Hence, for problems with large scale separation, the imbalance eigenvalues is even more severe than in the Poisson case. For this reason, training PINNs with gradient based optimizers to achieve an accurate approximation equations of the form (3) becomes increasingly difficult function of increasing problem scale separation.
solutions to multiscale
A common strategy to deal with the imbalances just highlighted is to assign the different loss components CPDE and LBC different weights APDE and λBC; this was proposed e.g. in Wang et al. (2021a), ; van der Meer ct al. 2022) and Wang et al. (2022). It is reasonable, for example, for the equations considered here to try and weight LPDE by APDE = e to cancel the 1/e² factor in the NTK Kún subblock. Although this of course changes the magnitude of the cigenvalues of the NTK subblock, it does not change their distribution, and hence cannot directly resolve spectral bias Wang et al. (2022), which may explain the lack of success using adaptive weighting strategies that was reported in both Leung et al. (2022) and Han & Lee (2023). The numerical examples in Section 4.4 beloo illustrate that even under this rescaling, the stiffness of the gradient flow dynamics described by remains, translating to poor PINN performance.
4 Numerical results
4.1 Scaling of the neural tangent kernel matrix
Theorems 3.2 and 3.3 state that the Frobenius norm of the neural tangent kernel (NTK) matrix K' associated to the PINN approximation of the multiscale BVP 15) should scale inversely proportional to N. 2 for
9

---

Published in Transactions on Machine Learning Research (10/2023)
gradient descent training with infinitesimal time steps. Here numerically reproduce this scaling for neural network trained in more practical setting, namely with the Adam optimizer (Kingma & Ba, 2014).
As in Theorems 3.2 and 3.3 we focus on the Kin subblock of the NTK matrix. For a neural network with one hidden layer d, the 2j-th entry of Kin is given by
a
a
a
a L u(z);;(())
(K
Σ
L'u(z;;00
f(t)
L' (((zz;((tt)
L' ('u(a, 0(t)) ow
awa)
aw (L)
aw
θ
L (1 0(1) 0b
L u (xx;;((1))
(note (1; (); = 0). For network width d = 50 and N. = 256 collocation points, Figure 2(a) shows the Frobenius norm || associated to (5) for the particular case that ba [-x, Uh The norm is computed at time ( = 0, i.e. at initialization, for a sequence
0, and
(xx)) = 1/(2.1 + 2 sin((2xz/e) of e values between 1/10 and 1/100. For every e value, the neural network parameters are initialized from the normal distribution with mcan zero and unit variance. Figure B(b shows the Frobenius norms after initializing with the Glorot distribution (Glorot & Bengio 2010) and then training for ten thousand iterations with the Adam optimizer at learning rate of rj 10 In both cases, the norm increases 1/e², consistent with the theory developed above.
10a
Frobenius norm
Frabenius norm
À
À
A2
A2
10a
101
M/f
10a
10a
101
10a
10a
10
10a
(a)
1le
(b)
1/c
Figure 2: Scaling of the Frobenius norm of the NTK matrix and the first three principal eigenvalues (a) initialization and (b) after training with Adam optimizer.
of square matrix is equal to the Euclidean norm of its singular values. Since K is symmetric matrix, its singular values equal the absolute value of its (real valued) eigenvalues. Hence, Theorems 3.2 and 3.3 imply that { do as vanishes, where Xf are the eigenvalues of K Figures a) and b) further suggest that, for the BVP just described, the first few principal eigenvalues (denoted λι A2, and Aa) also increase at rate proportional to 1/e", both at initialization and after training with the Adam optimizer. In contrast, the spectral radius of the Kis subblock (not shown) is essentially independent e.
Recall that the Frobenius
For the particular case e = 1/80, Figure shows the full spectra of K both at initialization and after training. AL initialization, there are a handful of large eigenvalues, but the rest are clustered around the origin; after training, there is a single large eigenvalue well-separated from the cluster near the origin. In both cases the ratio X X is quite large, indicating that (12 is stiff initial value problem.
4.2 Learning two-scale function
Consider next the two-scale function
u' (x) sin(2x) + csin(x/e) sin(T/c)
(27)
π
10

---

Published in Transactions on Machine Learning Research (10/2023)
At initialization
10.0
1e-17
After training
7.5
7.5
5.0
(m(A)
Im(A)
0.0
0.0
2.5
-5.0
-7.5
-7.5
-10.0
5-10"
10"
1.5-10'
2-10"
2.5
10"
3-20"
1-10'
2-10'
3-10"
4-10
5-10"
(a)
Re(A)
(b)
Re(A)
Figure 3: The eigenspectrum of Kin. (a) at initialization and (b) after training,
This function satisfies the Poisson BVP
& de (x) 4 (2i) + sin(x//),
1
x <
(28)
with homogeneous Dirichlet boundary conditions u²(- = = 0: it additionally satisfies the elliptic, Darcy BVP given by
d d vl (x)) (z), dx (o(z/c) de
(29)
for x E T), again with homogeneous Dirichlet conditions, where
a(0/e) (2.1 2isn(z/)))
and the expression for the right-hand side forcing (x) is listed in Appendix E
We describe three different neural network approximations to the two-scale function u' u'()) by:
(1) regression (i.e. supervised learning),
(2) PINN solver for the Poisson problem 28 and
(3) a PINN solver for the Darey problem 29
We first report the results for each case when e = 1/32 and then repeat the experiments in Section 4.3 at smaller values to observe the behavior cach approximation as a function of increasing scale scparation.
For = 1/32, the first two methods listed above generate satisfactory approximations (27 indicating that (a good approximation to) the target function lives in the "span" of the neural network used; the third method does not. All of the numerical tests implemented in the open-source DecpXDE package (Lu)
et al. 2021b;, and in each case, the neural network weights and biases are randomly initialized with the Glorot distribution (Glorot & Bengio 2010
Consider first approximating uf (z) with a noural notwork by rogression, which we denote as (x; For [a E [-TT, x], the network parameters are determined by minimizing
(0 Σ (u - (x:))²
(30)
We 1/32 and use neural network with hyperbolic tangent activation functions and four hidden layers, each of width d = 40. Using N [ 12m/(e/2) 403 equispaced training points, the network is trained
11

---

Published in Transactions on Machine Learning Research (10/2023)
with the Adam (Kingma & Ba, 2014 and L-BFGS (Byrd al. 1995 optimizers for ten independent trials. In each test, Adam is used for 20, 000 and 10, 000 iterations with learning rates of 10-2 and y 10-3, respectively, and L-BFGS is then used, on average, for approximately 3600 steps.
Let up (z; denote the neural network averaged over all ten trials. Taking {= to be a set of equispaced points from - to v, Figures maximum error is approximately
a) and (b) show the resulting approximation and error, respectively. The
max Jeal == max Ju (x:) - E un (x;; 0) = 2.65 10-3
and the average variance in the over the interval [- æ²] is about 3 1007
2.0
10-2
10-1
10-3
0.0
-0.5
10:3
-1.0
10
-15
10
-2.0
-3
-2
(a)
(b)
Figure 4: (a) For c= 1/32, the two-scale function u (x) given by (7 and its approximations by regression (E un), a Poisson PINN (E up) and a Darcy PINN (E up); in each case the neural network contains four hidden layers of width d = 40. (b) The absolute value of corresponding generalization errors en, ep and ep 1000 equispaced points from
Next, consider approximating the solution to the Poisson problem function
28 with a PINN up(x; 0). The loss
N (((0)= No
d² dx
up (2x;0) (4s)n(2x,) + sin(x;/()))
(up(-x:0)} (Tp(T)))²))
(31)
For 1/32 and same neural network before, train with N, 1024 equispaced collocation points for ten independent trials. The training schodule consists of 10, 000, 10, 000, 20, 000, and 20, 000 itcrations of Adam with learning rates of 1 = 10-4, 10-5, = 10-6, and y 10-7, respectively, and then L-BFGS for few thousand iterations, on average.
Figures 4(a) and (b) show the resulting approximations and the generalization errors of the neural net- works alter training. Taking again E Eu;(2;;) 0) be the network averaged over all the trials, the maximum generalization error is
max
lepl max Ju (x:) - E up (z:; 0)] Rs 7.20 10
while the average variance across the interval is about 1.12 10-1
29 where the loss function is Finally, consider the PINN approximation uD(a;;) of the Darcy problem
1 No d (a(as/e) wn (2)=)) + I (axe) + (Ind (in(-==0)] + (MD(TT:))} de
(0)
Σ
(32)
N. 2 drr 1=1
Once again, 1/32, and the network architecture and number of collocation points is the same as the Poisson case. The training schedule consists of 10, 000, 10, 000, 20, 000, and 20, 000 iterations of Adam with
12

---

Published in Transactions on Machine Learning Research (10/2023)
C
c [ (0)
|aarp tx
(epll100000)
|Vrr 1 (80)
1/64
1.53 10-
5.42.
10-2
4.21 10-3
1.67
10-3
1/128
1.60 10-1
9.83
10
4.00 10-3
5.26
10
Table 1: co-norm of the gencralization error for the Poisson PINN throughout the domain Ω = [-x,, x. T] and on the domain boundary an for E 1/128.. Also shown are the maximum variances in the generalization errors, both in Q and on the boundary 02.
learning rates of 17 10-3, 17 10-4, thousand iterations, on average. The maximum generalization error
4 n 10-5, and 1) 10 respectively, and then L-BFGS for few
max lepl := max (z,) - E Eup(z;;0) ~ 4.31 10-1
which is about 40,, and the average variance across the interval is nearly as large about 2.83 10-1, From Figures 4 a) and (b) it can be seen that the error contains both low and high frequency components.
4.3 Results for decreasing E
We now consider the results for the three different neural network approximations to 27 at the values = 1/64 and e 6= 1/128.
Consider first the case of simple L regression (i.c. minimizing (3 For both values of c, N [2 / (c/2) equispaced training points are used for a network with four hidden layers, each of width 40. Ten independent trials conducted where the network is trained with the Adam optimizer for 20,000 and 10,000 iterations at learning rate of 17 10 and n 10-4, respectively, The generalization error is tested at 1000 equispaced points in [].
The results (not shown) are consistent with the frequency principle, i.e. the known spectral bias of neural networks in the regression setting Xu et al. (2020). The oo-norm of the average generalization error is about 1.71 10-2 and 8.88 10-3 for e 1/64 and 1/128, respectively, while both cases the maximum pointwise variance relatively small, on the order of 10-6, 6 The Fourier transform of the error in both cases is essentially zero everywhere except at the 1/c frequency, expected. The techniques proposed in Yang et al. (2021) for example could be used to mitigate this issue in this setting, if desired.
For the Poisson and Darcy PINNs, we consider a neural network with eight hidden layers, each of width d = 40. The training schodule consists of 2000 iterations of Adam with a learning rate of n = 10-6, 5000 iterations with f = 10-7 10000 iterations with 7 = 10-8 and then L-BFGS. The pointwise generalization errors for the two cases e = 1/64 and 6 = 1/128 are shown Figures 5 and i in Appendix F for the Poisson and Darcy PINNs, respectively. The co-norms of the average generalization error for each e case are reported in Tables and for the Poisson and Darcy PINNs, respectively. Also reported are the maximum average errors on the boundary an (). as well as the co-norm of the variance in the generalization error throughout 12 and on 0n.
Neither PINN does particularly well, as expected from the theory developed here for the Darcy case and the analysis in Wang et al. (2021b: 2022) for the Poisson case. Although there are some similarities in the error profiles-for for example, the variance over the independent trials becomes larger for both PINN cases as decreases there are some important differences to highlight well.
For the Poisson PINN, the average generalization error is approximately the same at both values of and in particular, the PINN satisfics the homogencous Dirichlet BC to about three digits of accuracy for cach e. To ensure this, we had lo set the boundary weight λnc = 10 (this value was kept consistent for all the cases listed in Tables [ and 2 In contrast, the average error in the Darcy case increases for decreasing e, and particular, the error at the boundarics increases. These results are consistent with the discussion Section 3.2 i.c. that the oscillatory coefficient within the differential operator for the Darcy PINN worsens the imbalance between the two loss terms LTDE and Lue already known to exist for Poisson PINNs.
13

---

Published in Transactions on Machine Learning Research (10/2023)
<
[ (x(
[Varp t (0)
[eD]IL=
|Varp (00(0))
1/64
1.01 10-1
8.43
10-2
7.84
10-2
8.34
10
1/128
3.80 10
3.06
10
3.11 10-1
2.99
10
Table 2: oc-norm the gencralization error for the Darey PINN throughout the domain Ω = [- π and on the domain boundary an for e E (1/64, 1/128}. Also shown are the maximum variances in the generalization errors, both in Q and on the boundary 0N.
4.4 Rescaling the loss function
As mentioned in Section 3.2 it is reasonable to try mollify the imbalance between the two parts of the Darcy PINN loss function by rescaling the PDE residual portion LPDE by c², as this cancels the 1/2² factor in the NTK matrix K subblock. Here we test this idea at a sequence of values.
We consider neural networks with three hidden layers each of width = 70. As before, we use hyperbolic tangent activation functions with N: 1024 collocation points. The training schodule consists of ten thousands iterations each of the learning rates v 10-3, 10-1, 10-5 and 10-6, and then, on average few thousand iterations of L-BFGS. As before, conduct ten independent trials and then average the results.
In particular we consider lcarning the solution function
u' (x) e (x-1)² + sinn(x//)) + cx
where er and 02 are chosen such that of satisfies homogeneous Dirichlet boundary conditions the domain [-m, m]; with
(x) = 1.1 + cos(x/e) cos(xx/e),
we set the forcing function fr in the Darey problem such that o' indeed solves
d a (x)
dr u0)) (x)
d
(33)
dx
The pointwise generalization errors for three values of e, as well as the magnitude of their Fourier coefficents, are shown Figure 7 in Appendix F Compared to the Darcy PINNs without rescaling, the homogeneous Dirichlet conditions are much better satisfied. For the = 1/36, the average error on the boundary is on the order of 10-3, consistent with the Poisson cases in the previous section. For 1/72 and e 1/144, the accuracy is lower (on the order of 10-1), but the variance is quite low, in contrast to the huge variance observed in the unscaled Darey PINNs. In each case, the error in Fourier space has both low frequency peak concentrated in the first few wavenumbers and a peak at the 1/c frequency, on average the errors in these peaks increase for decreasing e. In particular, at the smallest value = 1/144, an additional low froquency peak emerges around k = 16.
As discussed extensively Wang et al. (2022 and mentioned Section 3.2 both the magnitude and the distribution of the eigenvalues of the NTK matrix subblocks are closely connected with PINN performance. As shown in Figure 3 there is large discrepancy between the principal cigenvalue of the K subblock and the rest of the spectrum, which is characteristie sign that the gradient flow dynamics 11] are stiff. Although the rescaling introduced by setting ÄTdE = e changes the magnitude of the eigenvalues of the K subblock, it of course does not change their distribution, and we attribute the increasing errors at decreasing values to the discrepancy in the spectrum.
5 Conclusions
Physics informed neural networks have demonstrated success in a wide variety of problems in scientific computing (Karniadakis et al. 2021), however, they can sometimes struggle to approximate solutions to dif- ferential equations with multiscale features. Athough PINN convergence theory guarantees that minimizers of (regularized) loss functions in principle well approximate solutions to elliptic problems, in practice,
14

---

Published in Transactions on Machine Learning Research (10/2023)
achieving those minimizers is difficult for boundary value problems with highly oscillatory coefficients present
in the differential operator.
We show that for a class of linear, multiscale elliptic equations in divergence form, the Frobenius norm of the neural tangent kernel matrix associated to the PINN becomes unbounded as the characteristic wavelength of the oscillations vanishes. Numerical examples illustrate that during training, the ordinary differential equation that governs the evolution of the PDE residuals becomes increasingly stiff & 0, translating to poor PINN training behavior.
The present work considers standard, fully connected neural network architectures; it may be possible in future research to develop alternative PINN architectures that are specifically adapted to the problem. For example, different activation functions with problem-tuned inductive biases can be used, as suggested in Ziyin et al. 2020) for periodic functions, Li et al. (2020a) also proposed smooth and localized activation functions for multiscalc clliptic problems the context of the Deep Ritz method Yu ct al. (2018). In this context, it may be beneficial to develop a theory for the Neural Tangent Kernel matrix for energy-minimizal based neural network methods for solving multiscale PDEs, the derivative of the oscillatory coefficient a does not appear, in contrast to PINNs.
Future work may also benefit from more sophisticated learning techniques beyond (stochastic) gradient descent, for example genetic algorithms (Mirjalili, 2019). Another possibility is to employ continuation methods (Allgower Georg 2003), i.e. "curriculum regularization" (Krishnapriyan et al. 2021
Recent results establishing "frequency principle" in machine learning motivated this study; during training, neural networks tend to learn their target functions from low to high frequency (Luo et al. 2021 See also the work by Wang et al. (2021b) for spectral analysis of physics informed neural networks applied to Poisson problems with multiscale forcing. An interesting avenue for future study is to characterize the spectral behavior during training of deep learning methods for learning operators; see for example (Li et al. 2020b; Khoo et al. 2021 Bhattacharya et al. 2021 Lu et al., 2021a Zhang ct al. 2022). In the context of the present work, this would be the nonlinear map from the problem data a" and ʃ in (3) to the solution u Finally, we note many other novel neural network approaches to multiscale elliptic problems
discussed here have recently been proposed, for example (Han & Loe, 2023 Leung et al. 2022 Fabian et al. 2022
We hope some of the results presented here might inspire insight into these exciting approaches.
Acknowledgments
The authors thank Andrea Bertozzi for organizing the UCLA Computational and Applied Mathematics REU site where much of this work took place. The authors also thank Chris Anderson and Michael Murray for helpful discussions.
References
Eugene I, Allgower and Kurt Georg. Introduction to Numerical Continuation Methods. STAM, 2003.
Alain Bensoussan, Jacques-Louis Lions, and George Papanicolaou. Asymptotic analysis for periodic struc- lures, volume 374. American Mathematical Soc., 2011.
Kaushik Bhattacharya, Bamdad Hosseini, Nikola B. Kovachki, and Andrew M. Stuart. Model Reduction And Neural Networks For Parametric PDEs. The SMAT journal of computational mathematics, 7:121-157,
2021.
William L Briggs, Van Enden Henson, and Steve F McCormick. A mulligrid tulorial. SIAM, 2000.
Richard Byrd, Peihuang Lu, Jorge Nocedal, and Ciyou Zhu. A limited memory algorithm for bound constrained optimization. SIAM Journal scientific computing, 16(5):1190-1208 1995.
Weinan E. The dawning of new era in applied mathematics. Notices of the American Mathematical Society, 68(4):565-577 2021.
15

---

Published in Transactions on Machine Learning Research (10/2023)
Bjorn Engquist and Erding Luo. Convergence of multigrid method for elliptic equations with highly oscillatory coefficients. SIAM Journal on Numerical Analysis, 34(6):2254-22)) 22733, 1997.
Kröpfl Fabian, Roland Maier, and Daniel Peterseim. Operator compression with deep neural networks. Advances in Difference Equations, 2022(1), 2022.
Olga Fuks and Hamdi A Tchelepi. Limitations of physics informed machine learning for nonlinear two-phase transport in porous media. Journal of Machine Learning for Modeling and Computing, 1(1,, 2020.
Xavier Glorot and Yoshua Bengio. Understanding the difficulty of training deep feedforward neural networks. In Proceedings of the thirteenth international conference on artificial intelligence and statisties, pp. 249- 256. JMLR Workshop and Conference Proceedings, 2010.
Jihun Han and Yoonsang Lee. neural network approach for homogenization of multiscale problems. Multiscale Modeling Simulation, 21(2):716 734, 2023.
Kurt Hornik, Maxwell Stinchcombe, and Halbert White. Multilayer feedforward networks are universal approximators. Newral networks, 2(5):359-33 -66, 1989.
Arthur Jacot, Franck Gabriel, and Clement Hongler. Neural tangent kernel: Convergence and generalization in neural networks. In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett (eds.), Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018.
Weiqi Ji, Weilun Qiu, Zhiyu Shi, Shaowu Pan, and Sili Deng. Stiff-pinn: Physics-informed neural network for stiff chemical kinetics. The Journal of Physical Chemistry A, 125(36):8098 8106, 2021.
George Em Karniadakis, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. Physics- informed machine learning. Nature Reviews Physics, 3(6):422-440, 2021.
Yuehaw Khoo, Jianfeng Lu, and Lexing Ying. Solving parametric pde problems with artificial neural net- works. European Journal of Applied Mathematics, 32(3):421-435, 2021.
Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412. 6980, 2014.
Aditi Krishnapriyan, Amir Gholami, Shandian Zhe, Robert Kirby, and Michael W Mahoney. Characterizing possible failure modes in physics-informed neural networks. Advances in Neural Information Processing Systems, 34:26548-2653 26560, 2021.
Isaac E Lagaris, Aristidis Likas, and Dimitrios Fotiadis. Artificial neural networks for solving ordinary and partial differential equations. TEEE transactions on neural nelworks, 9(5):987- 1998.
Moshe Leshno, Vladimir Ya Lin, Allan Pinkus, and Shimon Schocken. Multilayer feedforward networks with a nonpolynomial activation function can approximate any function. Newral networks, 6(6):861-888 11-67, 1993.
Wing Tat Leung, Guang Lin, and Zecheng Zhang. NH-PINN: Neural homogenization--aasel physics-informed neural network for multiscale problems. Journal of Computational Physics, 470:111539, 2022.
Xi-An Li, Zhi-Qin John Xu, and Lei Zhang. A multi-scale dnn algorithm for nonlinear elliptic equations with multiple scales. Communications in Computational Physics, 2885) 2020a.
Zongyi Li, Nikola Borislavov Kovachki, Kamyar Aviwadenesheli, Kaushik Bhattacharya, Andrew Stuart, Anima Anandkumar, al. Fourier neural operator for parametric partial differential equations. Tn International Conference on Learning Representations, 2020b.
Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, and George Em Karniadakis. Learning nonlinear operators via deeponet based on the universal approximation theorem of operators. Nature Machine Intelligence, 3(3):218-225 229, 2021a.
Lu Lu, Xuhui Meng, Zhiping Mao, and George Em Karniadakis. Deepxde: A deep learning library for solving differential equations. STAM Review, 63(1):208 228, 2021b.
16

---

Published in Transactions on Machine Learning Research (10/2023)
Tao Luo, Zheng Ma, Zhi-Qin John Xu, and Yaoyu Zhang. Theory of the frequency principle for general deep neural networks. CSIAM Transactions on Applied Mothemalics, 2(3):484- 507, 2021.
Stefano Markidis. The old and the new: Can physics-informen deep-learning replace traditional linear solvers? Frontiers in big Data, pp. 92, 2021.
Scyedali Mirjalili. Genetic algorithm. In Evolutionary algorithms and neural networks, pp. 43-55. Springer, 2019.
Dimitris C Psichogios and Lyle H Ungar. A hybrid neural network-first principles approach to process modeling. AIChE Journal, 38(1)):1499 1511, 1992.
Nasim Rahaman, Aristide Baratin, Devansh Arpit, Felix Draxler, Min Lin, Fred Hamprecht, Yoshua Bengio, and Aaron Courville. On the spectral bias of neural networks. In Kamalika Chaudhuri and Ruslan Salakhutdinov (eds.), Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pp. 5301-5310, PMLR, 09-15 Jun 2019.
Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informed neural networks: A deep learn- ing framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics, 378:686-700) 2019.
Yeonjong Shin, Jerome Darbon, and George Em Karniadakis. On the Convergence of Physics Informed Neural Networks for Linear Second-Order Elliptic and Parabolic Type PDFs. Communications in Com- putational Physics, 28(5):2042 2074, 2020.
Remco van der Meer, Cornelis W Oosterlee, and Anastasia Borovykh. Optimally weighted loss functions for solving pdes with neural networks. Journal of Computational and Applied Mathematics, 405:113887, 2022.
Sifan Wang, Yujun Teng, and Paris Perdikaris. Understanding and miligating gradient flow pathologies in physics-informed neural networks. STAM Journal on Scientific Computing, A3055: 2021a.
Sifan Wang, Hanwen Wang, and Paris Perdikaris. On the eigenvector bias of fourier feature networks: From regression to solving multi-scale pdes with phfsics-ifformed neural networks. Computer Methods in Applied Mechanics and Engineering, 139:1, 2021b.
Sifan Wang, Xinling Yu, and Paris Perdikaris. When and why pinns fail to train: A neural tangent kernel perspective. Journal of Compulational Physics, 449:110768, 2022.
Zhi-Qin John Xu, Yaoyu Zhang, Tao Luo, Yanyang Xiao, and Zheng Ma. Frequency principle: Fourier analysis sheds light on deep neural networks. Communications in Computational Physics, 28(5):1746- 1767, 2020.
Ge Yang, Amurag Ajay, and Pulkit Agrawal. Overcoming the spectral bias of neural value approximation. In International Conference on Learning Representations, 2021.
Bing Yu et al. The deep ritz method: a deep learning-based numerical algorithm for solving variational problems. Communications in Mathematics and Statistics, 6(1):1-12, 2018.
Lulu Zhang, Tao Luo, Yaoyu Zhang, Weinan E, Zhi-Qin John Xu, and Zheng Ma. MOD-Net: A Machine Learning Approach via Model-Operee Model--pperato--aal Network for Solving PDEs. Communications in Compula- tional Physics, 32(2):299-33) 335, 2022.
Liu Ziyin, Tilman Hartwig, and Masahito Ueda. Neural networks fail to learn periodic functions and how to fix it. Advances in Neural Information Processing Systems, 33:1583-115 1594, 2020.
17

---

Published in Transactions on Machine Learning Research (10/2023)
A Description of PINN approximation for multiscale Poisson problem
Here we give a more complete description of the example from Section of a PINN approximation to the Poisson boundary value problem 4
with homogeneous Dirichlet boundary conditions.
The PINN solution computed with a fully connected neural network with four hidden layers of sixty nodes
each and hyperbolic tangent activation functions. The training consisted of 100, 1, 000, 15, 000, 30, 000 and 10, 000 iterations with the Adam optimizer with learning rates of v) 10-2 and = 10-6, respectively. We note, however, using constant value of = 10-3 throughout training produced results qualitatively similar to those in Figure ! i.e. the observed frequency principle was fairly robust to changes in the learning rate. As mentioned in the description in Section 2 for legibility, the results shown in Figure are that of the rolling average (of 50 values) to the trajectories of the error in the Fourier coefficients.
10-3
10-4, 17 10-3,
Finally, mention that the end of the training process, the co-norm error between the true solution and the neural network approximation is about 1.26 10-2
B Proof of Lemma 3.1
Proof. Computing first the entries of the gradient of the loss function
8 gives
Oc
1 N.
i L' L'u(zz;0))
A
N₂
2
pdd (zy; 0)
rb
0)
u(s;:
(34)
00,
No
00,
No
00,
5=1
Next, compute
N₂
d pde "ddooo(; (x:; 0(1) (L'u(z00(1)) - f(x)))
d
d0, A
Σ
L' L'ula;; 0)
di
di
a0,
1=1
N2
Oc
L' ulu 0)
(35)
00,
00,
where the final equality comes from the assumption of gradient flow dynamics gives
11
A similar computation
N₂
d
d
d6: di
Li
(()))((()))
(u(s;; 0())) g(s:))
Σ
u
s;;
0)
di
002
1=1
Np
ac
ul
8; 0)
(36)
00,
30,-
1=1
Inserting
34
into 35
and rearranging the sums gives
pdd (x:; (24; 0(1))
-
T. [,;; 0)
θ I''(xxy; 0)
ppde ppde 0)
dt
N
at
00,
N
a
a u(s;; ())00)()00))))
L u'u
0)
N
20,
00,
1=1
while inserting 34 into 36
similarly results in
N
N
To 0()) =
u( 0)
dt
00,
00,
1=1
No.
a
u(s;; 0)
u(s;; 8) no(s);)))
N
ae
00
==1
18

---

L. (((.a:0(1))) Published in Transactions on Machine Learning Research (10/2023)
desired.
C Proof of Theorem 3.3
Proof. From 25 we note that W E Iidxi W e Ráxd-1 for 1= ..... and W(++I) E ERlxds while ((4) E Rdd for
I,..., A and b(A+1) E I. The fully connected neural network 25 can then be written as
u(2; 0)
w(++1) (A) (I) + b((+1))
Uk
kA=1
(A) (z)
W(Ai Kaka-1 "Ka-1 (1-))(()) (z) + b
UKA
1s S da.
kx-1=1
da-2
(A-1)(2) (x) =0 Σ
W(d-1)
ka ihs UA (x) + ((2-11))
1 1s KA-1 S d
hx-2=1
d1
2x (z) o
Σ
w (2)
Kahi 4ki
+ bb6))
1 1 S d2
=1=1
"k
(1)
(x)
(w x + (1))
1 h1 s d1.
By lincarity, note that
0
L' u(x; 0) (c)
(37)
ow
for any 1 IS S da. As before, consider diagonal entry in the K subblock of the NTK matrix 14
N2
1
à
à 00,
(Kar)an )aa (1)
L (1
N
00,
By
37 sufficient condition for the Frobenius norm of K to blow-up, then, is for
lim <10
[L u (A) (xa)) 00
for some collocation point Ta, some S y S dr, and t t T (here the dependence of u(A) on and O is implied). Since
= (a((a//0) de dr u (xw))
daz u (xa) ( (xa/E) d 2 (d((xa/e)
2
aa<xa/t) ()) d²
de (a)))
1
dx wn ((N) 3
+
(38)
and since is a periodic C'l function, the desired result will follow if (i) the second derivative of u.y is uniformly bounded in e and (ii) the first derivative is asymptotically larger than e, so that
u (x) = oc
lim <j0 dr
(39)
for any E T.
19

---

Published in Transactions on Machine Learning Research (10/2023)
The first derivative is after applying the chain rule once. Using Einstein sunmation notation for brevity,
da-1
d
(x)
2 (x)
Σ
( R (1) +
(A)
Σ Ax-1=1
(
(x))
(40)
dx
2ka-1
pha-1 do u
AA-1=1
40
simplifies to
d (x) o' (wal) - (x) b/2))
(W(Ay 7ha-1 dr (x)-)(())
d
dæ
Contimuing with the chain rule (and the Einstein notation), we get (x) + ((a)))
dr my (x) = o' WTTa
d
-haa- "hai (a-1)(
(W(N) 7ka-1 From here we can recurse downwards to compute the full derivative as
of (w Ra- 2 (x) + ((--11) (W/-D ka-ska-o dx K (z)
d (A-2)
(x) = (W ska-1
(A-1) (2) + b//n)) x
dæ Whan o' d (x) + ( w (1-1) ha- ka-2 Wha
K o' d'owW/l (wWl kn-skars UKa-s (z) + ((1-2))
(A-3)
w-(3) (3) kskk laki "k (x) (
of (w(2) (1)
(wel kskk d (Wo)) (1)
o'
10)
II]
(41)
By the uniform boundedness assumption 19 on the network parameters 0, 41 is uniformly bounded in from above for any I E [a, b]. As a consequence of the positivity assumption I8 199 further implies that
each instance of do dr (41) is bounded below by a positive number that is independent of Note that 4ll generally contains sums and products of da dr and network weights (entries in the W(!) matrices). For that
any E Y, the weights asymplotically larger than e//A,
C
lim <L0
=0
(2
((A-1) (](A) ha ika Kaka
b1
A2ki
for any combination of indices 1 k1 < d1, < 222 S d2, and so meaning that 339 indeed holds.
Ergo, as in the proof of Theorem 3.2 as long as e vanishes monotonically to zcro in such a way that a' ((ta/e) o any 6, (38) limits to posilive infinity as desired.
D Proof of Corollary 3.4
Proof. From elementary linear algebra, we know that the Frobenius norm of the square matrix K (t) equals the Euclidean norm of its singular values, and because K (t) is a symmetric matrix, its singular values equal the absolute value of its (real valued) eigenvalues. (AE(t) denotes the eigenvalues, then
No K (f)lF (E 1/2 S (Noo(((
giving the rosult.
E Boundary value problem forcing functions
The right-hand side forcing function (a) for the Darcy BVP 29 with homogenenous Dirichlet boundary conditions is given by 9d (x)/h((x)) where
9 (xxx) 10 20 n 168 me cos (x) sin (x) 20 2T 4n cos? esin(T/e)) cos(x/e)
+ 21T + 160 me cos (x) sin sin(z)) sin(x//e)
20

---

Published in Transactions on Machine Learning Research (10/2023)
and
h (z) 400 (x//)) + 840 me sin(2/e) + 441 me.
F
Error plots for decreasing
Poisson error, 1/64
Poisson error, == 1/128
0.2
0.2
0.0
0.2
-0.
(a)
(b)
Pointwise generalization errors for the Poisson PINN 31 trial runs, while the solid line with circle markers indicates the mcan.
Figure 5:
The dotted lines indicate independent
Darcy error, 1/64
Darcy error, 1/128
0.2
0.25
0.0
-0.25
-0.2
0.75
-0.4
-1.25
(a)
(b)
Pointwise generalization errors for the Darcy PINN 32 trial runs, while the solid line with circle markers indicates the mean.
Figure 6:
The dotted lines indicate independent
21

---

Published in Transactions on Machine Learning Research (10/2023)
t 1/36
== 1/36
0.2
SUD
10a
101
-0.1
10-
-0.3
10"4
-0.5
10
100
125
150
175
200
1/72
1/72
- EUD -
10
-0.1
10
-0.2
10a
-0.3
10-1
-3
-2
25
75
100
125
150
175
200
= 1/144
1/144
0.25
- e EUD
10'
0.15
0.05
100
10
-0.10
3
25
100
125
150
175
200
33 with the PDE residual part of the loss function weighted by λPDE = e2. The left columns display the pointwise generalization error on [-x, v), and the right columns display the magnitude of the first two hundred Fourier coefficients of the error. The dotted lines indicate independent trial runs, while the solid line with circle markers indicates the mean.
Figure 7: Results at three different, values for the Darcy PINN associated to
22


Reviews:
[{"review_id": "MZAi5U6Z3p", "reviewer": "Reviewer_c8m4", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "An article is a pleasure to read. It tells a coherent story, contains good references, and clearly explains the problem. It is possible to improve the presentation in several places, but overall the results are well-explained. In addition to theoretical results, the article contains a few simple experiments that support the claims made by the authors.\n\nOn the negative side, obtained theoretical results contain somewhat artificial assumptions. It is unclear to what extent the author's findings explain the pathological behavior of physics-informed neural networks applied to problems with oscillatory coefficients.\n\nBesides, I cannot agree with the central claim from the abstract \"We show that if the coefficient in the elliptic operator contains frequencies on the order of $1 \\big/ \\epsilon$, then the Frobenius norm of the neural tangent kernel matrix associated to the loss function grows as $1 \\big/ \\epsilon^2$\". The results do not seem to support it.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors consider the training dynamics of physics-informed neural networks. More specifically, they show that for elliptic equation (in flux-conservative form) with strongly oscillating coefficients $a_{\\epsilon}(x) := f(x\\big/\\epsilon),$ $f(y) = f(y+1),$ $f(y)\\in C^{1}([0, 1])$, neural tangent kernel contains entries of order $1\\big/\\epsilon^2$. Next, they conclude that ODE describing the evolution of residual becomes stiff for sufficiently small $\\epsilon$, and this is the main reason why training stagnates. Several numerical experiments support theoretical statements.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "**Questions:**\n\n1. Claim \"We show that if the coefficient in the elliptic operator contains frequencies on the order of $1 \\big/ \\epsilon$, then the Frobenius norm of the neural tangent kernel matrix associated with the loss function grows as $1 \\big/ \\epsilon^2$\" is not supported.\n\nFirst, Theorem 3.2 is more restrictive, since it applies only to the periodic functions of the special form $a_{\\epsilon}(x):=a(x\\big/\\epsilon)$.\n\nSecond, I am not sure that the claim as stated is accurate. For example, lets consider $a_{\\epsilon}(x) = \\epsilon \\cos(2\\pi x \\big/ \\epsilon)$. For any fixed $\\epsilon$ this term contains frequencies of order $1\\big/\\epsilon$, yet the neural tangent kernel is not going to contain terms of order $1\\big/\\epsilon^2$ since the absolute value of the derivative is bounded by $2\\pi$\n\nI suggest somewhat reformulating the sentence in the abstract or providing a stronger theoretical result.\n\n---\n\n2. Theoretical result is too narrow. It is not clear whether it is applicable in a broader context.\n\nAs I understand, the authors want to study PINN for elliptic PDEs with oscillating coefficients. Yet the theoretical result applies only to the specific class of PDEs and specific loss. Here are several suggested directions that can make the article applicable to a broader set of problems:\n\na. What if $a_{\\epsilon}(x)$ is such that the derivative is bounded? Is it going to affect the training or not? For example, we can consider a rescaled version of equation (28) from the article:\n\n$$\n-\\frac{d}{dx}\\left(\\frac{\\color{red}{\\epsilon}}{2.1 + 2\\sin(x\\big/\\epsilon)} \\color{black}{\\frac{d}{dx} u(x)}\\right) = f^{\\epsilon}(x) \\color{red}{\\epsilon}\n$$\n\nAuthors observed poor training with a version that is not rescaled. After the rescaling tangent kernel does not contain $1\\big/\\epsilon^2$. Will this simple rescaling improve training?\n\nb. Consider a non-conservative form of elliptic equation\n\n$$\na_{\\epsilon}(x)\\frac{d^2}{dx^2} u(x) = f(x).\n$$\nNow tangent kernel again is well defined for $\\epsilon\\longrightarrow 0$ as long as $a_{\\epsilon}(x)$ is uniformly bounded. Is it right, that the training will be efficient for strongly oscillating $a_{\\epsilon}$?\n\nc. What if we use the variational form (e.g., as in the Deep Ritz method)?\n\nIn this case, we do not need to consider the derivative of $a_{\\epsilon}$. Does it mean that the neural tangent kernel is well-defined? Are we going to see poor performance of PINN?\n\n---\n\n3. Experimental evidence is not sufficient.\n\nThe approach taken by the authors in Section 4.2 is quite elegant. Yet, the results are not convincing. I suggest repeating the same experiments (learning in $L_2$ setting, PINN+Poisson, PINN+Darcy) for a range of $\\epsilon$ and showing errors for different methods.\n\nBesides, I suggest adding at least rescaled version of equation (28) as explained in 2. above.\n\n---\n\n4. Minor suggestions.\n\na. In equation (6), it is not clear right away that the index $k$ refers to the coefficient with frequency $k$. Please consider clarifying this part.\n\nb. In Figure 1, a logarithmic scale in $y$ seems more appropriate.\n\nc. Neural tangent kernel should be symmetric positive semidefinite. Yet in Lemma 3.1 $K_{ub}\\neq K_{bu}^{T}$. It seems there is a typo in the coefficients.\n\nd. Page 6 \"Although it does not hold for the ReLU function  $\\sigma(x) = \\max(0, x)$, these are of course not suitable for NN-based solutions to PDEs with second order (or higher) differential operators.\" Perhaps, only for a strong form of PDE. Neural Networks with ReLU functions form a nonlinear approximation space of piecewise linear functions with variable positions breaking points. They cover the space of linear finite elements.\n\n---\n\nI encourage the authors to respond to the three main points above with additional experimental data or theoretical results, or rebuttal. In general, I believe that the author's contribution is valuable and suitable for publication in TMLR after the above issues are settled, one way or another.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "9egpAzFnnC", "reviewer": "Reviewer_DENr", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strength: No\n\nWeaknesses: The major weakness of the paper is a lack of sufficient novelty. The fact that GD is difficult to learn oscillatory functions with high frequency is known in the literature; see Xu et al. 2020 for regression setting and other references mentioned by the authors. In the PDE setting, Wang et al. 2021 and Wang et. al 2022 also studied the convergence of training process as well as the spectral bias from an NTK perspective. The only difference between the present paper and these earlier works lies in the investigation of the blow-up of the norm of the NTK matrix as shown in Theorem 3.2 and Theorem 3.3. The proofs of both results follow from straightforward calculations of the NTK matrix and do not seem to require new technicalities or ideas. \n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper investigates the pitfalls of Physics informed neural networks for solving multiscale elliptic PDEs where the solutions exhibit oscillatory behavior. The authors show that when the coefficient of the elliptic operator contains frequency at the order of $1/\\epsilon$, the Frobenius norm of the NTK matrix associated to the loss function scales like $1/\\epsilon$, whereby claiming that training the networks with gradient-descent based methods becomes increasingly difficult for highly oscillatory problems. Numerical experiments are provided to illustrate the stiffness of the problem. ", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "No", "audience_explanation": "", "requested_changes": "I do not have many comments on the possible changes. But I do want to add that it is certainly an interesting and important problem to design efficient deep learning methods for solving multiscale PDEs. Unfortunately, the present paper only points out the issue of PINNs and could be substantially improved if the authors were to propose new algorithm/methodology that can solve/alleviate the issue.  ", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "fv94MihMWR", "reviewer": "Reviewer_37N7", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The strength is the novelty of the idea that studying NTK can help us understand PINN and that the theory feels sufficiently sound.\n\nThere are a few weaknesses I think this work needs to address\n\n1. The main result states that the NTK change by the order ($\\epsilon^{-2}$). But I am a little confused by what this statement means. Does it mean that before and after training, the NTK will change by order $\\epsilon^{-2}$? Or it simply means that before training, the NTK will change because of the change in the loss function of order $\\epsilon^{-1}$\n\n2. The argument that NTK changes as $\\epsilon^{-2}$ \"implies that... training the neural network with gradient descent based methods to achieve an accurate approximation of the solution to the PDE becomes increasingly difficult\" is tenuous. I do not see why the claim is implied, the authors need to explain this in detail. After all, I think this claim cannot be supported by the theory of the paper. NTK is only relevant to the case when the learning rate is vanishingly small and is not relevant when the learning rate is large. In fact, this work might imply that it is better to use a large learning rate, a point that is not discussed in the main text\n\n3. it seems that the problem of training identified in this paper already has a solution. The authors should discuss the relevant work of https://arxiv.org/abs/2006.08195.  The problem with PINN is due to the fact that the assumed model activation function is monotonic (Eq. 18), whereas learning the right solution requires a periodic function. It seems to the audience that the problem simply comes from the fact that the inductive bias of the model does not match the problem. Including a discussion of the result in https://arxiv.org/abs/2006.08195 and performing numerical experiments to either validate or disprove what I state above will greatly improve the manuscript", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This work proves that the NTK of PINN increases as $\\epsilon^{-2}$ where $\\epsilon^{-1}$ is the frequency of the highest mode in the differential equation", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Please answer the weaknesses I raised and adapt the draft accordingly.", "broader_impact_concerns": "NA", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: Novelty was a concern in the paper, as pointed out by one reviewer, given existing studies that show the difficulty of training imbalanced multi-objective loss. However, the authors clarified that the present work focuses on elliptic problems featuring a more general class of divergence-form differential operators with oscillatory coefficients compared to related studies. This discussion has been added in the revision, making the position of this paper in the literature clear. Moreover, the authors have addressed other concerns raised by the reviewers. According to the evaluation criteria, evidence supporting the claim is more important than novelty and significance. The paper meets the acceptance criteria. Please take into consideration additional comments raised by Reviewer c8m4, which do not affect the decision, when preparing the final version.
Decision: Accept (poster)