Fast Convergence of Natural Gradient Descent for Over-parameterized Physics-Informed Neural Networks
Authors: Xianliang Xu
Wang Kong
JiahengM
Zhongyi Huang
Ye Li
Venue: ICLR 2026 Poster (2026)
Keywords: natural gradient descent
over-parameterization
physics-informed neural networks
neural tangent kernel
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)

Abstract: In the context of over-parameterization, there is a line of work demonstrating that randomly initialized (stochastic) gradient descent (GD) converges to a globally optimal solution at a linear convergence rate for the quadratic loss function. However, the convergence rate of GD for training two-layer neural networks exhibits poor dependence on the sample size and the Gram matrix, leading to a slow training process. In this paper, we show that for training two-layer $\text{ReLU}^3$ Physics-Informed Neural Networks (PINNs), the learning rate can be improved from the smallest eigenvalue of the limiting Gram matrix to the reciprocal of the largest eigenvalue, implying that GD actually enjoys a faster convergence rate. Despite such improvements, the convergence rate is still tied to the least eigenvalue of the Gram matrix, leading to slow convergence. We then develop the positive definiteness of Gram matrices with general smooth activation functions and provide the convergence analysis of natural gradient descent (NGD) in training two-layer PINNs, demonstrating that the maximal learning rate can be $\mathcal{O}(1)$ and at this rate, the convergence rate is independent of the Gram matrix. In particular, for smooth activation functions, the convergence rate of NGD is quadratic. Numerical experiments are conducted to verify our theoretical results.

Published as a conference paper at ICLR 2026
FAST CONVERGENCE OF NATURAL GRADIENT DESCENT FOR OVER-AARAMETERIZED PHYSICS-
INFORMED NEURAL NETWORKS
Xianliang Xu², Wang Kong", Jiaheng Mao', Zhongyi Huang' & Ye Lit.2.3
Coollege of Computer Science and Technology, Nanjing University of Aeronautics and Astronautics Schooo of Mathematics, Nanjing University of Aeronautics and Astronautics "MIIT Key Laboratory of Pattern Analysis and Machine Intelligence, Nanjing
'Department of Mathematical Sciences, Tsinghua University edu. cn, zhongyihhtsingeeee edw..in.nddcccc..
xuxl1 1 90ma i s.
[wkong, 1 1901 620583,22226632 yeli20 @nuaa. edu. en
ABSTRACT
In the context of pvereparaaettriizatinn, there is line of work demonstrating that randomly initialized (stochastic) gradient descent (GD) converges to a glob- ally optimal solution a lincar convergence rate for the quadratic loss function. However, the convergence rate GD for training two-layer neural networks ex- hibits poor dependence slow training process. In this paper, we show that for training two-layer ReLU* Physies-Informed Neural Networks (PINNs), the learning rate can be improved from the smallest eigenvalue of the limiting Gram matrix to the reciprocal of the largest eigenvalue, implying that GD actually enjoys a faster convergence rate. Despite such improvements, the convergence rate is still tied to the least eigenvalue of the Gram matrix, leading to slow convergence. We then develop the positive definiteness of Gram matrices with general smooth activation functions and provide the convergence analysis of natural gradient descent (NGD) training two-layer PINNs, demonstrating that the maximal learning rate can be O(1) and at this rate, the convergence rate is independent smooth activation functions, the convergence rate of NGD is quadratic. Numerical experiments are conducted to verify our theoretical results.
the sample size and the Gram matrix, leading to
the Gram matrix. particular, for
INTRODUCTION
In recent years, neural networks have achieved remarkable breakthroughs in the fields of image recognition He et al. (2016), natural language processing Devlin et al. (2018), reinforcement learning Silver al. (2016), and so Moreover, due to the flexibility and scalability of neural networks, researchers are paying much attention in exploring new methods involving neural networks for handling problems in scientific computing, especially in the area of solving partial differential
equations (PDEs) Müller & Zeinhofer (2023) Raissi et al. (2019); Yu et al. (2018): Zang et al. (2020); Siegel et al. (2023). Among them, the most representative approach is Informed Neural Networks (PINNs) Raissi et al. (2019). In the framework of PINNs, one incorporate PDE constraints into the loss function and train the neural network with it. With the use of automatic differentiation. the neural network can be efficiently trained by first-order or second-order optimizers.
First-order methods, such as gradient descent (GD) and stochastic gradient descent (SGD), are widely used in optimizing neural networks as they only calculate the gradient, making them computationally efficient. In addition to first-order methods, there has been significant interest in utilizing second-order optimization methods to accelerate training. These methods have proven to be applicable not only to regression problems, as demonstrated in Martens Grosse (015).. 5. but also to problems related to PDEs, shown in Müller & Zeinhofer (2033) Raissi et al. 2019
"Corresponding author: Ye Li <yel220@mm..chu.cm

---

Published as a conference paper at ICLR 2026
As for the convergence aspect of the optimization methods, it has been shown that gradient descent algorithm can even achieve zero training loss under the setting of over- aee--paraettriizatonn, which refers to situation where model has more parameters than necessary to fit the data Du et al. (2018) 2019);
Allen-Zhu et al.] (2019a b): Arora et al. (2019): : Li & Liang 2018); Zou et al. (2020); Cao & Gu (2019). These works are based on the idea of neural tangent kernel (NTK Jacot al. (2018). which shows that training multi-layer fully-connected neural networks via gradient descent is equivalent to performing a certain kernel method the width every layer goes to infinity. Despite these attractive convergence results, the learning rate depends on the sample size and the Gram matrix, so it needs to be sufficiently small to guarantee convergence in practice. However, doing so results in a slow training process. In contrast to first-order methods, the second-order method natural gradient descent (NGD) has been shown to enjoy fast convergence for the L regression problems as demonstrated in /Zhan et al. (2019); Cai et al. (219)), and PINN problems as in Müller & Zeinhofer (2023); Guzmán-Cordee et al. (2025). However, the convergence NGD in the context of training PINNs is still open problem. In this paper, we demonstrate that when training PINNs, NGD indeed enjoys faster convergence rate.
1.1 CONTRIBUTIONS
The main contributions of our work are summarized as follows:
For the PINNs, we simultaneously improve both the learning rate 1) of gradient descent and the requirement for the width m. The improvements rely on a new recursion for- mula for gradient descent. Specifically, our analysis yields a different step-size criterion 1) = ((1/Amma)) which empirically permits larger practical learning rates than the ((A)) requirement in Gao et al. (2023), see Remark 3.8. The requirement for the width m, i.c.
m = ②
(n1 (n²)"
can be improved to m = á (log( m1 (m)))=))))))))))))()))))) where indicates x
X
that some terms involving log(m) are omitted.
We present a framework for demonstrating the positive definiteness of Gram matrices for a variety of commonly used smooth activation functions, including the logistic function, softplus function, hyperbolic tangent function, and others. This conclusion is not only applicable to the PDE we have considered but can also be naturally extended to other forms of PDEs.
We provide the convergence results for natural gradient descent (NGD) in training over- parameterized two-layer PINNS with ReLU activation functions and smooth activation functions. Due to the distinct optimization dynamics of NGD compared to GD, the learning rate can be (1). Consequently, the convergence rate is independent of n and Ao. leading to faster convergence. Moreover, when the activation function is smooth, NGD achieve quadratic convergence rate.
1.2 RELATED WORKS
First-order optimizers. There mainly two approaches to studying the optimization of neural networks and understanding why first-order methods can find a global minimum. One approach is to analyze the optimization landscape, as demonstrated in Jin et al. (2017); Ge et al. (015). It has been shown that gradient descent can find a global minimum in polynomial time if the optimization landscape possesses certain favorable geometric properties. However, some unrealistic assumptions in these works make it challenging to generalize the findings to practical neural networks. Another approach to understand the optimization of neural networks is by analyzing the optimization dynamics of first-order methods. For the two-layer ReLU neural networks, as shown in Du et al. (2018), randomly initialized gradient descent converges to globally optimal solution at a linear rate, provided that the width m is sufficiently large and no two inputs are parallel. Later, these results were extended to deep fully-come fully-connected feedforward neural networks and ResNet with smooth activation functions Du et al. (2019). Results for both shallow and deep neural networks depend on the stability of the Gram matrices throughout the training process, which is crucial for convergence to the global minimum. In addition to regression and classification problems, Gao et al. (2023) demonstrated the convergence of the gradient descent for two-layer PINNs through a similar analysis of optimization dynamics. However, both Du et al. 2018) and Gao et al. (2023) require sufficiently small learning
2

---

Published as a conference paper at ICLR 2026
rate and large enough network width to achieve convergence. In this work, we conduct refined full-batch convergent analysis of the over- perameeeemtee PINN regime for GD and NGD, building upon Gao et al. (2023). There'r re contemporaneous work analysis concentrate stochastic Jin & Wu (2025) and non-overparrmettrrieed Nießen & Müller (2025) settings.
Second-order optimizers. Although second-order methods possess better convergence rate, they rarely used in training deep neural networks due to the prohibitive computational cost. As a variant of the Gauss-Newton method, natural gradient descent (NGD) is more efficient in practice. Meanwhile, as shown in Zhang et al. (2019) and Cai et al. (2019) NGD also enjoys faster convergence rate for the L regression problems compared gradient descent. Müller& & Zeinhofer (2023) proposed energy natural gradient descent for PINNs and deep Ritz method, demonstrating experimentally that this method yields solutions that are more accurate than those obtained through GD, Adam or BFGS. After observing the ill-conditioned loss landscape of PINNs, |Rathore (RathoeR et al. (2024) introduced a novel second-order optimizer, NysNewtonCG (NNCG), showing that NNCG can significantly improve the solution returned by Adam+L-BFGS Moreover, under the assumption that the PL'-condition holds, Rathore et al. (2024) demonstrated that the convergence rate of their algorithm is independent of the condition number, which is similar with result. However, although the PL.-connition holds for over weer-paammterrieed neural networks the context of regression problems Liu et al. (2022), it remains unclear whether this condition holds for PINNs. De Ryck et al. (2024) showed that operator-precondiiioonnn analysis establishes convergence for linearized PINN problems. In this paper, we provide the convergence analysis for NGD in training two-layer PINNs with ReLU3 activation functions or smooth activation functions, showing that it indeed converges at a faster rate.
1.3 NOTATIONS
We denote [n] 21, n2 for n E N. Given a set S, we denote the uniform distribution on by (S. We T(F} denote the indicator function of the event E. For positive functions (n) and f2(n), we use (n) = ((f(n)), f2()) = (n)) or fi (n) X (2( to represent (n) S C22((n), where C' universal constant C. A universal constant means constant independent of any variables. Throughout the paper, we use boldface to denote vectors. Given I1. ................................ .Td E R. we use (II,,... ra) or [_1,... , xrg] to denote a row vector with i-th component N;
for i E [d] and then (æ1, ... za)" E Rd is column vector. T
1.4 ORGANIZATION OF THIS PAPER
In Section 2, we provide the problem setup for training two-layer PINNs. We then present the improved convergence results of gradient descent for PINNs in Section 3. In Section 4, we analyze the convergence of natural gradient descent in training two-layer PINNs with ReLU3 activation functions and smooth activation functions. In Section 5, we conduct experiments to verify the theoretical results. The limitations are briefly discussed in Section 6 and we conclude in Section 7. All the detailed proofs and experiments are provided in the Appendix for readability and brevity.
PROBLEM SETUP
In this section, we consider the same setup as Gao et al. (2023), focusing on the PDE with the following form.
Ou Doo u(x) g(x), E (0) 0 [0, an,
Σ 02 (x) f(x), (0,T)) Ω
(x)
dx?
(1)
where C Rdd is an open and bounded domain, x (xn. 11,
x.a E Rd+1 and xn E [o, T]
is the time variable. In the following, assume that ||||| < 1 for E [0,T] l and f.g bounded continuous functions.
Moreover, consider two-layer neural network of the following form.
o( w. a) = 1
Σ
(w x)
(2)
3

---

Published as a conference paper at ICLR 2026
where (wi w " Rm(d+)), == (a1,... am)) r " and for r [m], w, R
is the weight vector of the first layer, or is the output weight and o(.) is the ReLU activation function.
Here, (x* 1) TR++2 is the augmented vector from and in the following, we write x for
for brevity.
Similar to that for the L regression problems, we initialize the first layer vector w,(0) ~ .o I). output weight a, Unif([ [-1, for [m]] and fix the output weights. In the framework of PINNs, given training samples {xp 1 and {y, ور that are from interior and boundary respectively,
we denote sp(w) and (w)) by
1 Vni
do Ozn
d
02
(n(w)
p w)
ppi
w) -((xp))
(3)
dr
and
hy(w)= =
(0(ys) w) (yy))
(4)
f02
Then the empirical loss function can be written as
1
L(w)= = ([]((w)) 112 + (/(ww))))
(5)
2
where s(w) = (w), Sn, (w)) )² R. and (ww) = (h1(w), Ang (w)) E R
The gradient descent updates the hidden weights by the following formulations:
w, l + 1) w,(k) - aL(w(k))
(6)
Ow,
for all [m] and & N, where 1j 0 is the learning rate. The Gram matrix H is defined as (ww) = JJJ", where
@s1(ww) dw
Osn, (w) ow
(1)(w) Ow
Öhng (w))
J
dw
IMPROVED RESULTS OF GD FOR TWO-LAYER PINNS
To simplify the analysis, we make the following assumptions on the training data.
Assumption 3.1. For p E [mi] and j E [m2]. |1>p12 I v2, b S /2, where all inputs have been augmented.
Assumption 3.2. No two samples [ampp== U y د are parallel. This is guaranteed because augmenting se with (æ: 1) ensures all samples are distinct.
Under Assumption 3.2, Lemma 3.3 Gao al. (2023) implies that the Gram matrix H Ew H (w)) is strictly positive definite and we let λ Amin (H.).. Similar to the case of regression problem in Du et al. (2018), H plays an important role the optimization process. Specifically, under over- paramterrizzaioonnn and random initialization,
have two facts that (1) at
initialization (0) - O(L/vm), The following two lemmas can be used to verify these two facts, which are crucial in the convergence analysis.
H l| = OL1 /m) and (2) for any itcration E N. (k) H (O)]l2 =
Lemma 3.3. Ifm = Ω
d4 X
log ma))t))) we have that with probability at least 5, | [H
H < Ag
4 and Amin (0)) 2 Ao-
Remark 3.4. Under the premise of deriving the same conclusion as our Lemma 3.3, the Lemma 3.5 in Gao al. (023) requires that Ω (m)+ (log())))) where some terms involving log(m) are omitte. In contrast, one hand, our conclusion is independent up logarithmic factors in n1 + 722, and on the other hand, our conclusion exhibits a clear dependence on d. Moreover, the method in Gao et al. (2023) involves truncating the Gaussian distribution and then applying Hoeffding's inequality, which quite complicated. In contrast, we utilize the concentration inequality for weibul random variables, which serves simple framework for this class of problems.
[ming]"
4

---

Published as a conference paper at ICLR 2026
Lemma 3.5. Let R E (0, .11. if w1 (0).... (0) are i.i.d. generated from (o0, I), then with probability at least 1 1 - -n the following holds. For any set of weight vectors W1, ... Wm E Re+1 that satisfy ||w, w, (,)ll for any E [mm) then
( H(O) < CMMR,
(8)
where M 2(d log (2m 2) and C is universal constant.
Remark 3.6. The Lemma 3.6 in Gao et al.
(023) shows that when l|,,, (,)(2 S R (lggg m)² holds for all E [m], then (u) H())l & Aa
And
4. In contrast, our Lemma
Ap 3.5 only requires - (log
to reach same result.
For the L regression problem, shown in Du al. (2018), the convergence of gradient descent requires that the learning rate 7/ - ((oonn") where vi is the sample size of the regression problem. It evident that this requirement on the learning rate is difficult satisfy in practical scenarios, since An is unknown and n² is too large For PINNs. Gao et al. (2023) follows the methodology of Du al. (2018), thus inheriting similarly stringent requirements on the learning rate. By investigating new decomposition method for the residual, we arrive at our main result.
Theorem 3.7. Under Assumption 3.1 and Assumption 3.2, if we set the number of hidden nodes
d
mmd"
n21 + 122
m =
log
log
and the learning rate 7 O initialization, the gradient descent algorithm satisfies
TH-TA then with probability least 1 - 5 over the random
17A0 / L(0)
L(k)= -
for all ls E N.
Remark 3.8. It may be confusing that Gao al. 2023) has used the method in Du al. (2018), yet it only requires = O(Ao). Actually, it is because that the loss function of PINN has been normalized. If we let n1 n2 and H be the Gram matrix induced by unnormalized loss function of PINN, then Amin [ 00 Amin H0))/n leading to the convergence rate similar to that of regression problem. At this point, due to the normalization of loss function, || H ||2 Amaz ) can be bounded by the trace = which is an explicit constant that independent of the sample size 1, 712. As Amin depends the sample size, it is expected that our 7 ((1/Aman) \moa)) an improvement over = ((Amin) in Gao et al. (2023). A practical computation for 1D Poisson equation is Amin 3.47 x 10 and 1/Amax 17 x 104) 5.78 x 10 suggesting that our analysis indeed improves the learning rate requirements.
4 CONVERGENCE OF NGD FOR TWO-LAYER PINNS
Although we have improved the learning rate of gradient descent for PINNs, it may still be necessary the learning rates to be sufficiently small. Because, although tr tr(Ho) is an explicit constant, depends on the form of the PDE. However, the loss function of PINNs has a much worse conditioning due the appearance of the PDE operator. So the ill-conditioning occurs when we move from regression to PINNs, brings strict restrictions for the learning rate of gradient descent for PINNs. This a central motivation for second-order and natural gradient methods. Moreover, the convergence rate 1 iAo also depends on λo. which depends on the sample size and may extremely small. Zhang et al. 2019) and Cai et al. (2019) have provided the convergence results for natural gradient descent (NGD) training over- wverparramteerized two-layer neural networks for L regression problems. They showed that the maximal learning rate can be ((1) and the convergence rate is independent of λo which result a faster convergence rate. However, the situation in PINNs is significantly different from regression due to the presence of derivative terms from the partial differential equations, which complicates the analysis. Müller & Zeinhofer (2023) and Guzzán- et al. (2025) studied the energy natural gradient descent (ENGD) for PINNS with practical Woodbury, momentum,
5

---

Published as a conference paper at ICLR 2026
and randomization techniques, demonstrated highly accurate solutions empirically. A theoretical convergence analysis, however, has not yet been established in these works. We note that the NGD in Zhang al. (2019) and the ENGD in Müller & Zeinhofer (2023) coincide up to the choice of the Moore-Penrose pseudoinverse or the use of the Woodbury matrix identity; therefore, we do not distinguish between them in this work. In the section, we conduct the convergence analysis of NGD) for PINNs and demonstrate that it results in faster convergence rate for PINNs compared to gradient descent.
In this section, we consider the same setup as described in Section 2. During the training process, we fix the output weight and update the hidden weights via NGD. The optimization objective the empirical loss function presented in (5), which is defined as follows:
L(w) 1 (I|s(w)) [[((w)
(10)
The NGD gives the following update rule:
wik + 1) = wkk J(=) (J())
-1
(8(k) h(k))
(11)
where
J(k) = (J)( +n, E R (na)mmd))
is the Jacobian matrix for the whole dataset and v) > 0 is the learning rate. Specifically, for p E [Mi].
( (Dsp(k) T
@sp(k) Twn
27
(k)
E IR Ixm((+))
(12)
and for j E [12]]
( (k)
T
(k)
J [k)+s(k)
E R Lxm(d+2)
(13)
dee
Own
We note that Zhang et al. (019 and Cai et al. (2019) have independently and concurrently established the convergence of NGD in the context of regression problems. The difference lies in the fact that Zhang et al. (019) focused on ReLU activation functions, whereas Cai let al. (211)) considered smooth activation functions and consistently set the learning rate to 1. Here, following Zhang et al. (2019). we refer to this approach as NGD. In Cai et al. (201)), the authors derived this method based on NTK kernel regression and termed it the Gram-CGaus-NNwwan (GGN) method. The extension of NGD convergence from regression to PINNs is challenging because of the complexity of the PDE residual loss.
Remark 4.1.
Remark 4.2. The classical wok - (J))) J(R))T
method Bonfanti et al. 2024) is given by w + 1)
((k))
Although this formula looks different from the NGD
(n(k)
update (11), the two coincide when 17 = 1 at the level of the Moore-Penrose pseudoinverse: J(b)+ = (J))) J(k)T JKk)" (J(+) -1 However, this equivalence is only algebraic. In practice the two updates behave differently because J (k) E (m, (n2)+ is highly rectangular and never invertible strictly, and different pseudoinverse representations apply in row-dependent or column-dependent cases. The computational cost are also different, as pointed in Guuzmnn-Cordero eet al. (2025) with Woodbury's 's Identity. The NGD's formula (11) in this work, originally adopted from Zhang et al. (2019), also coincides with the energy natural gradient descent (ENGD) proposed in Müller & Zeinhofer (2023); [Guzmán-Cordena
al. (0225) once the Moore- Penrose inverse or Woodbury identity is applied. A crucial distinction arises in the over-parameterized regime. The Gram matrix J ((k)T"(() E ( becomes extremely high-dimensional and typically singular as m grows, while the NGD (K)) e (R (n2)))))x(()))))) won't 't. This difference is key for both practical scalability and numerical stability.
For the activation function of the two-layer neural network
o( w. a) = a x)
1 Vm
Σ
(14)

---

Published as a conference paper at ICLR 2026
consider settings where ((-) is either the ReLU activation function smooth activation function satisfying the following assumption.
Assumption 4.3. There exists a constant c > 0 such that suPsER (=)] S c and for any 2, E R,
1or o (2 )l S c|2 - 211.
(15)
where & E {0, 1,2,3.. Moreover, (.) is analytic and is not polynomial function. We also assume that for any positive integer 2 2, lim a(1) (x)/ó(x) Cn # 0, where the function Φ(-) needs
(bx) =0 +00 6(x)
lim ((x) = 0, lim x
holding for any constant >
Lemma 4.4. If no two samples in p U
و are parallel, then the Gram matrix H is strictly positive definite for activation functioos that satisfy Assumption 4.3, i.e. do == Aman (H²)
Remark 4.5. Assumption 4.3 holds for various commonly used activation functions, including
logistic function (=)) = 1/(1 + (with ((z) = e--), softplus function (=) = log(1 + + (with
6(2) = ") hyperbolic tangent function o (z)) (e* e *)/(e* +)= (with ((z) = 2²), swish
function a (z) = z/(1 + e e-) (with o(z) = ze =e=²) +) and others.
Unlike the approach for gradient descent, Zhang et al. (2019) focus on the change of the Jacobian matrix for NGD rather than the Gram matrix. Roughly speaking, they show that when ||wr 2 (0)ll2 small, then [| (w) - (() is also proportionately small. However, this approach is not applicable to PINNs, because the loss function involves derivatives. Roughly speaking, the stability considered in Zhang et al. 2019) is more global in nature, whereas ours is local. In fact, the PINN loss includes first- and second-order derivatives of the neural output (see Eq. (3)-44), so each Jacobian block 0sp/Ou, and /h,/Ou, contains higher-order derivatives of the activation and of the weights. Consequently. small perturbation in weights may cause large variations in the derivatives, violating the Lipschitz-type condition required by Zhang et al. (2019). Since the subsequent conclusions require the boundedness of local weights, we do not
this stability. Moreover, from the Theorem 1 in Zhang et al. (2019), we can see that this stability imposes additional constraints on the learning rate. Therefore, we instead focus on the stability of J (w) with respect each individual weight vector w, in the following Lemma, which provides more targeted approach.
Lemma 4.6. Let R E (0, 1}, (0), ... um (0) are i..dd. generated ( I), then with probability at least 1 - P(,, m. R) the following holds. For any set of weight vectors Wi,' Wm E ddd+2 that satisfy for any r E [m]] l/w, - w,(() ||2 < R, then
(I) when (-) is the ReLU activation function, we have that
[ - -(0))l2 S CMMR, where C. is a universal constant, M 2(d + 2) log(2 (2m(2md 2)/5) and
(16)
P(8, m, R) = 6 + mme =m;;
(17)
(2) when (-) is the smooth activation function satisfies Assumption 4.3, we have that [J(w) J(0)l22 CdR
(18)
for m 2 logg((0 (1/5), where C is a universal constant and P(š,m, R) = 6.
With the stability of Jacobian matrix, we can derive the following convergence results.
Theorem 4.7. Let L(k)- = Lun(k)), then the following conclusions hold.
(I) When a a-) is the R activation function, under Assumption 3.2, we set
1 de
md
+ n22
= Ω
(1 n)² log
log
X1
8
and 1 E (0, 1), then with probability at least 1 over the random initialization for all k E N
T(k)) (1 .1)* I(00)
(19)
7

---

Published as a conference paper at ICLR 2026
(2) When (.) is the smooth activation function satisfies Assumption 4.3, under Assumption 3.2, set
de
md
Th1 + n22
m = Ω
log"
log
1
X
5
and 17 E (0, 1), then with probability least 1 over the random initialization for all k E N
L(k) S (1 n1)* L(O).
(20)
In Theorem 4.7, the requirements of m with ReLU and smooth activation functions exhibit different dependencies on Äg and d. The discrepancy is primarily due to the distinct formulations presented in (16) and (18) of Lemma 4.5.
Remark 4.8. We first compare our results with those of NGD for I,2 regression problems. Given that the convergence results are the same, focus shifts to examining the necessary conditions for the width m. As demonstrated in Zhang et al. (2019) and Cai et al. (2019), required that for ReLU activation function and m
d login/6) for smooth activation function. Clearly, our result has worse dependence d, which is inevitable due to the involvement of derivatives the loss function. Moreover, our requirement for m appears to be almost independent of n. primarily because our loss function has been normalized. With smooth activation functions, in addition to the dependence on d, Theorem 4.7 (2) only requires that m M(Ao However, Cai et al. (2019) demands a more stringent condition, requiring that m = (2
m -
max
A
Comparing with our results in Section 3, the requirement for in Theorem 4.7 (1) is the same as in Theorem 3.8, when we make 1) less close to 1. On the other hand, since 7) = ((1) and the convergence rate only depends 17, NGD can lead to faster convergence than GD.
Note that as approaches 1, the width m tends to infinity, thus, the convergence results in Theorem 4.7 become vacuous. In fact, when though m is finite, provided that (-) satisfies Assumption 4.3.
1. NGD enjoy second-order convergence rate even
Corollary 4.9. Under Assumption 3.2 and Assumption 4.3, set 1 1 and
Ω log²
d
(mr)
(md)
log
Th + 122
8
then with probability least 1 - 8, we have
( ( +1
CB
h(t++1) 11)+
S
V m
((t))
for all t E N, where C is universal constant and B / 2(d 2) log (2 2)/8) + Moreover, we can get a second order convergence result for regression problems with smooth activation functions functtin as follows.
3/2
u(l 1)112
ly (1)1)
Vm
Instead of inducing the convergence rate of the empirical loss function, as shown in Condition 1, we perform induction on the movements of the hidden weights as follows. Condition 1. At the t-th iteration, have [w, B and
C L(O)
w, (t) w. (0)l22 S
vmdo == R
for all r [m]] where is universal constant and B 2(d + log (d+(2+
+ 1.
With Condition I, we can directly derive the following convergence rate of the empirical loss function. Corollary 4.10. If Condiiiion 1 holds for = 0,-.. kk and H' < R and R²" < - No. then L(t) & (1 (1-1)) 1)) L(0),
holds for = 0,... ,k, where R is the constant in Lemma 4.5 and R" CM / is in (16) when the 3 activation function, R" CdR is in (18) when satisfies Assumption 4.3.
8

---

Published as a conference paper at ICLR 2026
5 EXPERIMENTAL RESULTS
We conduct a comparative evaluation of the NGD against existing optimizers for PINN training with respect to accuracy, computational efficiency, and alignment with theoretical analysis. The experimental implementations are listed in the Appendix A.
Table I: Relative -ererror of Different Optimizers.
SGD
Adam
L-BFGS
NGD
1D Poisson
1.28e-01 +4,31e-02
6.46e-02 + 144000
2.63e-04
+8.05-05
1.67e-05
9.007-06
2D Poisson
1.45e-01 +7,340002
5.32e-03
+9.76004
3.17e-03
+8.60-04
1.12e-04
+6.90.05
ID Heat
5,43e-01 19,98-02
6.91e-03
1.16-03
4.98c-03
1.83300
3.42c-04
7.522.05
2D Helmholtz
8,48e+00
26.774000
1.06e+00
1.l1e:0
3,35e+00
1.946400
6.67e-03
1,86003
T0D Poisson
1.35e-02 8.177-00
3.15e-03 +8.36004
9.91e-04
1,47004
Comparison to Existing Optimizers. We report the relative 1,2-error of the NGD optimizer the commonly used first order optimizers (the SGD optimizer, the Adam optimizer) and second order optimizer (the L-BFGS optimizer) in Table Here 'nan' means the training loss becomes infinity. We see that NGD performs best on all five equations.
Learning Rate Study. We report behavior of convergence at different learning rates, showing the strong robustness of the NGD method to hyperparameter selection. Table Table[ demonstrates that, unlike SGD and Adam which demand small learning rates for convergence, the NGD maintains stable convergence across wide range of learning rates without notable accuracy deterioration. This characteristic, which markedly outperforms conventional optimization approaches, clearly illustrates the strong robustness of the NGD method to hyperparameter selection.
Table 2: Relative L error Comparison Across Different Learning Rates 77.
learning rate 1)
1.0
0.5
0.1
0.05
0.01
0.005
0.001
SGD
nan
nan
nan
nan
1.19e-02
6.91c-02
7.36e-02
Adam
1.0le+00
1.00e+00
1.00c+00
1.0lc+00
1.64c-02
3.25e-02
1.49c-02
NGD
1.97e-03
1.18c-03
3.24e-04
1.87e-04
1.12e-04
1.22e-04
1.68e-04
Network Width Study. A comparative analysis of the model performance performed with progressively increasing network widths. Table demonsteesssssstes that increasing network width leads significant accuracy improvements. This trend validates that wider architectures exhibit enhanced function approximation capabilities.
Table Relative -e-rrr Comparison Across Different Network Width m for NGD,
20
40
80
160
320
640
1280
2560
error
1.59-03
7.21e-04
5.18e-04
3.8e-04
3.08e-04
2.76e-04
1.78e-04
7.05e-05
Fast Convergence Study, We report the training loss convergence results for different optimizers. We train SGD and Adam for 10000/20000 epochs with learning rate le-3, and the NGD for 100/200 epochs with learning rate 0. Figure Dempirecal demonstrates that the NGD converges much faster than commonly used SGD and Adam optimizers, which consistent with our theoretical analysis equation (9) in Theorem 3.7 and equation (20) in Theorem 4.7.
Empirical Convergence Rates Study. We continue to report the empirical convergence rates of the NGD in different equations. We compare the empirical training loss curves of the NGD when ? = 0.1 with the theoretical linear rates in our main Theorems 4.7. The theoretical decay follows L(k) R= C(1 r and the fitted experimental decay is L(k) O(k-1.55 for ID Poisson equation, O 1.92) for ID Heat equation and O 11.13) for 2D Poisson equation in Figure 2 For the heat
equation, convergence initially exceeds the predicted rate and later slows markedly. This is consistent

---

Published as a conference paper at ICLR 2026
877
24
877
Figure 1: Training Loss Decay Comparison for ID Poisson (left), ID Heat (middle) and 2D Poisson (right) Equations.
with known NTK decay and multi-phase behaviors in PINNs. Generally, the empirical loss of NGD roughly follows the predicted linear regime in early iterations, before entering slower phase usually observed across all optimizers.
thandted boow
30"
30
377
337
4p075
403070
spccss
Figure 2: The Experimantal Training Loss Decay and Theoretical Decay (Theorem 4.7) for 1D Poisson (left), 1D Heat (middle) and 2D Poisson (right) Equations.
LIMITATIONS
The computational cost NGD mainly on the JJ-JT)-1 with the Jacobian matrix J is of size p. where = 71 + is the training data size and m(d + 2) is the number of trainable parameters. NGD will quite expensive for large amount of training data. a result, several cost-effective variants have been proposed, such as K-FAC|Martens & Grosse (2015); Dangel et al.) (2024) 2025), ENGD Müller Zeinhofer (2023);
Cordero et al. (2025) and mini-batch NGD. We only proved the convergence results for the full-batch NGD in this paper, and it would be interesting to investigate the convergence of these methods for PINNs future works. On the other hand, while the aver-parameterized assumption enables the use of NTK stability for proving global convergence, the practical guarantee for arbitrary sampled projected gradient descent without assumption on the network size Nießenn & Müller (2025) address different framework, and the NGD analysis without over- parameterized assumption represent interesting complementary direction.
7 CONCLUSION AND OUTLOOK
In this paper, we have improved the conditions required for the convergence of gradient descent for PINNs, showing that gradient descent actually achieves better convergence rate. Furthermore, we demonstrate that natural gradient descent can find the global optima of two-layer PINNs with ReLU3 or smooth activation functions for class of second-order linear PDEs. Compared to gradient descent, natural gradient descent exhibits faster convergence rate and its maximal learning rate is O(1). In conclusion, the NGD offers three key advantages: 1) more relaxed learning rate requirements; 2) faster convergence rates independent of Xo; 3) superior empirical performance. Additionally, extending the convergence analysis to deep neural networks, stochastic version of NGD, and studying the generalization error of trained PINNs are important directions for future research.
10

---

Published as a conference paper at ICLR 2026
ETHICS STATEMENT
We acknowledge that all authors have read and commit to adhering to the ICLR Code of Ethics.
REPRODUCIBILITY STATEMENT
To facilitate reproducibility, have taken the following steps: (1) Source code and configuration
files for all key experiments are provided as supplementary material. (2) All theoretical claims accompanied by full proofs (in the Appendix) and assumptions are clearly stated. (3) All
hyperparameters used, neural network architecture details are provided either in the main text or in the Appendix.
REFERENCES
Zeyuan Allen-Zhu, Yuanzhi Li, and Yingyu Liang. Learning and generalization in overparameterized neural networks, going beyond two layers. Advances in neural information processing systems, 32, 2019a.
Zeyuan Allen-Zhu, Yuanzhi Li, and Zhao Song. A convergence theory for deep learning via parameterization. International conference on machine learning, pp. 242-252. PMLR, 2019b.
Sanjeev Arora, Simon Du, Wei Hu, Zhiyuan Li, and Ruosong Wang. Fine-grained analysis of optimization and generalization for overparameterized two-layer neural networks. In International Conference Machine Learning, pp. 322 332. PMLR, 2019.
Andrea Bonfanti, Giuseppe Bruno, and Cristina Cipriani. The challenges of the nonlinear regime for physics- -nnssies-nformee neural networks, Advances in Neural Information Processing Systems, 37: 41852-41881, 4181, 2024.
Tianle Cai, Ruiqi Gao, Jikai Hou, Siyu Chen, Dong Wang. Di He, Zhihua Zhang, and Liwei Wang. Gram-gause newton gewss- method: Learning overparameterized neural networks for regression problems. arXiv preprini orX5v: /665 2019.
Yuan Cao and Quanquan Gu. Generalization bounds of stochastic gradient descent for wide and deep neural networks. Advances in neural information processing systems, 32, 2019.
Felix Dangel, Johannes Müller, and Marius Zeinhofer. Kronecker- factored approximate curvature for physics- informmed neural networks, Advances in Neural Information Processing Systems, 37: 34582-34636, 2024.
Felix Dangel, Bálint Mucsányi, Tobias Weber, and Runa Eschenhagen. Kronecker-factoree approxi- mate curvature (kfac) from scratch. arXiv preprini arXiv:2507.05/27, 2025.
Tim De Ryck, Florent Bonnet, Siddhartha Mishra, and Emmanuel de Bezenac. An operator precondi- tioning perspective on training in physics-informed machine learning. In The Twelfth International Conference Learning Representations, 2024.
Jacob Devlin, Ming-Wer Wei Chang, Kenton Lee, and Kristina Toutanova. Bert: Pre-training of deep bidirectional transformers for language understanding. arXiv preprint a 048000 2018.
Simon Du, Jason Lee, Haochuan Li, Liwei Wang, and Xiyu Zhai Gradient descent finds global minima of deep neural networks. In International conference machine learning, pp. 1675-1688 PMLR, 2019.
Simon S Du, Xiyu Zhai, Barnabas Poczos, and Aarti Singh. Gradient descent provably optimizes over-parameterized neural networks. arXiv preprint arXiv: 110.00054 2018.
Yihang Gao, Yiqi Gu, and Michael Ng. Gradient descent finds the global optima of two-layer physics- informed neural networks. In International Conference on Machine Learning, pp. 10676-10707. PMLR, 2023.
11

---

Published as a conference paper at ICLR 2026
Rong Ge, Furong Huang. Chi Jin, and Yang Yuan. Escaping from saddle points online stochastic gradient for tensor decomposition. Conference on learning theory, pp. 197-842. PMLR, 2015.
Evarist Giné and Richard Nickl. Mathematical foundations of infinite-dimensionnl statistical models.
foundaation Cambridge university press, 2021.
Andrés Guzmán-Corden Felix Dangel, Gil Goldshlager, and Marius Zeinhofer. Improving energy natural gradient descent through woodbury, momentum, and randomization. Advances in Neural Information Processing Systems, 2025.
Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recognition. In Proceedings of the IEEE conference on computer vision and pattern recognition, pp- 770-778, 2016.
Arthur Jacot, Franck Gabriel, and Clément Hongler. Neural tangent kernel: Convergence and generalization in neural networks. Advances in neural information processing systems, 31, 2018.
Bangti Jin and Longjun Wu, Convergence of stochastic gradient methods for wide two-layer physics- informed neural networks. arXiv preprint wrXiv:2508.21577 2025.
Chi Jin, Rong Ge, Praneeth Netrapalli, Sham M Kakade, and Michael Jordan. How to escape saddle points efficiently. In International conference on machine learning, pp. 1724-1732. PMLR, 2017.
Arun Kumar Kuchibhotla and Abhishek Chakrabortty. Moving beyond sub-gaussianity high- dimensional statisties: Applications in covariance estimation and linear regression. Information and Inference: A Journal of the IMA, 11(4):1399-145, 2022.
Yuanzhi Li and Yingyu Liang. Learning overparameterized neural networks via stochastic gradient descent structured data. Advances in neural information processing systems, 31, 2018.
Chaoyue Liu, Libin Zhu, and Mikhail Belkin. Loss landscapes and optimization in weer-parameterized non-linear systems and neural networks. Applied and Computational Harmonic Analysis, 59: 85-116, 2022.
James Martens and Roger Grosse. Optimizing neural networks with kronecker-factored approximate curvature. In International conference on machine learning, pp. 2408-2417. PMLR, 2015.
Johannes Müller and Marius Zeinhofer. Achieving high accuracy with pinns via energy natural gradient descent. In International Conference on Machine Learning, pp. 2023.
PMLR,
Jonas NieBen and Johannes Müller. Non-asymptotic analysis of projected gradient descent for physics-informed neural networks. arXiv preprint arXiv:2505.07311. 2025.
Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics, 378:686-707, 2019.
Pratik Rathore, Weimu Lei, Zachary Frangella, Lu Lu, and Madeleine Udell. Challenges in training pinns: loss landscape perspective. In Forty-first International Conference on Machine Learning. 2024.
Jonathan W Siegel, Qingguo Hong, Xianlin Jin, Wenrui Hao, and Jinchao Xu. Greedy training algorithms for neural networks and applications to pdes. Journal of Computational Physics, 484: 112084, 2023.
David Silver, Aja Huang, Chris J Maddison, Arthur Guez, Laurent Sifre, George Van Den Driessche, Julian Schrittwieser, loannis Antonoglou, Veda Panneershelvam, Mare Lanctot, et al. Mastering the game of go with deep neural networks and tree search. nature, 529(7587):484-489, 2016.
Bing Yuet al. The deep ritz method: deep learning-based numerical algorithm for solving variational problems. Communications Mathematics and Statistics, 6(1):1-12, 2018.
Yaohua Zang, Gang Bao, Xiaojing Ye, and Haomin Zhou. Weak adversarial networks for high- dimensional partial differential equations. Journal of Computational Physics, 41 1:109409, 2020.
12

---

Published as a conference paper at ICLR 2026
Guodong Zhang, James Martens, and Roger B Grosse. Fast convergence of natural gradient descent for over parameterreee neural networks. Advances Neural Information Processing Systems, 32, 2019.
Difan Zou, Yuan Cao, Dongruo Zhou, and Quanquan Gu. Gradient descent optimizes over- parameterized deep relu networks. Machine learning. 109:467-492, 2020.
A EXPERIMENTAL IMPLEMENTATION
In this section, we provide several examples to demonstrate the superiority of the natural gradient descent(NGD) approach. The configurations used in these examples are listed in Table We report the relative -er-r of the NGD optimizer to the commonly used first order optimizers (the SGD optimizer, the Adam optimizer) and second order optimizer (the L-BFGS optimizer) in Table T The relative -eror is defined as follows:
E relative L error Li=1 [((x;) Mre (x:))
(21)
[xi)r(xx)
where denotes the predicted solution and Uyef represents the reference solution. To show the generalization ability of NGD, should note that the testing collocation points J are different from the training samples {
) and {y, 1-=1-
Table 4: Configurations of Different Equations.
N₂
No
batch size
hidden layers
hidden neurons
activation function
ID Poisson
500
2
100
128
tanh(-)
2D Poisson
1,000
200
100
128
tanh()
ID Heat
1,000
200
100
128
tanh()
2D Helmholtz
1,000
200
100
128
tanh
10D Poisson
10,000
1,000
100
128
tanh()
A.1 ID POISSON EQUATION
First, we begin with a toy example of the ID Poisson equation to display the performance of the NGD method. The equation is defined in the domain f [0, n.,
(x(x) fx), x E Ω u(x)= = 0. an.
(22)
The true solution is set as u(x) si) which allows us to derive the corresponding force term f() = sin. We randomly sample N, 500 points in the domain Q. For the neural network architecture, we employ single hidden layer model with 128 units and tanh(-) activation functions across all computations. The NGD optimizer trained for 100 epochs, while the LBFGS optimizer is run for epoch with maximum of 500 iterations per epoch. All other optimizers are run for 10, 000 epochs for comprehensive comparison. The relative 2 is 1.67e 05 for the NGD optimizer. Figure 3 shows the predicted solution for the ID Poisson equation alongside the reference solution. The prediction
in excellent agreement with the reference solution, highlighting the depicts superior performance of the NGD method. Figure deepict the loss decay during the training process, we can see that the NGD method achieves a quite small loss at the very beginning.
A.2 2D POISSON EQUATION
We consider a 2D Poisson equation in the domain Ω = [0, x [0, 1).
8 2
82
f(x, (x. v) E 2
Oy
(23)
u( = 0, (x. v) E 0n.
13

---

Published as a conference paper at ICLR 2026
1.0
0.8
0.6
0.4
0.2
Reference solution
0.0
Predcited solution
0.0
0.5
1.0
1.5
2.0
2.5
3.0
Figure 3: Reference Solution and Predicted Solution for the ID Poisson Equation.
The true solution given by u(x,y) sin( (y)), and the force term f(x, y) 2x sin(Tx) sin(Tr) is consequently derived.
=
We sample N₂ = 200 random points on the boundary 80 and N, = 1, 000 random points within domain Ω. We employ single hidden layer model with 128 units and tanh(-) activation functions across all computations. We run the NGD method for 200 epochs, while the L-BFGS method is trained for epoch with maximum of 5, 000 iterations per epoch. All other optimization methods are trained for 20, 000 epochs. The resulting relative T/-error 1.12p - 04. Figure 5 illustrates the prediction of the 2D Poisson equation, along with the exact solution and the absolute error between them. It is clear that the predicted solution closely matches the reference solution, further demonstrating the superior performance of the NGD method. Figure 6 shows the loss decay during training, demonstrating that the NGD method converges significantly faster than other optimization methods.
A.3 ID HEAT EQUATION
We consider the ID heat equation
Ou(t.x) )u(((.)) St
(1, E [0,
Drd
u(0,x) = sin TT))
E [0.1]
(24)
u(t,z) =0 O, (t,
( [0,1] (0.1).
The reference solution is analytically defined by u(t, exp( sin(TTr) We generate N,, 200 random sampling points for the boundary and initial conditions and N, = 1, 000 random points in the domain Ω to evaluate the PDE residual. The neural network used for all computations consists of hidden layer, each containing 128 neurons with tanh(-) activation functions. To train the model, the NGD method for 200 epochs and the L-BFGS method for 1 epoch with maximum of 5, 000 iterations per epoch, and other optimizers are trained for 10, 000 epochs. The resulting relative err-r is 3.42e 04. Figure provides visual comparison between the predicted and
14

---

Published as a conference paper at ICLR 2026
NGD
100
SGD
Adam
10-2
loss
10-4
10-6
10
0
2000
4000
6000
8000
10000
epochs
Figure 4: Loss Decay for the 1D) Poisson Equation.
08
00
64
05
C4
66
Figure 5: NGD Prediction and Analysis for the 2D Poisson Equation.
15

---

Published as a conference paper at ICLR 2026
102
NGD
SGD
101
Adam
100
10
loss
10
10-3
10
10-5
10-6
0
2500
5000
7500
10000 epochs
12500
15000
17500
20000
Figure 6: Loss Decay for the 2D Poisson Equation.
exact solutions for the ID heat equation, along with the corresponding absolute error distribution. The high degree of accuracy in the predicted solution demonstrates the effectiveness of the NGD method, showing its ability to capture the solution with remarkable precision. Figure shows the loss curve over the course of training for the ID heat equation. Notably, the NGD method rapidly reduces the loss, reaching low value in the training process, demonstrating its efficiency in optimization.
Buc: SAV731
Abodute
94
94
02
Figure 7: NGD Prediction and Analysis for the ID Heat Equation.
A.4 2D HELMHOLTZ EQUATION
We deal with the 2D helmholtz equation the domain Ω [0, [0, given by
8 + बي By u(I,y) = O. (x,y) E OΩ.
8hv + 12u( u) = [(x, y), (x. v) E Ω
(25)
16

---

Published as a conference paper at ICLR 2026
10
10
10-3
loss
10
10-5
10-6
NGD
SGD
Adam
10
0
2000
4000
6000
8000
10000
epochs
Figure 8: Loss Decay for the ID Heat Equation.
The reference solution for k: is u(æ; v) (rx))) sin(4ry), and the force term f(x y) can be easily computed. To evaluate the performance of the NGD approach on the 2D Helmholtz equation, generate N₁- 200 random boundary points an and N, = 000 random points inside the domain Ω. The neural network employed consists of 1 hidden layers with 128 neurons utilizing t-nh(-) activation functions. Training carried out for 200 epochs using the NGD method and epoch with maximum of 5, 000 iterations for L-BFGS. All other optimizers are run for 20, 000 epochs for comparison. The computed relative err-r 6.67e 03, which is orders of magnitude 9 illustrates the predicted solution along with
layer,
lower than those of the remaining optimizers. Figure the exact reference solution and the absolute error distribution. The results indicate that the NGD method effectively captures the oscillatory nature of the Helmholtz equation, achieving high level of accuracy. Figure I0] shows the evolution of the loss function during training for the 2D Helmholtz equation. In particular, the NGD method demonstrates rapid convergence, achieving low loss value at the end of the training process.
A.5
1OD POISSON EQUATION
We conduct experiments to show that the NGD can also perform better than SGD, Adam, and L-BFGS for high-dimensional PDEs, despite that all optimizers become more challenging as the dimensionality of PDEs increases. As an example in higher dimensions, we consider again the Poisson equation 10 spatial dimensions
Au = f(x), x xEΩ- [0, 1] u(x)
(26)
Ele sin(TTTx), x E 09.
We use the manufactured solution
10
u
10
IR I + Σ
sin (
(27)
17

---

Published as a conference paper at ICLR 2026
Pedirtan
50
C034
-09
02
Figure 9: NGD Prediction and Analysis for the 2D Helmholtz Equation.
104
NGD
SGD
103
Adam
102
101
loss
100
10
10
10-3
10-4
0
2500
5000
7500
10000
12500
15000
17500
20000
epochs
Figure 10: Loss Decay for the 2D Helmholtz Equation.
18

---

Published as a conference paper at ICLR 2026
hence f = = We sample N₂- 1, 0.00 random points the boundary ōn and N 10, 000 random points within the domain Q. We employ a single hidden layer model with 128 units and tanh(-) activation functions across all computations We run the NGD method for 200 epochs. while the L.-BFGS method trained for epoch with a maximum of 5, 000 iterations per epoch. All other optimization methods trained for 20, 000 epochs. The resulting relative error 9.91e 04. It is clear that the predicted solution closely matches the reference solution, further demonstrating the superior performance of the NGD method.
A.6 TRAINING EFFICIENCY COMPARISON
For training efficiency comparison among different optimizers, we present the computational time, memory usage, and error rates for both the 21) Poisson equation and the 10I) Poisson equation. Table and T demonstrate that training time and memory requirements increase for all four optimizers as the problem dimension grows. Despite this, the NGD method still achieves the lowest error while maintaining comparable computational overhead.
Table
Training efficiency comparison for 2D Poisson equation.
Optimizers
Ir
Epochs
Training efficiency
Training time
Max memory
Rel. 1.2 error
SGD
0.001
20,000
0.047 s/epoch
15min49s
14.62 MB
1.45e-01
Adam
0.001
20,000
0.054 s/epoch
18min2s
14.75 MB
5.32e-03
L-BFGS
200
0.51 s/epoch
Imin41s
41.53 MB
3.17e-03
NGD
0.1
200
3.67 s/epoch
12min
14.75 MB
1.12e-04
Table 6: Training efficiency comparison for 10D Poisson equation.
Optimizers
Ir
Epochs
Training efficiency
Training time
Max
memory
Rel. L2 error
SGD
0.001
20,000
0.92 s/epoch
2h33min
328.11
MB
1.05e-02
Adam
0,001
20,000
0.95 s/epoch
2h39min
328.11
MB
2.31e-03
L-BFGS
200
26.1 s/epoch
1h27min
349.17
MB
nan
NGD
0.1
200
37.1 s/epoch 2h4min
328.11 MB
9.91e-04
A.7
NUMERICAL EXAMINATION MULTI-LAYER PINNS
While our convergence proof is based on the two-layer PINNs for simplicity, the extension to practical multi-layer PINNs are missing. The restriction of two-layer is primarily technical: it enables precise control of the NTK evolution and allows us to rigorously establish Jacobian stability (Lemma 4.6) and global convergence (Theorem 4.7). Extending these results to deeper networks is indeed possible but significantly more involved, as it requires layer-wise coupling analysis of the NTK (as in Allen-Zhu et al., 2019). Nevertheless, we report the NGD for different layer PINNs on 2D Poisson equation, to show the NGD can converge as depth increases. We train NGD with learning rate 0.1, and each trained with 200 epochs. Table / shows the convergence trends remain consistent with our theoretical predictions: NGD maintains a small relatively L error for different layers, while the overall convergence slows moderately depth increases. Especially the NGD's memory requirements keeps almost the same as layers increased (note that the inverse of JJJ" is independent of the parameters), and the computational burden only increase almost linearly with the total parameters.
Table 7: Training comparison of NGD for different layers on 2D Poisson equation.
Hidden layers
Total parameters Training efficiency
Training time Max memory
Rel. 1.2 error
512
3.67 s/epoch
12min 13s
14.75 MB
1.12e-04
33, 280
8.29 s/epoch
27min37s
14.75 MB
3.41e-04
82, 432
17.53 s/epoch
58min25s
20.17 MB
4.29e-04
19

---

Published as a conference paper at ICLR 2026
A.8 GENERALIZATION LOSS EXAMINATION
While our theoretical analysis focuses on the optimization of the training loss, the behavior of the generalization error after training is not covered by our current convergence guarantees and is beyond the scope of this work. Nevertheless, we provide an empirical study to examine how the number of collocation points affects overfitting for differentt optimizers the over-parameteriae regime. We consider the 2D Poisson equation and vary the total number of collocation points N = N + No used train the physics-informed loss (5). NGD trained with a learning rate of 0.1 for 200 epochs to ensure stable convergence of the training loss. SGD and Adam are trained for 20,000 epochs with learning rate le-3, and L-BFGS is run for one epoch with a maximum 50,000 iterations. To approximate the generalization error, we evaluate the physics-informed loss (5) on very fine grid with N2 = 100, 000 interior points and N₁ = 1, 000 boundary points. As shown in Table g increasing the number of collocation points consistently reduces the generalization error, bringing close to the training loss. For the 2D Poisson problem considered here, using approximately N 5, 000 samples appears sufficient mitigate overfitting while maintaining small generalization emoor across all examined optimizers.
Table 8: Generalization error comparison using different collocation points.
Optimizers
Training loss
N = 100
N = 500
NV = 1,000
N = 5, 000
N = 20, 000
SGD
2.13e-03
5,80c-02
1.37e-02
1.03e-02
2.59e-03
1.68c-03
Adam
9.7Ie-06
1.03e-03
4.59e-04
1.24e-04
3.41e-05
1.39e-05
L-BFGS
7.74c-06
9.16c-04
4.31e-05
3.93e-05
1.02e-05
8.65e-06
NGD
2.86e-06
2.51e-04
2.04e-05
1.22e-05
2.78e-06
2.91e-06
B PROOF OF SECTION 3
Before the proofs, we first define the event
Agr = (Ew ||w w, (0)l2 S R, [( x, 2 0} + I(w,,(0 2 0}}
(28)
for all E [n].
Note that the event happens if and only if (0) < ||x thus by the anti-concenntratio inequality Gaussian distribution, we have
2R
P(Arr) = P=..ioo.leee (121 < R) P=P-MNo.N (121 < R) S
(29)
2m
Let S, = (r E [m] T(A+)) = 0} and S.- = [m]]S..
Then, we need to recall that 0sp(w)
o (w). up Wrp + a' (ww xp))
1 oo"" "w, (w, xp))lumm 11:c0
20 20"
20"(w,,
0
Ow,
Vmn
0s41
w
(30)
and
@hy(w) Ow, / /mnz
ar
a' (w yubyy
(31)
B. PROOF LEMMA 3.3
Proof. In the following, aim to bound (0) - as (0) - [H(0) - I| Note that the entries H (0) H have three forms as follows.
(s:(w(() Ōs, (wo(1))
Os:(w(O)) as (w)(w))) Owo
Eww(o)
(32)
Owo
Ow,
Owo
Ds, (w) ahy (w)())
Ds. (s((w())) (wo
E Ewy(o)
(33)
Owr
Ow,
Our
Ow,
20

---

Published as a conference paper at ICLR 2026
and
@h,(w(()) (2(())()))
(h((w(0))
Thy (h,(2)(0))
>}
L
(34)
Sw,
Ow,
Ow,,
ow,
For the first form (30), to simplify the analysis, we let
(1) w. (0) (0)x; to
Oc+1
(0, 21,0) (0)] 20 (0))(0))=)) (o)
0
and
(ii)) (Z, Z,))
then
(sp(((())
)s,(w
Osp (sp(2((0)
as, (s)(2(())) ow,
Er
m Σ (i,( Ex,(
Ow,
Ow,
Ow,
Note that [x, S
(0)
thus
|
21+ 1|],,
< 11+ |1+1|/,,00)13
11 - d.
Here, for more details on the Orlicz norm, see the remarks after Lemma D.I. For the centered random variable, the property of by quasi-norm implies that
X (ij) - E[X, (c)llle S (1,)l + |||| (1)ll ddd.
Therefore, applying Lemma D.1 yields that with probability least - 8,
m
$ d
+ d2
log (11 )))²
(1 EX,(1))
log
m
m
8
m
which directly yields that (w)(w())) 0s, (w)(w))))
Osp(ww(0)
0s, (w)(w)))) ow,
d²
1
d²
1
(o)
log
+
log
ow,
Ow,
Owr
m
8
nim
à
(35)
Similarly, for the second form (31) and third form (32), we can deduce that (00 (h,(w(()) Ow,
< (s((w(()
(h,
d
Exx())
Or. Ow, and
dw,
1 m2
K Oh, (n,(()(0)) Ohz(wr(0)) Ow,
K
(Ohs(er(O)) ahy(w(()))
>
d²
Eur(o)
ow,
ow,
Dw,
n2m
Thus applying Lemma D.I yields that with probability at least
1 d,
T
(s((w(O) 0h Ow,
Xy
Os (s((w)())
Thy (h((w())))
de Vnnnz m
(
dc
Eu) and with probability at least 1 - 5,
log
log
Ow,
Ow.
Ow,
v2m
(36)
@l:(w)()) (w,(w(()) ) Ew)(o)
(h,(2r(0))
ah (h,(w((0))
d²
de
5
log
+
log
Ow,
Owr
Ow.
Ow,
12 R
72m
8
(37)
21

---

Published as a conference paper at ICLR 2026
Combining the above deduce that with probability at least - 6,
(0)
H
x0 1
[H H 11
22
log
711 + 112
+ m
dº
log
n1
+
182
6
6
d
111 + n22
log
5
Thus when
d4
log
n
t
S Ag i.e.
d'
112
m ?
log
A0
we have Amin Amin(H(()) (H())) 2 AAA
B.2 PROOF LEMMA 3.5
Proof. We first reformulate the term (sp(kk Ow.
in (28) as follows.
@sp(w) Ow, /mni σ (' xpp)
( 1 (0tt1) o" (w
w.o0p0 Wy0 - 2w,1)
+ (w xp)
xw)" wrr 11
suffices bound ((w) (O)lE which can in turn allows us bound each entry of H H(w) (0).
For E [m1] and j E [m1]. we have that
Hig (w)) (s,(w) Ow,
m
Ow,
+ +o' o'uw ( (0tt1) o" (w,"" xx:) 11 11:0; 2w,1
σ" (wf xii)
Wr0020
( (0t+1) - (w, (;;))
(w. x.))
Wr0Di0
o (w. xx)))
Wro0211 2w,1
After expanding the inner product term, we can find that although it has nine terms, it only consists of six classes. For simplicity, we use the following six symbols to represent the corresponding classes.
o 1,,00,0 σ 0,00
For instance, o represents
a (w." xx)
W00000
(w x;)
1 o (w) x;)
Wro0io Wy000,1 2w,
Wro0n1 2w,1
0a+1
0
In fact, when bounding the corresponding terms for Hiy(w) (0)( the first four classes can be grouped into one category. They are of the form fi ((w)) (w)), where for each i (1 S 4), fi(w) is Lipschitz continuous with respect
||| and [f:(w)] X ||wl|| (Note that ()) = (a ).))²) On the other hand, when [|w w2ll2 S R S 1, we can deduce that
(w>)) S R( w1 112 + 1). Thus, for the terms in Hyg (w) He, (0) that belong to the first four classes, we can deduce that they are less than CR(] (0) + 1), where is a universal constant.
Ifi (ww))ff(( (etee//((wwww) (w:)fs(w))))) (f. - fi
For the classes o σ and o o they are both involving o" that is not Lipschitz continuous. To make it precise, we write the class o a explicitly as follows.
σ" (. )n) (w, x 1
Wy0Cin
Wroa011
22,,1
Dy-
22

---

Published as a conference paper at ICLR 2026
Note that when ||,,- w. (,) < R, we have that
lo (w xbj) (w, (0)) 2 = [ [1 æj 2 03 (w,(()) x 2 of {(Ar)
where the event Air has been defined in (36).
Thus, we can deduce that for the terms in (w)(w) (0) that belong to the classes of a and they are less than
C' (( (0)) + 1) + R( w. (0) + +1)]
where C' is a universal constant.
Similarly, for the last class o² that are of the form
o (w az)o " (w w æy,
we can deduce that
lo "" xu)o "al)""""("ii (w, xz) 1, 2, (w,(0) 22)0 (w,, 10)))
S S(Ar How,, (0)11 + R(| w. (0) 1).
Combining the upper bounds for the terms in the six classes, we have that
1
1
m
[Hl((w) (1,
R
w. (0)
(1
+ H(A,) ((A0)) +w,,
R
m1
m
m
1
1
(0)1
Σ
M
[ 144
m
m
(38)
where the last inequality follows from that (0)]1 2 n0, (0)13 + due to Young's inequality for products.
Now, we focus on the term
20,
Since
8 )2) <
P (0)= 2 2 log
and then
2(( +)) ))
[|w, (0)113 2 2(d + 2) log
S
This implies that
221(d ( )
3r E [m], lw (0)11) 2 2(d + 2) log
s 8.
(39)
8
Let M = 2(d + log 2m(d+2) then
1(4, H| w. (0)1
I{A (10,0)140| < M) + Σ HA } w. (0) > M}
r=1
2 m
Σ [[Art Iw (0)) w. (0) M}.
m
Applying Bernstein's 's inequality for the first term yields that with probability at least e em
1 < 4R.
23

---

Published as a conference paper at ICLR 2026
Moreover, from (39), we have that with probability least 1 - 6, the second term I (0)] M = holds for all r E [m].
mnt mme--m, Thus with probability at least -
we have that for any [ni] and j [n1].
1
[1s(w) (1)(
RRM + RM R
11
1
X
M
For i E [ni].j E [n1 + 2, m2] and i E [n1 + 1, n2l, E [122], from the form of 0hy(w)) in0.
i.e.,
ah5(w)
ar
o' (w,, y;)y
Ow,
Vi2mm
we can obtain similar results for the terms
Osi Oue
0h, aur
and
Olss 8he
8h:
Oro
With all results above, we have that with probability at least 1 - 8 - -mR
mn
(H(w) HO)l $ M
B.3 PROOF OF LEMMA B.1
Indeed, the stringent requirement of the learning rate in Du al. decomposition method for the residual. Specifically, in Gao et al. residual in the (k 1)-th iteration is same as the one Du al, (2018), i.e.,
(0018)
stems from an inadequate
(202)) the decomposition for the
sk + 1) (bk+1)
(((R))
s(k)
+ ( s(k + +3)
(m(
(k)
(40)
h(k))
h(k))
which leads to the requirements that = ((o)) and m = (oly(n1, n2, 1/6). Thus, requires a new approach achieve the improvements for y and m. In fact, we can derive the following recursion formula.
Lemma B.1. For all is E N, we have
s(k + +1)
((k)
h(k + +1)) = = -I (H:))
h(k))
(l,),
(41)
where
(k) = (A), ... , tm (k:)) E tm2
and for [n1].
(k)) sp(k 1) -p((k) - ( [ksp(k)
wi( + 1) - -20(k))
( Ow
(42)
for [mz],
/Oh5(k) Ow
(k) = h15k + 1) -1;(k)
w + 1
w wo(k))
(43)
In the recursion formula (39), (k) serves as residual term. From the proof, we can see that [| (k:]l2 (1/vm)) and thus, m becomes large enough, only the term I I- (k) is significant. This observation is the reason for the requirement of 17.
Proof. First, have
- (sp(())
wik + 1) -w(k)) )) +
(Osp(k)) [sp((k)
8p + - sp(k) = sppkk (k: + 11 - -kp(k)-
Ir w + 1) zw(k))
Oro
== (k)) + (k)...
(44)
24

---

Published as a conference paper at ICLR 2026
For the second term (k), from the updating rule of gradient descent, have that @sp(k)
(k)
w 1) 20({))
dw
(kp(k)
θL(k)
Ow
Oue
(k) OL(k)
Ow. Sw,
(sp(k)
sk( (s)(k)
(h1(k)
(45)
+
(k)
Ow.
==1 n2
n1
(k)(p) (k)(k)
@sp(k) [h2(k) )
!
se(1:) +
Sne, Ene, (k)(k)
Oen
Sen,
(s(k) (2(k))
(k))
where (k)p denotes the p-row of H (k).
Similarly, for h(k), we have
[hh((k)
w(k + 1) - -2(kk)) )1
(hy(k)
(1 + 1) (A)) = 1gy(k + 11)
w + 1) 22)(k)
Ew
Ow
1 (k) + 12. (k)
(46)
and
(s()) h((k)
+ (
(47)
Combining (42), (43), (44) and (45) yields that
s(k + 1) (s()))
hik (((+1)) (1:))
[1) + I2())
(s(k))
(k) mk(
(h(k)
A simple transformation directly leads to
s(k + = (I (k))
((k)
h( h(k+1))
h(k)
+ (1s),
which is exactly (39), the new recursion formula we need to prove.
B.4 PROOF OF THEOREM 3.7
Similar to Du al. (0018 and Gao al. 2023). we prove Theorem 3.7 by induction. Our induction hypothesis is the following convergence rate of the empirical loss and upper bounds for the weights.
Condition 2. At the t-th iteration, we have that for each E [m]. [0), S B and
1700
I(t)= <
L(O),
(48)
2
where B
2(d 2) 2m(d+2)
+ 1 and L(k) is an abbreviation of L(w(k)).
From the update formula of gradient descent, we can directly derive the following corollary, which indicates that under overrparametteiizaton,, the weights are closed to their initializations.
25

---

Published as a conference paper at ICLR 2026
Corollary B.2. If Condition 2 holds for 0, then we have for every [m]]
||,, (I: + 1) u. (0)]l2 CB2 R,
(49)
Vmdo
where C' is a universal constant.
Proof Sketch: Assume that Condition 2 holds for 1 O,... ,, it suffices to demonstrate that Condition 2 also holds for t = k +
From the recursion formula (40), have that
s(k + 1) h(k + 1)
[ (k)))
(((k) (((k)
I (1:)
(50)
(s())
2
(1(kk) (((k) (B)]l2
< I-- (
+ (ITT
2 - n (b)1l2
(((k))
where the inequality follows from the Cauchy's inequality.
Combining Corollary B.2 with Lemma 3.5, we can deduce that when is large enough, have
(k) H())|22 < oo/4. Thus, Amin (H(k) 2 2n/2 and I (k) is positive definite
when 77 = 0(1/ Hz))2)
l|2). On the other hand, with Corollary B.2, we can derive that (k) =
OknV/(k)/ /m). Plugging these results into (48), we have
skk + 1) h(k +
1120
2
1
((k)) (kkk)
+
O
0
(51)
m
mm
n
8 (A:)
2
he
where the last incquality holds when m is large enough.
Now come to prove Theorem 3.7.
Proof. Corollary B.2 implies that when m is large enough, we have [l,, then (k, + 1 1|2 < B. Thus, in induction, only need prove that (46) also holds for = k: 1, which relies on the recursion formula (49).
(k+++) 1)- - vr,(0))2 - 1 and
Recall that the recursion formula is
s (k + (1k + +1) (I- (k))
(s())
(k).
h(k))
From Corollary B.2 and Lemma 3.5, taking CM < Amin (H())) Amin (H()) Aa 2 z and
Ag 7
in (8) and R¹ < R in (47) yields that
An
No S H +
A 2
H S H (0) 2 +
H
4
Therefore, if we take 7) S 3 2 then I (k) is positive definite and - (k)
ndo
Combining these facts with the recursion formula, we have that
s(k: 1)
h(k + 1)
(sk))
2 + [1))(1)]13 +2 ( (I- (HH)
(s(k)
II nH
(h(k)) _I((k)
(52)
(((k)
2
< (1-
n7d00>
)
2
(k))
22-11-11 mdo
(s(k)
[1) +2
[1)
2
(1k)
26

---

Published as a conference paper at ICLR 2026
Thus, it remains only to bound []]()]]22
For I (k), recall that I1 (k) (k), ,1' 11 (())),)))))),,,, ,, , E mutr and for /Osp(k) Sw wokk + 1) - 20(k)
p E [mi].
(k) sy(k + 1) - -kp(k)
for j [n22].
(b) h1 + 1)) { [Oh((k)
(1g(k))
:w(:: + 1) -wo(t)
Our
Recall that
5p(k)
1
Σ a (w,(k)? (.))T(())))))) 0,- (w) ap))1 w/l (1)113 - f(zx))
(k))T)))llll
Vm
and
(sp(k)
ar /nim
( (x))so(k)))) +o (ur(k)? xp)) ( 1 (0t+2) - -o"" (w,,k))) xp)) /. [p))11rr11 1Bar
σ " (m0,
Ow.
-20" " (w,, xxp
w (w))((t))
Define Xpn (k) == o' (w,( (1)k))=x)ww)ap))) (k))) and (k) == o" (w, (x,)) (k))1. ie., X (1)) and
(k) (k) are related to the operators de Then define
Ou and Au respectively.
(xpr Ow.
(k) = (k +1) -x (b:)
.,, 1) - 1) (1,(k)
and
(k) Ow.
(k)(k) = Xpr (k: (x)
w. (k. + 1) wk)(k)
At this time, we have
(R))
mmm Σ ar
[1 (k)
X (R:)
The purpose of defining X,r (k) and (k) in this way is to enable us to handle the terms related to the operators dt and Au separately.
By
We first recall some definitions. For p E [n1].
Ap. (Bw: ||w w. (0)l22 S R, l( x, 2 20 + 1(w 2 20}
and S, = (r E [m] I [Ap.r = 0}}, S. = = [mi](Spp
In the following, we are going to show that 1x] (R:) = OOl w. 11) 1 - w. (k))))) for every r E [m] and = W[ (k- + 1) w. (k)) for E Sp. 18 (k:) = O [|,,,, (k + 1) w, w,,(k)) (k) for
E Sp. Thus, we can prove that [k)]2 O ( (VL(K)
Then combining with (69) leads to the
vm
conclusion.
For X.r (k),
(k), from its definition, we have that
Xbr (k) = o' (w, (k + 1)T woo (k + 1) o' (k)) xx wrokk)
(w + 1) - w, (k)) xp)0 (k)² (0) - (wro + 1) - wro(k))o (w, (k) xp))))) (a' (w, 1)T ax,)) a' (x)" (w)) 1) (w, + 11) -m,,(k), x,)0 (w, (R:)² wp)))
27

---

Published as a conference paper at ICLR 2026
From the mean value theorem, can deduce that there exists ((k) E R such that o (wr(k + 1)a xp)) - o' (w, (k)) xp) o" (s(K)) (ks +) (r)(k). xp)
and
" (k)) xp)p 1 (S(k)) - w. (A:) 10 ({(+))) - o
/ 1. (k- 1) w. (N:)l22
Then, for Xpr (k), we can rewrite it as follows.
Ror = o (5(k:) (kn, + 1) - w. (b)), up) + 1) (w, (k: + 1) w. (k), ap)d" (wor ay)wro(k
= =I(o" (C((6)) (((())))"""(w,(
(w, (k + 1) - w, (k)) xp w (k + 1) + [(o + w. (k), 0 (k)) xp)) (1rr + 1) w)))
This implies that
X (1) Bl (k+ + 1] - wr
For Xpr (k). we write it as follows explicitly.
Xpr = o (w, (k+ + 1) x,)) 1 (k + 1)12 - o. (w, (A:) T xp)) 20,1 (1)112
(w, + 1) - w. w,) xp) (w, (k)² xp) w (1)
(53)
2/w,1 + 1) - w (k), w. (k)) (w. (k)))
Note that for the term o" (w.((k)''u,)))) ,,,)| (1: we can rewrite it as follows.
" (k)² 2p) (1))
= (w, (k) T xp |0, (k) wr (k- + 1) + w. (k- + 1)
=o" (w, xx ) (k:) w.1 (k 1)1l + ||0,1 (k + 1) 2(w,1 (k + 1) wx)(k),, wri 1)91,
(54)
where the first term o² (w, (k))#>>))lpp))))))) up (k) - w,, (k- + 1) = O lr, +1) (k)1)).
Plugging (52) into (51) yields that Xfx (k) = (o" (w, (k + 1)" mp)
(w, (w,(k)) x. 1 (k- 1)112
- (w, + 1) - w. (k). ),)0 o (k, 2,,)1 || w 11:)11
+ 2(wri + 1) - 1 (k), w (k: + 1) - w., (k)) (w, (k) T xx)) + + (BBlu,, + - w. (1) + +1) xp) ())) σ' (w)" xp) (w + 1) 1 (k),,, xp)o (w,,)" (w,(k) xx)p2 (k)pTT))1)) w,) (A: + 1)112
+ (w, 11- (k), ap)o (w, (k) x>)xxxppp100))> xz,) (1 + +)111 w (1)
+ +OBBlu, 1) (k))) (k)
Io" (w, (ks + (1) xp) - o (k) T xp)) - (20, + 1) wr)(),, (w,,(t)Tzp)) || (k + 11112
+ +O[Bl w. (k + 1) w. (k)
(55)
Thus, we only need to consider the term
o" (k + 1)T=,, ax,)) σ² (w,, (k:) xxp)) (w, (k: + 1) (k)(k), (w,,
(w, (w) xr))
(k)))ar)))
For Sp. since || w. (k + 1) - w. (0) 2|2 < R, |wo, (k) w.)ll S R. we have that I kw, (k 1)th x 2 0} ((u, (k) x, 2 0). which yields that
o" (w, + (
= [( + ow"
- oo (w, x2,) (w, (k- + 1) (,),),, 2 " (k) x,))
2 0} - (w, (k) (p) (k) x 2 0]]
- + [u, x, 2 0}
(56)
[[w + 1) x, 2 0] (w,(k)
x 2H1
(w.
+ (kw
(z), (k)" x, 2 0)
O.
28

---

Published as a conference paper at ICLR 2026
For E Sp' the Lipschitz continuity of o" implies that
(k+))))))xx))>->>> (w,(k)² ( -k. tp)) 0 (w, (k)) - (k,1)-
(57)
Combining (53), (54) and (55), we can deduce that for r E Sp:
(k) & Bllee, 1) w. (1:)11
and for r E Sp.
1X 12)% (K)] & Bl| (lc 1) (k) B² 1|,,, 1) (x)kll2.
With the estimations for Xpr (k) and Xpr (k), we have
m
[P(()
(1 pr (A:] 1xp (k)))
(58)
Bl (k + w. (A)122
Σ
B² B2,, ww
+ w,(()ll2.
mm
Vim
rES,
For j E [12]. we consider 7 [k=+)()), which can be written as follows.
(k)+((k)= = hyglk + 1) - (k; (o(k +) - wkk), (h)
dw
m
0r
oww (k + 1)) (k) 41) ) (ww, (k + +1 -., ysdo (wr(()?
12271
From the mean value theorem, we have that there exists ((k) E IR such that
aww, (k + 1)T yj) (w,,(k)) y,) = a' (un + 1) -0, -.,(())) (R:). ms)
and
10' (C(K)) o' (w, y/)1 2B)((() (k)(k))
s 22v2BBlw,, (k + +1) (, (A:)l2.
Thus.
loww + 1) yj) - (k) y;) - (w, (I: + 1) - u. (k),, 3/)00' (w= y3)l wk),), 81) (k), (k)² yソ
(w, (k: + 1) w. (k),, uo3o (k) v;)
= 1(0' (C(K)) - (w, (k) (w)) (k+ + 1) w, (R:). y
Bllw, (k + 1) w. (k))l2
Therefore, for j E [n2]
1 (k)+
Iw + - 2
(59)
From the updating rule of gradient descent, we can deduce that for every r E [m].
al(k) Ow,,
BB²
|w (1: + 1) w,,(k]l2 =
L(k).
(60)
Vm
29

---

Published as a conference paper at ICLR 2026
Plugging (58) into (57) and (56), can deduce that
B
m
B²
[P((R] x
+ 1)
(k))1
||w, + 1)
(1)
rest
B
m
B²
B
BB
L(k)
L(k))
7mm
m
Vm
res,
(61)
L())
BB
L(k)
1
1rr E $
Vnim
Vn1
m
BB BZ V /OO L(k)
nBA
[(k))
Σ Arr E 5
Vnm
v n1
and
B
m
(
(1:)
wr (k 1) w. (k)]l2
/n2m
17 Ba
L(k)
(62)
v m2m
n B
(0) v L(k)
/n2m
Note that
2R
P(Ap.)
2 Sp fr [m] [{Ap-) = 0}.
Thus, from Bernstein's 's inequality, we have that with probability least 1 e-R
1 s m Σ 1(pr) 4R.
172 Σ Hre S
Then the inequality holds for all p E [na] with probability at least 1 - Nie m mIl Plugging this into (59), we can conclude that for every p E [ni]
BB LLO) (OO]LR (T(K)
nBt LT())
+
R.
(63)
11
Combining (60) and (61), we have that
T1
1 (k)
R(KP + Σ
(k)
B
L / L(K)
+ nB* L(k)R.
Plugging this into (50) yields that
(k +
1)
|
hkk
1)
1 1A
(s(A))
11
(A:)
2
+ (1)]11) + (1- ndO
(xxte)
(((k)
|
(b)1l
(k(k))
2
(((k)
2
nAo
1T L[())
2
11² /((B)
s(k) 2
+ C
+ +nBBR B
+ 20
+ nBBRR
2
m
m
h
nu
((k)
I
(1:))
30

---

Published as a conference paper at ICLR 2026
where C' is a universal constant and the last inequality requires that
12 LT())
S n
nB R S ndo-
mm
Recall that we also require CM²R
for R in (8) and
CB L()) < R
R
vmdo
for R in (47) to make sure (k)) H (0)22 S 4.
Finally, with R O Ao (x)))) and the upper bound of I((0), m needs satisfies that
m M+BBL(0))
d
log
6
md
log
n1
12
A
1
C PROOF OF SECTION 4
C.I PROOF OF LEMMA 4.4
Proof. Recall that
(w) dw
0sn
@hi(w)
Ōhag (w) dw
=(w)=DD) D
Ow
Iw
and
H (o H (w)).
We denote so w) = o' (w x: x)w0 - o" (wTTx)]lwill [ where w = (wo, w w0 E R, wi E Rd, then
@sp(w) 1 ar @p(ap:) ow, Vnn vm Ow. Similarly, we denote (yy w) o (wT y). then
(hy(w)
a, Oey w.) Ow Vnz vm Ow.
With the notations, we can deduce that
aep(ppi w)
ap(ay w))
E (
1 I m,, S $ n1,
dw Ow Osp(pp w) Ow Owo Orp(yp; w)
H
N(O.I)
15pSm, +1 1 mm +
aupyy; w) Ow
E (
Th + S m + n2, m + +1 < $ m + 122
12
Ew
where is the (ph entry of 0 p.j
The proof of this lemma requires tools from functional analysis. Let 7. be Hilbert space of integrable (d + -idmennniol vector fields on IR4+2 i.c., f H. if (l < co. The inner product for any two elements f. g in H. Ew-N(0, ((f()))) positive definite is equivalent to show that
g( (ww)). Thus, proving o strictly
@ Ewo
ap(an w)
Jup(yng:
E H
Owo
Ow
Ow
linearly independent. Suppose that there are an anj 31,
Brg E R such that
Üp((11 + + On,
dp(an,;w) + B1 dw
00 + .... Bng
(pzi w)
an
= 0 in H.
dw
dw
dw
31

---

Published as a conference paper at ICLR 2026
This implies that
0p([): 01
Ocp(ann w)
@ub(yii))
Irb (r w)
an
B1
+
+ Bug
(64)
Ou
Ou
Eue
holds for all w IR++1 (.) is smooth.
First, we compute the derivatives of şp and da with respect to w. For the k-th derivative of (by; w) with respect to w, we have
0 bur ((0)(w' y)y
Dw
where & denotes the tensor product.
For (0 w), let (a;; w) = a (wo awwn and yo (x;; w) o" w)m for Sis d. Then
( w)) = (0 w)
(x;; w.).
For the k-th derivative of (o(æ; w) with respect w. analogous to the Leibniz rule for the k-th derivative of the product of two scalar functions, we have
0 po w) Ow
((++1) T w0 (k))
x (k-
(00
E0
(65)
= [k++)) (w= x woe (k) (k)
(w² x
x ((A)
where (1,0, .0) E Rd+1 and in the second equality, e denotes that is placed at the i-th position.
Similarly, for (p;( un) where sis d, taking i example, we have
0 (xx; w)
((++)) x w 2kw o (k+1)
xx) 00 0(k-1)
Owk
(66)
k:k 1)0 (w x)
(k=22)
iSKjjk
where e, E Rdd+1 a vector with the (i + 1)-th component equal to and all other components equal (4) indicates that is placed at the i-th position.
to 0, and ei
By combining the derivatives of (ar( un) ... paaa;; w) from equations (65) and (66), we can compute the k-th derivative of woo (w- x) Σ wfo" (w"")) as follows:
I @ky((:; Owk
Ok po(l;; w)
d 2h v)
Ow
Ook
= -ho (k+1) ( (2==) xx) (1k) + o (K) (w0 x)
x x 0(k
(67)
011) - 2kw: (1+1) (wo x))
e (2)
(kk-1)
+k(k (8) (x)
Σ
el0 8 _el 2 (k)
(1)
1Siejsk
Note that when any two points in (x1,.... Fni ... Yn, } nar-pllll,el, the tensors
8(m .... m
+22) yi
((mltn2) y
S(n +22)
32

---

Published as a conference paper at ICLR 2026
linearly independent (see Lemma G.6 in Du al. the (k 1) derivative of both sides of equation (64) with respect to w, yielding
(2019)) This observation motivates to take
0h (()1;)) +... + ani 0 ( w) Ow
Ih ( w)
2 s ( w)
+ B1
0. (68)
Owk
Ow
Since this equation holds for all w E gd+1, we specifically consider w = (wo. Oa), where wo is to be determined. Under this condition, equation (68) becomes
W0
Cpp
o (ktt)) (wazh))
+ p=1
(0))) (woa)))))) + Bg (k) (w0v5)))))) (01k) = O,
Σ ==1
p=1
(69)
where the tensor Zp defined as
2p
X ((R - lok 1)
Σ Σ
ei
_ @ 2
(70)
t=1
Isijjsk
((m)(x))
By assumption, for any positive integer n 2 O we have
lim
On + 0. We first consider
210000
[(x)
the case where all input components satisfy we 2 condition, equation (69) simplifies to:
yh
0. Under this
woo (k+ 1) (wna)
Σ
ap (0)))
+ o (8) (wna)
Σ ap
+ to (woaa)
Σ
3,9,
(K)
= 0. (71)
Dividing both sides of equation (71) by ((waa) yields
(wot)) (wwa)
(k) (wna)
(k) (w)a) o(woa)
wo
a
@(k)
o(waa)
Σ xp
+
3
(k)
= 0. (72)
(A) (wae) converges to non-zero ((woa)
Now taking the limit as we tends to positive infinity, observe that
(m)a) diverges to infinity. This asymptotic behavior leads to the following constant, while wo d(woa)
conclusions:
Σ
pp (k)
O,
2p
3
0(k) 0.
(73)
1=1
By the linear independence of the tensor products (established earlier), we can deduce that op for all p 1,.., Ri, which subsequently implies
122 Σ B,y, ((k) and thus Bj 0 for all
j=1
j 712
When our previous assumption does not hold that is, when x9. ... 2 v. ... vn are not neces- sarily all equal-we we proceed as follows:
Case 1: All elements are strictly positive.
Let b = minf2].- ... 20, vl. ... v }. Dividing both sides of (68) by d(wob), we observe that for any x > b,
(wnoa)
lim ++00
d(wyb) = 0.
Thus, the problem reduces to the previously considered case where all inputs are equal (to b). Due to linear independence, the coefficients Op, and B3 corresponding to the minimal b must vanish. Repeating this process iteratively, we conclude that all ap and 33 must be zero.
Case 2: Some elements zero.
33

---

Published as a conference paper at ICLR 2026
Since the zp interior points (and thus non-zero), any zero- zero- valued inputs must correspond to boundary or initial conditions. Let assume without loss of generality that vi. . . are all zero.
1. If or (0) 0, then following previous method,
conclude that the coefficients corresponding to non-zero inputs vanish. Returning to equation (68), we have
Σ 3,0 (K) (ww リ,ソ,, ((k) O.
(74)
By the independence of yl, Mog. we can deduce that 3 (w0" yj) holds for all j E [m2]]-
From the assumption, we can select w such that a (4) (w" yj) # O, and consequently, By = 0 holds for all j [mz]]-
If a (k) (0) O. let b be the smallest strictly positive value among = (71) by ((wbb/)). Since all other positive terms decay to zero as + +oc, we obtain:
In, Vi w/e Divide
fig
lim
3 (8) (0)g, (kk) /o(wab/2) O.
(75)
100
(1)) This implies that Σ B,y, = O. By linear independence, all [j = 0.
Remark C.1. The key point in the proof lies in the fact that the order of the PDE in the interior is higher than that of the initial and boundary conditions, allowing for natural extension broader classes of PDEs. For general PDEs, we may focus solely the interior and boundary, assuming the interior is of second order and the boundary of first order. Suppose the second-order interior term
is taken at Zo. i.e., it has the form De?
&y and the first-order boundary term is also taken at zo. Since we assume that all cn-commonents
translate the coordinates, without loss of generality, positive.
(2) x)) yields that For the interior, taking the k-th derivative of w (2)
k
w ()+ (w= x) m())))))))) + 2kwng (k+1) ( x)
e s0
0k--1
k(k - + (N
i
x (k)
ISissk
For the boundary, taking the k-th derivative of (w To x) yields that
k
I (k+ ((+1) ()) x)wa + (w
@ a (0)
As before, we set w = (wo, 0). Then, the equation (69) becomes
n1
w Σ Op
[ (k(+1)
(woz)))
+ we
Opp
[ (k) (w0orr))))))
+ Σ ap (o)t
n2
12
Bg [ (k) (1))(woo)
((k)
Σ By (k) (woB)) ( 12]
+ wD
=0.
j=1
where zp, =}: 23 are tensors of similar form Zp in (70), whose explicit definitions are omitted for simplicity. Dividing both sides by myo(wg=() reduces it to the form considered earlier. We can therefore conclude that the Gram matrix is strictly positive definite. Indeed, since the orders of the interior and boundary terms in the partial differential equation differ, we relax the conditions in Lemma 4.4 to simply requiring that no two samples in (xI ppm are parallel and no two samples in {y, } are parallel. In brief, we can set = na - and = n2 1 in equation (68), and then use the method described above.
nn
34

---

Published as a conference paper at ICLR 2026
Remark C.2. For the activation functions sin(x) and cos(z). in equation (69), may assume that o (1++1)()) (x) sin(x), (z)) cos(r). Dividing both sides of equation (69) by wy and letting
Wp + +00, we can obtain
lim
[(k) = 0. sin (e (w)z)) p
p=1
We express the general form of the components of the tensor above as
Cppp sin(wpz,)) = Σ a (in(b))))),
where the b, > O are distinct and Cp denotes the components of the tensor app (k), For simplicity, we
denote wy as æ. To prove that apTp @(k)
0, we need to show that any component of this tensor
is zero, i.e., its general form satisfies
CpCp = 0. This is equivalent to proving Σ 0g = O.
Let f(x) Σ a, sin(b,x), note that dividing both sides of equation (69) by wp yields that f(x)
(1/zz). Thus, consider the average energy of (z) over the interval [T, L). i...
L 1 T+L
((x)dr.
Expanding this, we obtain
2
a sin(b:z) dx
L
T
1
TTL
de sin 2 (b:z) + Σ ajay sin(b,a) sin(bbat) dx
L
TT
1=1
i
18
1
sin (2b (T I.)) (2
0
L
4b;
sin (b: b;) I.)) sin ((b: b;)7)
a;aj
2(b: bU)
1
(b; + by) + L)) sin ((b; + b,))))
a,j
2(b; by)
Taking the limits L + too and T + too in the above equation, the right-hand side tends to 1 Σ a?
1=1
Regarding the left-hand side, recall that f(x) (1/xx), thus for any 0, there exists To such that for æ > To. [f(x) < e. Therefore, for > To and any L, we have
T+L
T. T (xx)do=se
By the arbitrariness of e, we can deduce that
1 T+L
lim
(x))dx= O.
T.
JT
n
"
Hence, we can deduce that 1 Σ a² = 0. which implies that Σ a, = 0, i.c. p=1
Σ
apC = 0. Finally,
we obtain
(kk 0.
a
p
35

---

Published as a conference paper at ICLR 2026
Applying the same approach as before, conclude that dp = 0 for all p E [ni 3j 0 for all
j E [n2]-
C.2 PROOF OF LEMMA 4.6
Proof.
Recall that
Ösp(w)
(
o (ww xp)wroop + +o (wf xp)
o"" (w x 12,1 12
Ow,
Oz+1
0 -2 xp)) (u0)2))
and
(1,(w)
ar ow, Vrgm
o' (wf yo y,)
(1) When of-) is the ReLU activation function. From the form (sp(w)
Dw. we can deduce that
Ösp(w)
Osp(0)
Ow,
Our,
[R[[w, (0) + 1) [ xp 2 0} 1(w,(0)²x;; x, 2 0}((1 w. (0)1 + 1)] (76)
Vnmm
[[(]w, + 1) + ((r [1 (|]l,(())) (0)113 + 1)]
where the second inequality follows from the fact ||w w,) R < 1 and the definition of Apr in (28).
Similarly, we have that
0hy(w) (h) Sw. 2 V22m Combining the above equations, deduce that (w) - (0)
R( + 1).
(77)
Ow,
[J(w) J(OI]
Bi 41
Σ
[|:(ww 10)1
(w)
(0) Oun,
122
(h,
(h)
Sw
Ow,
Ow, 2
na
(R(I (0)l22 + + 1 (4 pr { (1, (0))3 10013++0 + 1))² 2 +
>
(ill2w [1)]2, R2)²
mmm
nmm
j=1
R2
(0))1 1)
[1 21 (1
1
R2
w| (lw (0)l3 1)
1=1
(Apr)
< + (1 > M + 1)
mm
R2
M
w. (0)11 11)
[App)]
w. (0) w. (0)]3 M),
m
Ti
m
p=
36

---

Published as a conference paper at ICLR 2026
where M = 2(d + 2) log(2 2)/6). Note that from (39), we have
P 3r [m], (0)l 2(d + log (2m(d+ + 2) 2) 8.
8
On the other hand, applying Bernstein's inequality yields that with probability at least 1 N1e n mm
m
Σ 1mp <
4R
holds for all p [n.]. Therefore, we have that
( - (O)] S MR + R² - M M²R
holds with probability at least - 6 n1e mff
(2) Note that when a satisfies Assumption 4.3, o o and o " are all Lipschitz continuous and bounded. Thus, we can obtain that
(w) (sp(O) 2 Ow, Ow,
R(I w, (0) 2 w (0,,2 & Vnim
R(w (0]1, (78)
where the second inequality is from Young's inequality.
Similarly, we have
(w) (hz(0) Sw,
2 R(]),,(0)]lz + 1).
(79)
ow,
12
/n2m
Combining the above equations yields that [J(w) J(0)
@sp(w) (sp(() Ow,
Öhy(w)
0hz(0)
+ Σ
Ow,
Ow,
Ow,
1
1
(Rllw, + R2) +
[H]]w,(())1 + R2)
n2m
R2
(ww,(()) 1)
d²
d²
1
2A
R2
log
log
Vm
m
where the last inequality follows from the fact that
|||| [||,,, (0)
S d2 and Lemma D.I.
C.3 PROOF OF THEOREM 4.7
For the sake of completeness in the proof, we restate Condition and Corollary 4.11 from the main text, and label them as Condition and Corollary C.3, respectively.
Condition 3. At the t-th iteration, we have
[w (t)l - B and
w. (t) w. (0)l2 < CB2 L(O) RRR
Vmdo
for all [m], where C' is a universal constant and B - 2(d 2) log
(2m(d+2)
+1. 1.
37

---

Published as a conference paper at ICLR 2026
Corollary C.3. If Condition 3 holds for t 0,... and R R and R" 2 vI Au, then L(t) S (1 (1-11) L(0),
holds for = 0,...,.. where R is the constant in Lemma 4.5 and R CMM.RRRR R in (16) when a is the ReLU activation function, R = CdR. in (18) when a satisfies Assumption 4.3.
Thanks to Corollary C.3, sufficient to prove that Condition also holds for = k: + 1. For readability, we defer the proof of Corollary C.3 to the end of this section. In the following, we are going to show that the Condition also holds for Corollary C.3 leads to Theorem 4.7.
1, thus combining Condition 3 and
(s()) First, let utt)
Sketch Proof of Theorem 4.7.
then from the updating formula of NGD
(h(t)
(11), we have
u(t + 1) - u())
= u (w(t) -nJ (t)² H (t) ((t)))=(((tt)) u(z((t))
(il(w())) JJJ [t) (t) (t))(((((((()) ds
Ow
(u((w(()
(80)
11 J(1) (1)) ((( ds
Ow
+ @r((w(t)) @u(w(())
(t) H (t) (t)))((())))))))))))))))()))(( ds
Ow
Ow
== (t) + I2(t),
where the second equality is from the fundamental theorem of calculus and w(s) swi( + 1) (1 - s)w((() = wr((= (1)" (())-=u((1.
In the proof, assume that Condition 3 holds for 0, ... kk. Then from Corollary C.3, to prove Theorem 4.7, it suffices to demonstrate that this condition also holds for tort- = k+1. Here, we primarily explain the process from Condition 3 Corollary C.3, while other content is placed in the following full proof of Theorem 4.7.
Note that Pu(w(()) (t), thus (t) (t). Plugging this into (80) yields that u(l 1) (1 (() + (2)..
(81)
From the above equation, we can see the difference between NGD and GD. Recall that the iteration formula for GD
ull = -1 (1-1)(()))(() + (1)
Precisely because of this, the convergence rate of GD is inevitably influenced by Ào whereas that of NGD is not.
From the stability of the Jacobian matrix, we can deduce that [[2(()]l2 ((nllu(t)la///m Plugging this into (81) yields that
| 1)112 < [1 11)-n())11111
[+) 2(1 - ()) [2(1)
(1 (1- 2
1 m 2(1
1)0
(m
-
[[](
(82)
< (1 - (1)))))
where the last inequality holds is large enough.
Full Proof of Theorem 4.7. Recall that we let R" CM /R in (16) when o is the ReLU activation function and let R" CdR in (18) when satisfies Assumption 4.3.
38

---

Published as a conference paper at ICLR 2026
First, we set R' R and R" 6m since R² s VI /VAo
Then from Lemma 4.5 we have
[1](1) J(0)l22 <
A thus
v3A0 v3Ao
v3
Omin (J(1)) 2 min (J(0) - [1](1) - J(0)l22 2
2
3
and then Amin (L)) 2 Aa From the updating rule of NGD, we have
3 for = 0,... Ai, where ominn)) denotes the least singular value.
s(t) ht(t)
w. (t + 1) (t) 17 [J(t) (H()))
where
[J(t) (ts(t) Dsn, (t) (h)(t)
Uhn (t)
Sw.
Ow. Ow,
Owr
Therefore, for 1 0,... and any E [m], we have
(t, (t+ 1) (t)
7| [J(t) 11 H (t) L(t)
31
Ào
3//
[J(1)",, l1 L(t)
$ Ao
37
0sp(t)
n2
(ty
V((I)
(83)
No
Oae
Ow.
j=1
B+ +1 1
V((t)
Ao
BB²
mmAo L(t))
nB2
(1 - -22 VL(O),
vmAo
where the last inequality is due to Corollary C.3. Summing from 0 to k yields that
||,, (k: + 1) - (0)l2
(t, + 1) uv
C BB²
(1 1/222 /)(1)2/((0) L(0)
mmo
CB² L(O)) /TO)
Vmdo
where C' is a universal constant.
Now, when R' < 1, we can deduce that w, + 1))l2 < B, implying that Condition also holds for ts + 1. Thus, it remains only to derive the requirement for m..
CB²VL(O and R" -I-TV/
Recall that we need m satisfy that R'
(1) When a is the R activation function, in Corollary C.3, R = CMVR S vI v1-1//0. CB²
implying that R S
(1-m)Ae M
Then H
/L(O]
< R implies that
Vm
M M L(0)
m =
1)
X
39

---

Published as a conference paper at ICLR 2026
From Lemma D.4 for the estimation of L(0), i.e.
n1 + 122
L(0) S df log
8
we can deduce that
ds
mdd"
Th + 12
Ω
log
6
log
1)) Xp
8
(2) When a satísfies Assumption 4.3, we have that
v
n
CB²
V [TO)
R
R
R.
d
vmdo
From Lemma D.4, we can deduce that
m =
ds
log²
2
(dd)
log
n1 + 122
1 - A
8
Proof of Corollary C.3. Similar as before, when R S R and R"
<
we have minn (J(t)) 2
v and then Amin H())) 2 An for L = =0...,, k.
3
(t)))
Let u(() = h(t))
then
(a(l + 1) - (t)
u (wr(1) J()) ((()))(((((()))) ve(w(())
L. (u((w(s)
Jn (t)² (t) (u(())))(()(t)) ds
dw
L @u(w(t))
(84)
(t)² H(t) uu(w(t)) ds
Ow )u(w(()) )u(((s)) [n) Htt (t))))((((tt)) ds
+
Dw
dw
(() I2((),
where the second equality is from the fundamental theorem of calculus and w(s) = swo(1 + 1) + (1 - (L) = (1) J() H(1--1u(1))
Note that (t)) J (1), thus [() (t)( Plugging this into (84) yields that Ow
u(t + 1) (1 ())(() + (t)..
(85)
Therefore, it remains only to bound [2)
[a(1))2 1."" )u((w()) @u(w(s))
nn H() (u((((()))) do
Ow
Iw
/'' [[(((())
$ [|) (t) "ll""lll((((())l / [J(w(t)) - J (wo(s)
(86)
VV((t)
[((w(()) ((w
Vxo
VL(()
(IJ(w(t) J(0))l2 + ([(((ss) ((0)l2)d
V
n L(t)
R"
40

---

Published as a conference paper at ICLR 2026
where the last inequality follows from the fact that
(s) w. (0)l22 S sl|w, (4 + 1) ww (0)l2 + (1 (2,(((1))))) w, S Ri s R
and Lemma 4.5.
Plugging (85) into the recursion formula (84) yields that
|u((t + 1) [|(1 u)) + [2)))
(1 7)) (t) + (1) 2((1 ()) [2))
S (1 1)) + 2 +(111 ((1-1))
17) n))l((()) (t) 12 [2())l2
(R 2(1 17) Rn
1))
||u(t)
A
Vo
where C' is a universal constant.
Then we can choose R² such that
Cnt)(t)
(2012 S Vx
S
=
where C' is a universal constant and C1 is a constant to be determined. Thus, we can deduce that
|/ 1)112 < [(1 7)2 (C1n))++(( r)] [1((
[(1 - 7) + c 2 - C) + 7 - 1)] (1))1
< (1 - 17) (t)
where in the last incquality is duc that we can choose C1 such that n0 + 7) - 0.
Note that since 77 E (0, 1), the quadratic equation 7n²² + 2(1 77 + 17 - 1 = 0 has one negative root and one positive root, denoted as zo and respectively. Therefore, the condition C₁ S is sufficient to satisfy the requirement. The explicit form 21 can be written as:
2n 1) + 4(1 n)² 4n(1- 1
VI- 7 1+ /I
2m
Thus, C1 is sufficient to satisfy that n + 21) +
150.
From this, we can deduce that
R" S'SOIA S VI 7-
Therefore, conclude that (1n( S 7) (1))11 holds for f 0,.... k.
C.4 PROOF OF COROLLARY 4.9
Proof. In the proof of Theorem 4.7, we have proved that Condition holds for all t E N. Thus, is sufficient prove that Condition can lead to the conclusion in Corollary 4.9.
Setting 7 = 1 in (85) yields that
ul1 1) = II(().
We have that
[t)(2) VL(()
Voo S' ([]((((t))) J(w(s)))ldds..
(87)
Since wrs)) = swit + 1) + (1 s))((),, then for any E [m]] we have [s)l2 s s sllw, + 1)]|2 + s)]lm,, B.
41

---

Published as a conference paper at ICLR 2026
When (-) is smooth, can deduce that for any [m],
Ös (sp(w(s))
0sp(wo(t)
1
2
B B
(t)l)l2 I (s) w. B2+1)lwo,ss)
(B²+1) (t+1)) w. (t)ll2
Ow
ow,
2
mmm
Vnim
and
Th (w(s))
ahz(w((1)
I (B+1) (++1))
ow,
Ow,
7m
Vimm
We know that for any E [m].
B
Vmno L(t).
Thus for any E [0, 1), we have
J (mn(s) ((w(())
Üsp (wr())
2 + Ösp (sp(w(())
ah (h)(()()))
(h,(2r(t)
Ow,
Iw.
ow,
Brr
(1+* + 1) -,
(t)13 + (B² + +1) 1), 1) 20 (t)12)
Vmd L(()) )
B²
$ B²
Plugging this into (87), we have
/L(t)
/ [[((v(() J (ur()))ll2ds
Vo
VL(t) Ba
v L(I)
VAo
Vmdo
B²
L(t).
Vm
Combining with the fact ult 1) = (t) yields that
( s(1 (s((+1)) +1
CB
(0(()
((() (((1))
h( (n((+1))
S
mx
holds for E N, where C' is universal constant.
In the proof above, we only require that R' < R and R" = CdR < v3Aa leading to the requirement for that
d
md
M1 + 112
m =
13
log²
log
8
D AUXILIARY LEMMAS
Lemma D.1 (Theorem 3.1 Kuchibhotla Chakrabortty (202 If X1. Xn are independent
mean zero random variables with [[X,lls. < 3o for all < S n and some > O, then for any
vector a (a1, .nn) E TR.", the following holds true:
X 2 2cC(a))bl2vt 2eI; (et)))bllee+22,;))))))))121 (Bblst))) s 2t for all t 2 0,
42

---

Published as a conference paper at ICLR 2026
where b (a1l|Xillvar lilu.
R",
C(a) == 1ooxfv2,22/// (112 /a)
if 1,
4e + 2 2) ifa21.
for ((a) = oc when a S and B(a) = a/(a 1) when a 1
(a))) = 41 [C(a), i if
V2
4e,
In the following, we will provide some preliminary information about Orlicz norms.
Let f [0, oo) + [0, oo) be a non non-decreasing function with f(0) = O. The f-Orlicz norm of real-valued random variable X is given by
| == inf(C> > D:
$ 1}.
C
If [| < oc, we say that X is sub-Weibull of order a > O, where
(x)) Ex
Note that when a 1, || ll., is a norm and when
< 1, 11 11.s. is a quasi-norm. Moreover,
since (|a) 16) I |a|a + [b]- holds for any a, b E R and 0 < < 1, we can deduce that
21x 1/2
1/2
1 C1"
HVI 10p0
1x|a 1719 10 1eoelet elcl"
2111 1010
Ee
Ee
Re
Re 101"
Re
This implies that
|X + Ylo. I 21/a max (X 1l.. s (1I + [p.))
Furthermore, for p.q > 0, we have || = ||| x 1//0 And in the related proofs, we may 1%
frequently use the fact that for real-valued random variable N 1), we have || llv I /6 and || |1.s 1| S
Lemma D.2. If || [Yllee < oo with a, B > 0, then we have |XYI,,, S || || where y satisfies that
1
1 1
7 α
3
Proof. Without loss of generality, we can assume that || ||lles = 1. To prove this, let us use Young's inequality, which states that
I2 +
ry
Sor x, y 2 0.p.9
p
9
Let p= a/7,9 = -//7, then
|
X 7p
Y 79
E [lexp((XXYP))
< E
exp
[X|
=
E
exp
exp
P
eep((
XX (T)
exp((x13)
< E
p
9
<
2,
where the first and second inequality follow from Young's inequality. From this, we have that
|YYl,, S [1| ||
43

---

Published as a conference paper at ICLR 2026
Lemma D.3 (Bernstein inequality, Theorem 3.1.7 in Giné & Nickl (2021)) ) Let X, 1 i v be independent centered random variables a.s. bounded by c oc in absolute value. Set o 1/n Σ E and S,, = =1/n> Ei=1 X5. Then, for all 0.
2
ct
Sn
SE
n
3n
Lemma D.4. For 0 < < 1, with probability at least 1 - ò, we have that when m 2 log" (=
I(0) (#:0)) -0 log
nh
+ 112
8
Proof. Recall that for p E [n1].
sp(0)
1
m dr
(o'
(0), wp)o o" (w.(0)" 2,,) (0)]l) f(xp)
m
and for j E [122]
h₂(0)
1
Σ
a o (w, (0)) - ((v;)
T22
m
Then
n2
1 (h,
L(O)=
(0))²
Σ
2
f=1
1
R1
(o' (0)² xx)0o00))) (000))-o"'w,00000000
"(w,(0)" (00)(0))=pp)llea)) 2,)0)),)))))
2
n1 Σ fzpp)
ar
11
171
p=1
ng
n2
a (,)
g (3;).
112
m
n12 j=1
Note that
(0) - o" and (w, (0)" y;)
Jar (o' (w, xp)wro
(0,,(0))=p))l ,,,)1,,,100013)) $ (0,,
< w,, (0)l w. (0)
Since ||| (0)] (d) and (0) y pl (0) mpllv 0(1), from Lemma D.2, we have that
I0, [0, w. (0' (0"" O(d), /w, (0)" y,lll = O(d).
Applying Lemma D.I yields that for fixed p E [ni] and j E [na] with probability at least 2e -f,
m
ar (o' (00 - a" (ww0))00"""00""""""""0" (0133) S <dvi+ d
m
mm
and with probability at least 1 - -t
d
Σ a a (w, y))
S ddvi
/
44

---

Published as a conference paper at ICLR 2026
Then taking union bound for all p E [n1] and j E [n2] with 2(n) n2)e = 8 yields that
((0) 2 dvi +
Vm 1.3
d2s
2 d +
d log
R1
+ 122
log
n1 +
8
6
+ 1/2
x d² log
6
since m 2 logg (m)=te)
45


Reviews:
[{"review_id": "fYmH3Ti05L", "reviewer": "Reviewer_h7cP", "summary": "This paper theoretically analyzes the global convergence of natural gradient descent (NGD) for training overparameterized physics-informed neural networks. The motivation is strong since NGD enjoys faster convergence, as shown in experiments, compared with SGD and Adam. Compared with several previous relevant works, the paper adopts advanced concentration inequalities and a loss decomposition method to further improve the bound.", "questions": "Do you think NGD is the next optimizer for PINNs? Can you discuss your insights?", "limitations": "", "rating": 6, "confidence": 5, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "(1) The presentation is clear and easy to follow.\n\n(2) Experiments show the superior performance of NGD over L-BFGS and Adam, which motivates the study in this paper on the convergence of NGD.\n\n(3) The theory is correct and informative.", "weaknesses": "(1) Can you extend the work to stochastic versions of algorithms? I think there are some works extending GD for PINNs to SGD. Can you similarly do the extension? I think the concern for NGD is its size dependence on data samples, which is unfavorable. Therefore, the convergence (in expectation) for stochastic algorithms may be more attractive.\n\n(2) A previous work (The Challenges of the Nonlinear Regime for Physics-Informed Neural Networks; NeurIPS 2024) shows that the Newton’s method enjoys the fast convergence as in your theorem, where the convergence rate is independent of the least eigenvalues of the Gram matrix. I think the algorithm is quite similar to yours, with (JJ^T)^{-1}. Can you comment and discuss on this paper? What are the advantages of your algorithm and analysis, compared to this paper?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}, {"review_id": "TbrbKDpVMv", "reviewer": "Reviewer_TUZa", "summary": "The paper studies optimization dynamics for over-parameterized two-layer PINNs trained on a class of second-order linear PDEs in the NTK regime. \n\nFirst, it refines NTK-style analyses for gradient descent (GD) on PINNs by: \n\n- deriving a new residual recursion that improves the admissible stepsize from $\\mathcal{O}(\\lambda_0)$ to $\\mathcal{O}(1/\\|H_\\infty\\|_2)$\n\n- weakening the network-width requirements via concentration for sub-Weibull variables, yielding linear convergence with rate $(1-\\eta\\lambda_0/2)$ (Theorem~3.7). \n\n\nSecond, it establishes that full-batch Natural Gradient Descent (NGD) achieves convergence with step-size independent of the NTK kernel which could be very ill-conditioned as sample size increases.\n\nExperiments on Poisson, Heat, and Helmholtz equations (1D--10D) compare NGD to SGD, Adam, and L-BFGS and illustrate faster loss decay and lower relative $L^2$ error.", "questions": "- Line 206: \"In contrast, on one hand, our conclusion is independent of n1 and n2\". I am not sure I understand this, isn't the statement dependent  through m (logarithmic dependence).\n\n- In remark 3.8 it is said that the trace H^{\\infty} is independent of the sample size. However H^{\\infty}  is a matrix of size n time n, so a priori its trace depends on the sample size, unless it has a particular structure. Could the authors provide more evidence for why would such quantity be independent of sample size. \n- In line 262, it is said: \"we fix the output weight a and update the hidden weights via NGD.\" Why fixing the last layer and only update the hidden weights? This seems unusual? \n\n- In line 320 on the approach of Zhang et al. 2029, it is said that: \"However, this approach is not applicable\nto PINNs, because the loss function involves derivatives.\" Can the authors elaborate on why this is not applicable here?\n\n\n- In line 1823 of the proof of lemma 4.6, there is a reference to equation 46, but I don't see how this equation leads to a control on the probability of the weights being larger than M. Perhaps a typo in cross referencing? \n\n- In line 1865, did the authors mean lemma D.1 instead of D.4? Actually, in the statement of lemma D.1, the paper introduces two objects L_n and L_n* but uses only L_n*. This is confusing especially since the expression of L_n* simplifies a lot when expressing it directly. I suggest the authors directly define L_n* explicitly.  \n\n\n- In line 1011 of the proof of lemma 3.3, I don't understand why the second term was discarded, it decays faster in terms of m, but not in terms of log(1/delta) (in the limit when delta goes to 0).", "limitations": "", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "**Significance.** \nPINNs are widely used yet optimization pathologies are common; providing conditions under which GD can use practical step sizes and NGD can enjoy $\\mathcal{O}(1)$ step sizes---and even quadratic behavior for smooth activations---can influence both theory and practice in scientific machine learning. I highly appreciate this work. \n\n**Originality.**\nThe work extends NTK-based optimization theory from regression to PINNs with PDE with improved rates for GD compared to prior work (Gao et. al. 2023) and provides the first analysis of NGD for PINNs that provably shows its effectiveness. \n\n\n**Quality.**\nThe technical development is careful: improved GD bounds come from a sharper recursion (Lemma B.1) and stability of $H(k)$ via sub-Weibull concentration; NGD uses Jacobian stability tailored to individual neurons (Lemma 4.6) rather than global matrix stability. \nTheorems~3.7 and 4.7 and Corollaries 4.9--4.10 are stated with explicit dependences on d, $\\lambda_0$, and $m$.\n\n**Clarity.** \nThe paper is well organized, with explicit assumptions, side-by-side remarks comparing to prior work, and worked proofs in appendices. \nFigures and tables (p. 9--17) clearly support the claims.", "weaknesses": "**Comparison with existing results**: Although the paper compares the results in terms on conditions on the number of neurons, sample size, etc. It particular, the authors should explain why the approach in Zhang et al 2019 for natural gradient for regression is not applicable here as well. In particular, how the dependence on the derivatives makes such approach inapplicable. \n\n**Experimental evaluation**: Some information is missing. For instance, the authors should clarify what they mean by the L^2 error, is it the error on the training samples (colocation points), or do they consider test samples. In case, they evaluate on training samples, that's fine since the results are for the training, but it would be also interesting to evaluate on a test samples to get a better sense of generalization.\n\n\n**Minor**: The abstract contains notations (lambda_0 and H^{\\infty}) that are not defined. That might give a bad first impression to the reader, please fix this.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}, {"review_id": "BqtfxNIsiR", "reviewer": "Reviewer_jYGf", "summary": "The manuscript provides training guarantees for physics-informed neural networks in an overparametrized setting with NTK-scale parametrization, both for first-order gradient descent as well as a second-order natural gradient descent. For gradient descent, the allowed step size of $O(1/\\lambda_\\max)$ is different compared to prior works with $O(\\lambda_\\min)$ step size on optimization in PINNs. For natural gradient descent, a convergence rate independent of the condition number of the problem is derived. Additionally, computational experiments comparing first-order to second-order optimizers in PINNs are provided.", "questions": "+ Can you elaborate why you regard your $O(1/\\lambda_\\max)$ as an improvement over $O(\\lambda_\\min)$? \n+ Can you elaborate on the relation of the natural gradient defined in Section 4 to previously proposed methods like energy natural gradient? In case of equivalence, can you adapt your presentation accordingly, see also my discussion above. \n+ Can you comment why you make the overparametrization assumption? Also, a discussion of the overparametrization assumption inside the manuscript seems reasonable. \n+ Can you comment on the different condition numbers of supervised learning and PINNs? Can you elaborate why you only expect a bad condition number for complex PDEs (see line 251)? This could be added in the manuscript as a motivation for second order and natural gradient methods. \n+ Can you validate the theoretical convergence guarantees within the computational experiments? In particular: \n   + Compare the loss curves to the linear convergence guarantees from Theorem 3.7 and 4.7. \n   + Compare the generalization error for different network size and optimizers. As the convergence guarantees are merely given with respect to the loss function, it is informative to contrast it with the relative L^2 error. \n   + Improve readability of the plot. \n\n**Minor comments:**\n1. In the abstract, it is stated that, *However, the learning rate of GD for training two-layer neural networks exhibits poor dependence on the sample size and the Gram matrix*. Note that the term learning rate usually refers to the step size in machine learning, not the convergence rate of the learning process. \n2. It is stated in several places that the learning rate can be O(1). However, it is shown that the convergence rate is O(1). \n3. In the limitations section on scalability of natural gradient methods in these scenarios, the following directly relevant works are not referenced: \n+ Kronecker-Factored Approximate Curvature for Physics-Informed Neural Networks, Dangel et al. 2024\n+ Improving Energy Natural Gradient Descent through Woodbury, Momentum, and Randomization, Dangel et al. 2025", "limitations": "", "rating": 4, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "+ Overall, the paper is well written and easy to navigate and the main contributions are clearly described. \n+ The paper addresses a timely and important question of the convergence of optimization physics-informed models. \n+ The guarantees go beyond the use of first-order methods and treat natural gradient, which seems to be the arguably most efficient optimizer in PINNs. \n+ The theoretical analysis uses timely methods from supervised learning and rigorously transfers them to PINNs.", "weaknesses": "+ Discussion of related works: \n    + Use of *our NGD*: At multiple places in the manuscript, *our NGD* is used to refer to the natural gradient method that is being analyzed. Also, in Remark 4.2 it is claimed that the method studied is different from previously proposed natural gradient or Gauss-Newton methods for PINNs. The iteration of energy natural gradient (Müller and Zeinhofer, 2023), which reduces to a Gauss-Newton iteration for PINNs, is given by $$ w(k+1)=w(k)-\\eta (J(k)^T J(k))^+J(k)^T \\binom{s(k)}{h(k)}, $$ where $A^+$ denotes an arbitrary pseudo-inverse of $A$. By choosing the Moore-Penrose inverse and considering an overparametrized setting in which $J(k)J(k)^T$ has full rank, this specializes to $$w(k+1)=w(k)-\\eta J(k)^T (J(k)J(k)^T)^{-1} \\binom{s(k)}{h(k)},$$ which agrees with the iteration studied in the manuscript. Hence, from my understanding, the manuscript is studying the convergence properties of the optimizer proposed by Müller and Zeinhofer (2023), however, this is not mentioned at any place. Rather, the impression is conveyed that the manuscript proposed a new variant of a natural gradient method, see Remark 4.2 as well as line 408 and Table 1. This \n   + There are certain works on optimization in PINNs missing in the discussion of related works in Subsection 1.2., e.g., \n      + AN OPERATOR PRECONDITIONING PERSPECTIVE ON TRAINING IN PHYSICS-INFORMED MACHINE LEARNING, Tim De Ryck, Florent Bonnet, Siddhartha Mishra, Emmanuel de Bézenac; in particular, Theorem 2.3 shows that for the linearization problem with preconditioning achieves a converge rate given by the condition number and that the natural gradient preconditioning achieves an optimal condition number of 1. As such, this seems to be a linearized (hence easier) version of the main result in the manuscript. \n      + Convergence of Stochastic Gradient Methods for Wide Two-Layer Physics-Informed Neural Networks, BANGTI JIN AND LONGJUN WU, 2025; provide an exponential convergence rate of SGD under overparametrization\n      + Non-Asymptotic Analysis of Projected Gradient Descent for Physics-Informed Neural Networks, JONAS NIESEN, AND JOHANNES MÜLLER, 2025; provide a sublinear convergence guarantee for arbitrary sampled (S)GD without assumption on the network size \n+ Overparametrization assumption: The manuscript is making the global assumption of overparametrization. Where I understand that this allows the use of an established machinery, I do not believe that this is the setting that PINNs are used in. In particular, note that in PINNs, the data points are synthetically sampled integration points of the computational domain rather than data points like in supervised learning. As such, the problem in PINNs is rather an optimization than a statistical problem. In practice, new data points are sampled continuously throughout optimization. Hence, it is not clear how practically relevant the setting of overparametrization is. However, this assumption is not uncommon and the only work I am aware not making an assumption on the size of the PINN is by Niessen and Müller (2025). \n+ Difference of PINNs to supervised learning: In the introduction, it is mentioned that the reason why NGD is not used in supervised learning is due to its high computational cost. However, I believe that there another, arguably more important reason: The loss function of PINNs has a much worse conditionining due to the appearence of the PDE operator. It is stated in line 251 that the conditioning can be really bad for complex PDEs. From my understanding, this can also be the case for simple PDEs. \n+ Experiments: The experiments compare natural gradient to SGD, Adam and L-BFGS. However, this comparison has been made at several places, in particular, by Müller and Zeinhofer (2023). I think, it would be much more informative to not repeat this comparison, but to see, how the empirically observed convergence rates relate to the theoretical guarnatees. Further, in relation to the overparametrization assumption, the influence on the network size and optimizer on the generalization error is not studied empirically.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}, {"review_id": "3PHdMYI29X", "reviewer": "Reviewer_yT77", "summary": "Authors provide a proof for the convergence rate of a certain type of neural network that depends on the norm of the Hessian, which seems to be a new result (or an extension of a known result to a broader class of problems). The paper is clear and well-written, and overall seems like a good contribution to the field of theoretical guarantees for neural network training.", "questions": "See weaknesses.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The paper is clear, of a high-technical quality and of a fair significance. It is an extension of a result known for regression problem to more general PINNs (although still restricted to two-layer networks).", "weaknesses": "A few things come to mind:\n1) The type of network it is applied to is quite restrictive (two-layer ReLU PINNs). While I understand that this allows you to construct the proof, could there not be a way to extend the result to more general architectures? This would make the result and impact much stronger.\n2) It would be good to numerically validate the main result of the paper, where you take a simple problem where you can calculate $1/|H^{\\infty}|_2$ exactly and verify that the convergence rate matches it.\n3) A natural question that comes to mind is what happens when you have deeper networks. It would be interesting to at least numerically study this question, by replicating your experiments with three and four-layer networks perhaps.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)