Dual Natural Gradient Descent for Scalable Training of Physics-Informed Neural Networks
Authors: Anas Jnini
Flavio Vella
Venue: Accepted by TMLR ()

Abstract: Natural-gradient methods markedly accelerate the training of Physics-Informed Neural Networks (PINNs), yet their Gauss–Newton update must normally be solved in the parameter space, incurring a prohibitive $\mathcal{O}(n^{3})$ time complexity, where $n$ is the number of network weights. We show that exactly the same step can instead be formulated in a generally smaller residual space of size $m=\sum_{\gamma}N_{\gamma}d_{\gamma}$, where each residual class $\gamma$ (e.g. PDE interior, boundary, initial data) contributes $N_{\gamma}$ collocation points of output dimension $d_{\gamma}$.

Building on this insight, we introduce Dual Natural Gradient Descent (D-NGD). D-NGD computes the Gauss–Newton step in residual space, augments it with a geodesic-acceleration correction at negligible extra cost, and provides both a dense direct solver for modest $m$ and a Nyström-preconditioned conjugate-gradient solver for larger $m$.

Experimentally, D-NGD scales second-order PINN optimisation to networks with up to 12.8 million parameters, delivers one- to three-order-of-magnitude lower final $L^{2}$ error than first-order (Adam, SGD) and quasi-Newton methods, and—crucially—enables full natural-gradient training of PINNs at this scale on a single GPU.

Published in Transactions on Machine Learning Research (10/2025)
Dual Natural Gradient Descent for Scalable Training of Physics-Informed Neural Networks
Anas Jnini
anas.jninn@emmiinn.t
Deportment of Information Engineering and Computer Science University of Trento, Trento, Italy
Flavio Vella
flavio, vella @unitn it
Deportment of Information Engineering and Computer Science University of Trendo, Trento, Italy
Reviewed on OpenReview: https: // owtpr/veeerrvvee...mm net/ formm? id= GDHVRy6SDd
Abstract
Natural-gradient methods markedly accelerate the training of Physics-Informed Neural Net- works (PINNs), yet their Gauss-Newton update must be solved in the parameter space, incurring a prohibitive Onn)) time complexity, where is the number of network trainable weights. We show that exactly the same step can instead be formulated in generally smaller residual space of size m = Ey Nydy, where each residual class (c.g. PDE interior, boundary, initial data) contribules Nr collocation points of output dimension dy.
Building on this insight, we introduce Dual Natural Gradient Descent (D-NGD). D-NGD computes the Hlass-Newtt step in residual space, augments it with geodesic acceleeraooooooii correction at negligible extra cost, and provides both a dense direct solver for modest m and Experimentally, D-NGD scales second-order PINN optimization to networks with up to 12.8 million parameters, delivers one- ofroeorrde-oo--maannuuee lower final error I than first-order (Adam, SGD) and quasi-Newton methods, and crucially enables full natural gradient training of PINNs at this scale on single GPU.
conjugate-gradee gradient solver for larger m.
1 Introduction
Partial Differential Equations (PDEs) Partial Differential Equations (PDEs) form the backbone of mathematical models used to describe wide array of physical phenomena-ramme from fluid flow and heat transfer to the behavior of advanced materials. Conventional discretization-basse techniques, such as finite element and spectral methods, often demand highly refined meshes or basis expansions to attain the desired level of accuracy. This refinement drives up computational costs, especially in engineering scenarios that call for numerous simulations with varying boundary conditions or parameter sets. In recent years, machine learning approaches- most notably those employing neural networks have emerged as promising complement substitute for these traditional solvers, offering potential gains in efficiency and flexibility
Raissi al. (2099 Li et al. 2021 Jnini & Vella
Inforie- Neural Networks (PINNs) PINNs are machine learning tool to solve forward and inverse problems involving partial differential equations (PDEs) using neural network ansatz. They have been proposed as early as Dissanayake & Phan-Thiee 1994) and were later popularized by the works Raissi et al. (2019); Karniadakis et al. (2021). PINNs are a meshfree method designed for the seamless integration data and physics. Applications include fluid dynamics Cai al. (2021): Jnini et al. (2025a b). solid mechanics Haghighat al. (2021) and high-dimensional PDEs Hu (2023) name but a few areas of ongoing research.

---

Published in Transactions on Machine Learning Research (10/2025)
Despite their popularity, PINNs are notoriously difficult to optimize Wang et al (020 and fail to provide satisfactory accuracy when trained with first-order methods, even for simple problems Zeng et al. (2022)
Müller & Zeinhofer (2023). Recently, second-order methods that use the function space geometry to design
gradient preconditioners have shown remarkable promise in addressing the training difficultics of PINNs Zeng
et al. )022); Müller & Zeinhofer (2023); Ryck et al.] (2024); Jnini et al. (2024p; Müller Zeinholer (2024)
(GN) principles, can harness curvature information for improved convergence, their canonical forms entail (nn") per-iteration complexity and Onn²)
While second-order optimizers, such as those based on
memory for n parameters, rendering them impractical for large-scale deep neural networks. Matrix-free methods have been proposed to compute Gauss-Newton directions without explicitly forming the Hessian
Martens 2010); Schraudolphi (2002); Zeng et al. (2022); Jnini et al. (2024). Despite reducing computational costs, these methods suffer from ill-conditioning, leading to slow convergence for large nelworks without efficient preconditioners.
To address these challenges, this work proposes Dual Natural Gradient Descent, a novel optimization framework incorporating the following key features:
We propose a primal-dual viewpoint of the Gauss -eeusn step: instead of solving the usual n-dimensional normal equations in parameter space, we work with their dual in the m- dimensional residual feature space, where m = Σ, N,dy, with each residual class y (e.g. PDE interior, boundary, initial data) contributing N, collocation points of dimension dy, and since m & n in practical PINNs this dual system is far cheaper to assemble, store, and solve than the original paramoter-space system.
We propose an efficient geodesic-acceleration correction within the same dual frame- work. A second-order (geodesic) term is obtained by solving one additional linear system with the same left-hand side and reuses the factorisation already built for the primary dual operator, adding negligible overhead while improving step quality.
We propose a low-rank Nyström spectral preconditioner for large batch sizes. When m is too large for direct factorisation, we propose using a Hessian-Free iterative method with a preconditioned conjugate gradient, we provide an efficient Nyström preconditioner based on column sampling for the dual problem.
We demonstrate scalability and accuracy on several PDE benchmarks. The resulting dual natural-gradient method trains PINNs with tens of millions of parameters and dimensions in the hundreds of thousands, consistently outperforming first-order and quasi-Newton baselines by more Ihan an order of magnitude across several representative problems. To the best of our knowledge, our contribution the first to extend Natural- Natdret-Gradient methods to PINNs of this scale.
Related Works
Second-order Optimization in PINNs The challenge of elfectively training PINNs has spurred signif- icant research, with growing consensus underscoring the necessity of second-order optimization methods. Recent literature highlights this trend: approaches leveraging an infiniteddmmensional perspective, for in- stance, have demonstrated the potential to achieve near single-precision accuracy Zeng et al. (2022 Müller & Zeinhofer (2023): Ryck et al. (2024); Jnini et al. (0224); Zampini al. (2024); Schwencke & Furtlehner (2024); |Mckay et al. (2025). However, the practical application of these methods is often constrained by their high per-iteration cubic computational cost, particularly when scaling to larger network architectures due to the need to solve substantial linear systems in the parameter space. Exploration into more scalable alternatives includes quasi-Newton methods, which Kiyani et al. (2025) evaluate for their efficiency and accuracy across stiff and non-linear PDEs by leveraging historical gradient information. Complementing this, Wang et al. (2025) provide a theoretical framework for gradient alignment in multi-objective PINN training, demonstrating how sccond-order information, through Hessian preconditioning, can resolve direc- tional conflicts between different loss components, thereby demanding more sophisticated curvature-aware approaches.
2

---

Published in Transactions on Machine Learning Research (10/2025)
Hessian-Free Curvature Approximation Matrix-free methods have been proposed to compute Gauss- Newton directions without explicitly forming the Hessian Martens (2010); Schraudolph (2002 Zeng et al. (2022]; Jnini et al. (2024). Despite reducing computational costs, these methods suffer from ill ill-conditionng, leading to slow convergence for large networks without efficient preconditioning JJnin et al. (2024) Our algorithm addresses this by proposing an efficient preconditioner for the dual syslem, signilicantly improving the inner solver convergence. This idea of cutting large gaps within the leading eigenvalues of the Hes- sian spectrum is also aligned with recent advances in preconditioning techniques, such volume sampling Rodomanov & Kropotov (2020g, polynomial preconditioning Doikov & Rodomanov (2023). and spectral preconditioning Doikov et al. h0024: Frangella et al. (002)).
Efficient Second-Order Methods via Primal-Dual Formalism Our work introduces primal-dual formalism to scalably solve the regularized Gauss-Newton (GN) least-squares problem that arises in PINN training. This primal-dual viewpoint, which moves the main computation from the large parameter space to smaller dual space, has been previously linked in other fields to the push-through identity Henderson & Scarle (1981). For instance, this connection is foundational in kernel methods and Gaussian Processes, where it is the basis for the kernel trick. Despite its effectiveness, this approach has remained largely unexplored within the Natural Gradient Descent (NGD) and, in particular, the PINN communities. Related algebraic tochniques have appearod, for example in the work. Benzing (2022), to create a scalable step for the Fisher Information Matrix.
Geodesic Acceleration in Natural Gradient Methods Beyond the natural gradient step, higher-order corrections can improve convergence speed and step quality. The concept of using a geodesic acceleration term, which accounts for the curvature the manifold along the gradient path, was notably explored by [Song iet al. (2018); Bonfanti et al. (2024) and can be interpreted as an approxination of the Riemannian Euler method. Our primary contribution in this context is to demonstrate that this acceleration term can be computed with negligible overhead within our dual framework. By reusing the existing factorization of the dual operator for a second linear solve, we make the geodesic correction practical and scalable, integrating it seamlessly into our optimization workflow without significant performance penalty.
Connection to Operator Learning While our work focuses single-instance PDE solutions learned through PINNs, the scalability of D-NGD is also highly relevant to the more general task of operator learning. Architectures like Physics-Informes Neural Operators (PINOs) Li et al. (2021) and PIDeepONets Wang et al. (2021h learn mappings between function spaces and often require models with vast parameter dimensions (n), rendering traditional O(n3) second-order methods impractical. Crucially, these frameworks incorporate a PDE residual term, creating a nonlinear least-squares objective analogous to that of standard PINNs. Because D-NGD's complexity scales with the residual dimension (m) rather than the parameter count (n), it provides practical pathway apply curvature-awar, optimization these large-scale models, a domain that has thus far been largely reliant first-order methods.
2 Preliminaries
2.1
Physics-Informed Neural Networks
For a given domain Ω Igd (or Or 1 x Ω for time-dependent problems, where I time interval), consider general PDE of the form
Lu = in Ω
subject initial and boundary conditions, collectively denotod as = 9 on 0Ω (where Nn here generally refers to the spatio-temporal boundary). PINNs approximate the solution w of the PDE using a neural network ansatz us, parameterized by 0. The loss function is defined as:
No
Nos Σ
L(0)- = 2N
Σ n=l
(Cuo(En) (((n)) + 2N n=l
(uo(Ix) - 2)
(1)
3

---

Published in Transactions on Machine Learning Research (10/2025)
where {En, E R are the collocation points in the interior of the domain and {ax E @0) are the points which initial and boundary conditions are enforced.
2.2 Gauss-Newton method for PINNs
Residuals and Jacobian. Let i ,,.... No and I,.... Non. Define the discrete residual map to
r: (0) + R to be
r(0) (ra(0); ron E R""
m = Nada + Nondan,
with
ra) VNg
(L f(x:)) E Rd0, ran(O)) VNon (m0(2))-g(z))
g(x5)) E Rdoo
Its Jacobian is [(0) 0gr(0) Rmxx and the loss is
L( (0))1
Gauss Newton natural gradient descent. First-order optimizers, such as gradient descent and Adam, often fail to provide satisfactory results due to the ill-conditioning and non- convexitx of the loss landscape, as well as the complexities introduced by the differential operator C Wang et al. (2020). Instead, function-space-
inspired second-order mothods have lately shown promising results Zeng ct al. (2022) For the remainder of this paper, we adopt Gauss Newlon Natural Gradient Descent (GNNG) Jnini et al. (2024). Linearizing the residual map in function space and pulling the resulting Gauss-Newtr operator onto the tangent space of the ansatz yields the Gauss-Newer Gramian:
G(0) = J(B)
(2)
and update rules at iteration A:
Dut1 = 02 G(Ox)) (OL(
k
(3)
It was shown that the Gauss Newton direction in function space corresponds exactly to the Gauss Newton
step for L(0)= = 11) parameter space when the same quadrature points are employed Jnini et al. (024 g Müller & Zeinhofer (2024).
Least squares characterisation.
Equivalently, the Newton increment is the solution of the linearized problem
~ arg min ABER" [[(r())) ((xx) A0
(4)
whose first-order optimality conditions give the normal equations
T(O) J(Ok) Ao = J(Ox) r(
(5)
Given that VoI((e)) = J((B))Tr(O.)) we recover the update equation 3
3 Optimization in the Residual Space
In this section we propose a primal- viewpoint on update can be cast either as an n n linear system in parameter space or as an m m system residual space (where in typical cases m & n), allowing the residual-space cestem-white can be solved by dense factorization matrix-free iterative method-tee scale to problems with large batch sizes. In Section 3.3 detail the iterative solver and the low-rank Nyström preconditioner that accelerates it. We incorporate the classic Levenberg Marquardt damping À > 0 Levenberg (1944); Marquardt (1963) into the Gauss Newlon model equation defining the Tikhonov- problem
updates for PINNs: the same parameter
min FA(40)) (1(,(0))) + J(Ow) A0M1 + ||00.13 A0
4

---

Published in Transactions on Machine Learning Research (10/2025)
Proposition 3.1 (Primal normal equations). For any parometer vector 8k E IR", define the primal Gramian G(Ox) J(O))" J(Oz).
The unique minimiser 20g E I of the damped least-squares subproblem then salisfies
(G(Oa) + AI)) A0% = VoL((k).
Here, VoL(Ok) ((0)) (((0)) is the gradient of the unwwigghed least-squares loss g1lr(0=)13
Proof Outline: The proof follows by standard differentiation of the objective function F'(40s) with respect to A0, and setting the result to zero. The objective strictly convex for À > 0, guarantecing a unique minimizer. The detailed derivation is provided for completeness in Appendix A.1.1
The linear system given in Proposition 3.1 corresponds to the classical Levenberg Marquardt update in the n-dimennional parameter space. Forming and factorizing the matrix J(0x)) J(0s) AI, in equation 3.1 requires O(mn² + n²) time and Onn²) memory, which becomes impractical when n is large.
3.1 Dual Formulation via KKT in Residual Space
An alternative formulation to the system equation 3.1 be derived by working in the m-dimensional residual space.
Definition 3.1 (Residual Gramian). Let the residual Gramion matriz be defined as
K& ((kk J(O%) E R
We introduce an auxiliary variable == J(as) A0; E Rm and enforce the constraint yk = J(O.. A0e using Lagrange multipliers VA E R The Lagrangian for the subproblem of minimizing F)(40)) (from Proposition 3.1 becomes
L(A0,, yk 14) = 1 (1m)())) + y* + I + UI (yk- - J(Ok) A0k).
Applying the Karush Kuhn Tucker (KKT) conditions for optimality provides the dual formulation.
Theorem 3.1 (Dual Normal Equations). Let yk := (0x) Ad, E R" be the predicted residual decrement associated with a parameter step 20g E R". Applying the KKT conditions to the minimization of F2(00)) yields the dual linear system for the optimal yt:
(K+ Am)) u = ((0s) VoL(Ou).
(6)
The parameter update can then be recovered by
A0t = - (J(Ok) Di VoL(Ox))
(7)
Conversely, any pair (Ark. Vse satisfying equation 6 and equation constitutes the unique primal-dual opti- mum of the Levenberg-Marquordt subproblem.
Proof Outline: The KKT conditions are derived by setting the partial derivatives of the Lagrangian L L(AO,, yk. v)) with respect to Abk, yk, and ia to zero. Algebraic manipulation of these conditions yields the dual system and the recovery formula for the parameter update. The full proof is available in Appendix A.1.2
Primal vs Dual Perspectives. We have derived two equivalent formulations for the LM step: the primal system in the n-dimensional parameter space, involving the n X n primal Gramian G Xl, (with G = J J),, and the dual system in the m-dimensional residual space, involving the m x m matrix K₂ + A4n (with K2 J(O,) ((x)). In practice one solves the smaller system directly-using the primal when n < (cost O(n²)) or the dual when m & (cost O(m²))-or or resorts an iterative, matrix-free Preconditioned Conjugate Gradient solver with preconditioning for large-scale problems.
5

---

Published in Transactions on Machine Learning Research (10/2025)
3.2 Geodesic Acceleration
When interpreting the Levenberg Marquardt step velocity along geodesic in parameter space, one can improve the update by adding a second-order correction that accounts for an acceleration dg along that same geodesic:
Pk-1 = 8% + Uk + ax.
where dg satisfies
(0k) an -fevs
and
d r0p + ttxx)
fve == dta
E IR
1==0
is the second directional derivative of the residual along uk. Computing for involves two Jacobian-vector products, capturing curvature beyond the standard LM step. Below we show that this geodesic-acceleration term can also be computed using our proposed primal-dus formalism.
Definition 3.2 Accileertiooo Subproblem). At iteration k, the acceleration correction a is the minimiser of the damped system
min aER"
]11](0s) a Soell + 111a
where fev is defined above.
Proposition 3.2 (Primal and Dual GA Charactorizations) Let KR ((B)) ag of equation admits both a primal and dual formulation:
The unique solution
(Primal) (J(Ok) J(On)) + AIn) a₁ (0s)) fev:
(9)
(Dual) (KCx + Al Ya.k - Kr Jons Ox = -(0x) (ya,k + Smn).
(10)
Remark 3.1 (Implementation Cost of GA). Solving for a via the dual system equation 10 reuses the same operator KC + Am as the primary LM step. Hence, once factorization or preconditioner for the dual solve is available, the only extra costs are for computing Jw (one Hessian vector product two JVPs) and solving additional linear system with the eristing operator. If using an iterative solver, the preconditioner may be shared; with direct method, only extra buck- suck-uubttitution is needed.
3.3 Hessian- Bossian-FFree Solution of the Dual System
We aim to solve the dual linear system, previously stated in equation For moderate residual dimension m, one assemble Kj and solve equation t with a dense Cholesky factorisation. When the total number of residual evaluations m (which can be thought of as batch size in this context) is very large, however, building Sonjugate-Gradient Kj explicitly is not tractable. In such cases, we turn to the Conjugate- Gradient (CG) method (Hestenes & Shewchuk 1994 Saad 2003). CG only requires repeated applications of the lincar operator 0 0 (K] Alm) and the right-hand side vector.
Stiefel 1952
Definition 3.3 (Kernel Vector Product). For v E m the action of the operalor is oblained three
difommaticdiiferenntttioo passes:
Reverse mode: u=
2. Forward mode: H = J J===JJuee
Tikhonov shift:: v0 Av.
We denote this composite map by w kvp(w).
A CG iteration therefore invokes kvp(-) once, allowing us to solve equation 6 Hessian-free: neither J nor K. is ever formed explicitly, yet the exact solution is recovered through operator evaluations alone.

---

Published in Transactions on Machine Learning Research (10/2025)
3.3.1 Low-Rank Nyström Spectral Preconditioner for the Dual System
Both the primal and residual Gramians share the same eigenvalues. In typical PINN application, they often show rapidly decaying spectra and large spread between the top and bottom eigenvalue Ryck et al. (2024). which lcads to severe ill-conditioning. To address this, we propose a spectral proconditioner. Inspired by the Nyström ideas of Frangella al. (2023); Martinsson & Tropp (2020), we approximate the kernel K through mank-runncted eigendecoomposiiion K. s UAU Here, the columns of U E Rm"x'' are approximate cigenvectors and E R x contains the corresponding approximate cigenvalues, where l' the effective rank determined by the Nyström procedure (lypically S e, the initial nunber of landmarks). We then build the operator
p-1 = - + 1 (4 UUU)
(11)
which damps the leading Nyström modes by the regularised factors (As + A) and acts as > on the orthogonal complement. We employ p-1 as left preconditioner Conjugate Gradient, one application of the preconditioner requires only two dense natrix-vector products with U or UT. Empirically, this spectral preconditioner compresses the spread of eigenvalues and reduces CG iterations.
4
Algorithmic Implementation
This section outlines the implementation of the dense and itcrative dual solvers. Detailed algorithms are provided in Appendix A.2
4.1 Dense Dual Solver
For problems where the residual dimension m is considerably smaller than the parameter dimension n (m < n), direct approach to solving the dual system (Equation equation [ is feasible. This involves explicit formation of the regularized residual Kernel, K, and Cholesky docomposition.
Residual- Kernel Assembly Let the residual vector introduced in Sestion 2.2 be
r(0) = (ra(0), TT)
m Nada + Nandan,
with Jacobian split J(0) = (Jn; Jan). The Gramian required by the dual formulations (Sections
3.1
3.3
Kon Kana
Knaq Kausn)
K= =
Kann = =KJon-
The Kernel K JJT is built on-the-fly, block-wise, and in parallel (Algs. 3 2 Vector-Jacobian products (VJPs) provide Jacobian components for IC entries; these components are then discarded, avoiding storage of the full m x Jacobian J. Through Jax Bradbury et al. (2018) just-in-time compilation, The Accelerated Lincar Algebra (XLA) compiler optimizes these on-the-fly computations; its operator fusion capabilities en- hance speed and reduce memory overhead by preventing the explicit materialization of complete intermediate Jacobian arrays before their consumption. Moreover, symmetry is exploited for diagonal blocks. The final matrix used by the dense solver is the regularized residual Gramian K2 Kr. Alm-
With the fully assembled and regularised Gramian Kg in hand, we perform a single Cholesky factorisation and two backssubttitutions to obtain both the Gauss Gauss-Newton Newton velocity and its geodesic correction. The procedure is detailed in Algorithm
Dense Dual Solver Steps
Overall cost. Assembly O(m2n) (assuming JVPs/VJJP for kernel entries), Cholesky 3m² and each triangular solve Omm?); memory is O(m²). This remains efficient for up a few thousand.
7

---

Published in Transactions on Machine Learning Research (10/2025)
4.2 Iterative Dual Solver
When the residual dimension m becomes too large for the dense solver outlined in Section 4.1 forming and factorizing the m m residual Gramian K. (Definition [.1] becomes computationally prohibitive. In such scenarios, we resort. to iterative method to solve the dual linear system presented in equation 0 Specifi- cally, we employ the Preconditioned Conjugate Gradient (PCG) algorithm. For efficient PCG convergence, we employ the Nystrôm-based spectral preconditioner detailed in Section 3.3.1 This approach relies on the approximation K. x UAU and was proposed in Arcolano (2011), the practical construction of which is outlined next. First, a small subset < m residual components are selected as landmarks. In the context of PINNs, a landmark corresponds to a specific scalar residual value, evaluated at a single col- location point and for single output dimension of the residual map. Let ff1......4 ie) denote the set of indices for these landmark residuals.
Using these landmarks, we form two key submatrices of the full Gramian K:
.KIr ReX, The Gramian matrix computed between the landmark residual components themselves.
Ker (mm-()xx; The Gramian matrix computed between the non-landmark residual components (indexed by C' [l....,m. m) \ XI) and the landmark residual components.
An eigendecomposition is performed on the (lypically small) landmark Gramian: Krr QA&QT, where Q Exx contains orthonormal eigenvectors corresponding to the positive eigenvalues in the diagonal matrix AQ E grxr (where r S l is the effective rank of K(r)..
An extended, non-orthogonal basis E R for the Nyström approximation is then constructed. The rows of Ü corresponding to the landmark indices are set to Q (i.e., Ü, + Q), while the rows corresponding to non-landmark indices are computed as ÜC + KciQAa- The Nyström approximation of the Gramian is then given by k Ü UAQUT
To obtain an orthonormal eigendecomposition U for this approximation R which is beneficial for constructing the preconditioner, we proceed as follows. Define M = -' E mmrr Performing Singular Value Decomposition (SVD) on M yields M V.WT, where V g has orthonormal columns, Σ E ggrx is diagonal with singular values, and W R is orthogonal. Then, the Nystrõm approximation can be expressed K TM (VEWT)(VEW)) () = VV Thus, the desired orthonormal eigenvectors are U + V, and the corresponding eigenvalues Å D The complete procedure for constructing U and Å is detailed in Algorithm 5
This Nyström-based spectral approximation R R UÂUT is then utilized to build the preconditioner p-1 (as defined in Equation equation II for the PCG algorithm. The full PCG-based solution of the dual system, incorporating this preconditioner and relying on kernel-vector products (Definition 3.3), is outlined
in Algorithm 0
4.3 Optimization Workflow
Given a partial differential equation (PDE) and a neural network ansatz ug. we minimize the loss equation 1 using the described Dual Natural Gradient Descent Framework. The overall procedure, including step computation via either DeNsEDUALSoINe PCGSTEP and line search, is detailed in Algorithm
8

---

Published in Transactions on Machine Learning Research (10/2025)
Algorithm Dual Natural-Gradient Descent (D-NGD) Workflow
1: Input: initial parameters 0o: residual function Tin collocation sets (Xn, Xan) loss L(0); time budget maxx Levenberg- Saaaaberr-Marqurrttt rule for dp: Nyström rank & CG tolerance E; max CG its mcG maxi dense threshold tcense_threshi thresh' flag use use_densee lee__dnse_GG___ flag.
to f- --no((
while now() In max do
9 g VoL(O) {Compute J at all points} Let A0" be the computed parameter step.
6: m + [xa] + |Xanl
7:
if m < Cdeuse thresh then
8:
A0* + DENSEEDUALSLLVE(,, 9, Tio (Na, Xon), X, use srve_dens___ee BeanseGG__se flag)
clse
10:
A < PCGSTEE(()) 9, Fin {r,}, A, l, E, mcG_may
11: end if
12:
arg (001) L(0 + n A0*)
13: f 0 +
14: end while
15:
16:
return
4.4 Time and Space Complexity
We summarize per-step time/memory for the primal, dense dual, and dual PCG solvers.
Primal (parameter-space GN/LM). The primal approach forms = JJ and factors an n system; per-step time Omm?² + n²) and memory O(n²). This is preferable when m and is small enough to factor.
Dense dual .residual-space). The dense dual approach forms K. K =J JJJ and factors per-step time O(m²n + m²) and memory O(m²). This is preferable when > and m is small enough to factor; as a rule of thumb, if the hardware can hold m² entries at the target precision, use the dense dual Cholesky, otherwise use the iterative solver below.
m system;
Dual PCG (Hessian-free). The dual PCG approach avoids forming K. Each inner iteration costs = (mm) +0 so solve with l iterations costs ((() and uses memory linear in m, n plus O(E2). The Nyström preconditioner built once per epoch and reused within that epoch's inner PCG solves. Total wall-time is the per-step cost times the number of outer steps 8, and both s and grow with the effective conditioning x (typically = O(VRe)), reduced by good preconditioning).
5 Applications
To comprehensively assess the capabilities of Dual Natural Gradient Descent (D-NGD), we conduct series of experiments across diverse suite of benchmark partial differential equations. All solvers are implemented in JAX 0.5.0, compute the derivatives using in-house implementation of th Taylor-mode Automatic- Differentiation Bettencourt et al. 2019) unless specified.
Every run is executed on a single NVIDIA A100-80GB GPU. Unless stated otherwise, we track the relative I?? error with respect to an analytic or DNS reference solution and report the median value over ten independent weight initializations. Instead counting iterations, unless specified otherwise each optimizer is allocated fixed wall-clock budget of 3000 seconds for all experiments.
We compare the following Optimizers.

---

Published in Transactions on Machine Learning Research (10/2025)
ADAM (ingmaa & Ba. 2017] with learning rate 10-3, (81, (0.9, 0.999, c 10-8, and no weight decay.
SGD with Nesterov mnmentun 0.9 and the one-cycle schedule of Smith & Topin 2018 using peak learning rate of 10-2 and initial/final learning rates of 10-4
BFGS (Jaxopt Blondel al. (2021) implementation) with history size 300, tolerance 10-6, max- imum 20 iterations per line-search step, and strong- Wolfe backtracking line search.
Dense D-NGD the Dense Qual-Newton natural gradient solved by a dense Cholesky factorisa- tion (Algorithm with and without geodesic acceleration (suffix *+GA").
PCGD-NGD the matrix-free free CG variant (Algorithm 6 preconditioned by a Nyström approx- imation (Algorithm 5
Across all benchmarks, we employ standard Multi-Layer Perceptron (MLP) architectures with tanh activation functions. This choice is motivated by primary goal to develop architecture--gnosti optimization framework and opt to use one of the most commonly used architectures in the PINN literature Wang et al. (2023). While specialized architectures can undoubtedly improve the conditioning of the loss landscape and may lead to even lower final errors Jnini et al. (2025a), our results demonstrate that a powerful, curvature- aware optimizer can achieve sfate-of-the-art accuracy even with a vanilla network structure.
For the dense NGD solver, we report results both with and without the second-order geodesic correction (the former denoted by the "+GA" sullix). All hyperparameter settings listed in Appendix A Below, we introduce PDE benchmark and experimental setups used for the selected test cases.
5.1 The 10+1-Dimensional Heat Equation
We consider heat equation in 10-dimensional spatial domain plus time, defined for E Ω == :00==00,]] 1), and E [0, 1], with diffusion coefficient 1.4 The temperature field u(t,x) evolves according to:
)u((, (100)xA-... xAzu(t,)) 0,
u(0, Σ sin 2Tz,,
ulaa
(12)
The analytic solution Mex (t, x) = = 4x - Σ sin 2Tz, serves as ground truth. The PINN employs tanh-
MLP with layer widths 11 + 256 + 256 - 128 - 128 + 1, amounting to 118,401 parameters (n R 118k).
Each optimiser step uses Ng = 10a interior collocation points and Non = 10 boundary and initial condition collocation points, resulting in a residual dimension m 6k.
The performance of D-NGD on this problem is illustrated in Figure 1 (Left) and summarized in Table As observed, Dense D-NGD+GA achieves the lowest median relative L error of 8.52 x 10-6, closely followed by Dense D-NGD at 1.24 x 10-5, These results represent substantial improvement over the best-performing baseline, L-BFGS (9.82 x 10 10-5), by more than an order of magnitude. Compared to Adam (1.45 x 10-3) and SGD (3.48 x 10-3), the D-NGD variants are over two orders of magnitude more accurate.
5.2 Logarithmic Fokker-Planck in 9+1 dimensions
The Fokker Planck equation governs the evolution of probability densities under stochastic dynamics. We map the density to its logarithm y = logg and study:
1 1 Ao - 2 (00,.:) = lo (1)))
0gg
Vg-v
(13)
2
t E [0. 1], x 5 PINN formulations for this type of problem have been explored in Dangel et al. (2024h; ; Sun et al. (2024); Hu et al. (2024); Schwencke & Furtlehner (2024). For drift = _ and diffusion σ V21, the exact density remains Gaussian, and de x),x) log (). is available in closed form. The PINN
tanh-MI.P (10 256 256 - 128 128 1; 118,145 145 parameters). The
10

---

Published in Transactions on Machine Learning Research (10/2025)
initial condition g(0,z)) (x()) is satisfied by construction using the transformation ug(t, 2)) vo 2) the (0, + (o(2), where tig IR R IR is the core MLP output. Consequently, training focuses on the PDF residual, using Nn 3000 space-time time interior points sampled per iteration. Figure (Right) and Table summarize the performance on the Logarithmic Pokker-Plandee equation. Dense D-NGD+GA again leads with a median error of 2.47 x 10-3 followed by Dense D-NGD at 3.27 x 10-3. These D-NGD methods significantly outperform Adam (4.80 x 10-2) and SGD (5.48 10-2) by more than an order of magnitude. To the best of our knowledge, this sets new benchmark for PINNs in forward mode for this problem.
DNGD DNGD GA
Adam
L-BFGS
DNGD
Adam
L-BFGS
SGD
DNGD GA
SGD
Relative 1.2 Error
10a
Relative L2 Error
10
10
10-2
10-1
10-1
10-2
0
1000
2000
3000
1000
2000
3000
Simulation Time (s)
Simulation Time (s)
Figure 1: Left: Heat equation in 10+1 d; Right: Logarithmic
in 9+1d. We plot
the median relative L error across all runs, with shaded bands indicating the interquartile range (25th-75m percentiles) for all solvers.
5.3
Kovásznay flow at Re = 40
The steady two-dinensional Kovásznay solution is a classic benchmark for incompressible flow solvers Ko-
vasznay (1948). On the domain = [-0.5, 1.0] x [-0.5, 1.5], we solve:
(u- V Ju Vp v = O, V-u=O,
V
40
(14)
Boundary conditions are prescribed from the known analytic solution. The PINN employs tanh-MLP with four hidden layers, each with fifty neurons (7,953 parameters), similar to Jnini al. (2024). For training, 400 interior and 400 boundary collocation points sampled per iteration. Performance for the Kovásznay flow is depicted in Figure 2 (Left) and Table Dense D-NGD achieves exceptionally low median error of 5.23 x 10-7, with Dense D-NGD+GA performing similarly at 5.49 10-? L-BFGS reached 9.48 x 10 10-5, while Adam and SGD achieved errors of 4.50 x 10-3 and 6.94 x 10-3 respectively. To the best of our knowledge, this sets a new benchmark for PINNs in forward mode for this problem.
5.4 Allen-Cahn Reaction-Diffusion
The Allen-Cahn equation models phase separation and is a challenging benchmark due to its stiff reaction term and potential for sharp interface development, often causing difficulties for standard PINN training. We consider the equation on (t,x) e [0, [-1,1]] 1]:
u 10 Aura + 5
5u 0,
(15)
u(0, x) = x cos NT,
u(t, -1) = u(t, 1), u (t (t,
This involves a diffusion coefficient of 10-4, cubic reaction term, and periodic boundary conditions. The neural network is an MLP with an input layer (3 features: t, and æ after periodic embedding with period 2.0), four hidden layers (100 neurons each), and one output neuron with tanh activation leading to
11

---

Published in Transactions on Machine Learning Research (10/2025)
30,801 trainable parameters. Each training step uses N2 4,500 PDE interior points and Nan 900 boundary /initiiniiilal condition points. Results for the Allen-Cahn equation are shown in Figure 2 (Right) and Table Dense D-NGD+GA achieves the best performance with median error of 9.13 x 10 followed by Dense D-NGD at 1.21 x 10-5. This problem particularly highlights the deficiency of first-order methods and L-BFGS, which all failed to converge recover the physical solution the allocated Lime budged. Geodesic acceleration led to noticeably faster convergence for this problem.
DNGD
Adam
L-BFGS
DNGD GA
Adlam
L-BFGS
DNGD GA
SGD
DNGD
SGD
Relative 1.2 Ervor
Relative L2 Error
109
10-1
10-2
10
10-1
10
10-6
1000
2000
3000
1000
2000
3000
Simmlation Time (s)
Simulation Time (s)
Figure 2: Left: Kovásznay flow at Re = 40; Right: Allen-Cahn reaction-diffusion. We plot the
median relative L error across all runs, with shaded bands indicating the interquartile range (25th-75th percentiles) for all solvers.
5.5 Lid-driven cavity at Re 3000
High Reynolds number lid-driven cavity flow is a demanding benchmark where standard PINNs often strug- gle to achieve high accuracy beyond Re x 1000 Karniadakis et al. (2023). We consider the steady-state
incompressible Stoker equations in the domain I = [0, 1J2 with a kinematic viscosity " 3000 corresponding to Re 3000. The flow is driven by a lid moving with the profile:
u(x, 1) cosh[Co(z - +)l cosh(;CO)
1xx =0 C. 10.
(16)
No-slip boundary conditions (u = 0, = 0) are applicd on the other three stationary walls. To stabilizo training at this high Reynolds number, for our D-NGD method, we adopt the curriculum learning strategy described by |Wang et al.) (2023), involving training at progressively increasing Reynolds numbers (Re 400, and 1000) for 50 iterations each before transitioning to the target Re 3000; this approach shown to help avoid getting stuck in poor local minima. For the bascline optimizers, we employ standard curriculum approach involving 50,000 warmup iterations. The work contains about 6.6 x 10 parameters. At each step we sample 104 interior collocation points and 2 x 1 boundary points, resampling every iteration. Each run limited to 9000 s of wall-clock time. On this problem PCGD-NGD attains a median error of 3.59 x 10-3 over two orders of magnit below L-BEGS al 3.85 x 10-1 and far betler than Adam at 6.93 10-1 or SGD at 1.10. which fail to accurately reconstruct the solution field. To the best of our knowledge, this sets new benchmark for PINNs in forward mode at this Reynolds number.
5.6 Poisson Equation 10 Dimensions
We consider 10D Poisson equation, -u((x) = fx) for x E [0, 11 The source f(z) is from the analytical
solution u² (z) = Σ Ø24-1 224, so the problem becomes - -u((z) (z)) = 0 with Dirichlet boundary conditions from u (xx on 0(0,, 1). The PINN is an MLP (10 inputs; four hidden layers, 100 neurons each; output; Tanh; 411,501 parameters). Training uses m 10, 000 residual points (Ng 8. 000 interior,
Nao 2, 000 boundary). Due to the large residual dimension
10, 000, the iterative PCGD-NGD is
12

---

Published in Transactions on Machine Learning Research (10/2025)
employed. As shown in Figure (Right) and Table 1 PCGD-NGD achieves median error of 2.74 10 This is approximately twice as good as T-BFGS (5.47x10 10-1) and significantly better than Adam (3.51 x x10-2) and SGD (6.26 x 10- by about two orders of magnitude.
PCG-DNGD
SGD
PCG-DNGD
SGD
10
Adam
I-BFGS
10
Adam
L-BEGS
Relative L. Error
10"
10a
Lg Error
10
10-1
10-2
2
10-3
10-3
10-1
10-1
2000
4000
6000
1000
2000
3000
Simulation Time (s)
Simulation Time (s)
Figure 3: Left: Lid-driven cavity at Re 3000; Right: Poisson Equation in 10 Dimensions (with PCGD-NGD). We plot the median relative L error across all runs, with shaded bands indicating the interquartile range (25th-75th percentiles) for all solvers.
5.7 Poisson Equation in d = 10sh Dimensions
To gauge the large-scale behaviour of dual NGD embed the inseparable two-body test function
Mex (z)= -111|<11
Σ c sin v; + cos 1081 + I141 cos cos zi),
C ~ N(0,1.
(17)
into the unit ball B10 and pose
Au = in 1
u 0, f f== - Mpx-
(18)
A tanh-MLP with four hidden layers of 128 neurons (> 12.8M parameters) outputs d(x), and define (z)x)= (1- ||=]²²) d(x), which enforces the homogeneous Dirichlet boundary by construction.
To address the prohibitive cost of evaluating every second derivative we adopt the Stochastic Toylor Derivative Estimator (STDE) Shi et al. (2025): for each interior collocation point æ a random index set J C (1.....d) drawn and the Laplacian estimated as [J] Σ5eJ ( (Z(( 0). Each term Pu(x)(e,,0) is obtained in single Taylor-mode AD pass with the 2-jet (z, ej, 0). Training proceeds with 100 Monte- Carlo interior points per step. For this 10 Poisson problem, the performance is detailed in Figure 4 and Table [ The GA-enhanced varian tachieved a median error of 1.14 x 10-5, This result, benefiting from an approximate 3x improvement due to geodesic acceleration over Dense D-NGD (3.30 x 10 10-5), is nearly 16x better than the best-performing baseline (Adam, 1.87 10-4) and substantially lower than crrors from L-BFGS and SGD. Figure To the best of our knowledge, our work is the first to extend Natural- Gradient methods to PINNs of this scale.
13

---

Published in Transactions on Machine Learning Research (10/2025)
DNGD
Adam
L.BFGS
DNGD GA
SGD
10a
Relative L.2 Error
10a
10-2
10-4
1000
2000
3000
Simulation Time (s)
Figure 4: Poisson Equation in d 100 Dimensions We plot the median relative L error across all runs, with shaded bands indicating the interquartile range (25th-75th percentiles) for all solvers.
Quantitative comparison
Table collects the final median relative 1,2 errors for all six benchmarks.
Table 1: Median relative 2 error after the task budget (ten seeds). Best performance for each experiment is highlighted in bold. Method
Heat (10+1d)
Kov. (2d)
Poisson (10d) mR=42k
Fokker-P P.
A.-Cann (1+1d)
Cavity
Poisson
(9+ 1d)
(2d)
(10
n=118kk
n=ak
n=1188k
n=31k
n=#66k
1=133
wk=6k
m=0.8k
m= 10k
m=3k
5.66.
m==50k
100
SGD
3.48 10
6.94
10
6.26 10
5.48 10
9.93 10
1.10
3.67
Adam
1.45 10-3
4.50
10-3
10-2
4.80 x 10
10
10-1
1.87 x 10
L-BEGS
9,82 10-5
9.48
10-5
5.47
10-4
9.30 x 10
9.92 10
10-1
1.78 x 10-3
Dense D-NGD
1.24
10-5
5.23 x
3.27 x 10
1.21 x 10-5
3.30 x 10-3
Dense D NGD+GA
8.52 x 10
5.49 x 10-7
2.47 x 10
9.13 x 10
1.14 x 10
PCGD-NGD
2.74 x 10
3.59 10-5
Discussion Across all seven benchmarks (Table Figures 1 the different variants of D-NGD delivers the highest accuracy and the most reliable convergence and achieve state of the art accuracies for PINNs in several of the considered benchmarks. By swapping the intractable parameter space for an m X m solve in residual space, the method keeps the per-step cost proportional to the number of residuals rather than the number of weights. This single design choice lets us run many more curvalure-inforeed iterations within the fixed budget and, equally important, frees practitioners employ wider and deeper networks whose expressivity would otherwise be impossible exploit. Furthermore, we have shown that geodesic acceleration (GA) provides a consistent refinement negligible extra cost. By reusing the existing factorization and adding only one Hessian-voctor product per step, GA yields 25-65% lower final errors on four of the five dense benchmarks and can speed up convergence speed without degrading performance.
Gauss Newton solve in
6 Conclusion
Training high-fidelity PINNs at scale has long been hamstrung by the prohibitive cost of socond-order optimisation. By revisiting the Gauss Newton slep through a primal dual lens, Dual Natural Gradient Descent moves the heavy linear algebra into residual space, where is dramatically cheaper to assemble, store, and precondition. single Cholesky factorisation-or a handful of preconditioned CG iterations now suffices deliver curvature-ifformed updates even for networks with tens of millions parameters. The same dual operator supports geodesic-acceeeration term at negligible extra cost, further boosting step
14

---

Published in Transactions on Machine Learning Research (10/2025)
quality without hyper-paramcter tuning. Future work will explore adaptive residual sampling driven by curvature information, automated damping strategies based on stochastic spectral estimates.
Acknowledgments
A.J. acknowledges support from a fellowship provided by Leonardo S.p.A.
References
Nicholas Francis Arcolano. Approzimation of Positive Semidefinite Matrices Using the Nyström Method. PhD thesis, Harvard University, 2011. UMI Number: 3462777.
Frederik Benzing. Gradient descent on neurons and its link to approximate second-order optimization, 2022.
URL https://aaav..org/ab///
12250
Jesse Bettencourt, Matthew J. Johnson, and David Duvenaud. Taylor-mode automatic differentiation for higher-order derivatives in JAX. In Program Transformations for ML Workshop al NewrIPS 2019, 2019.
URL https://openreey: http:/rrwwweeaeeeeeeoooorrroaaas wet.//orruriidissxssnnnn
Mathieu Blondel, Quentin Berthet, Marco Cuturi, Roy Frostig, Stephan Hoyer, Felipe Llinarcs-T6pez Fabian Pedregosa, and Jean-Philippe Vert. Efficient and modular implicit differentiation. arXiv preprint 1.Xi::2105.1118,, 2021.
Andrea Bonfanti, Giuseppe Bruno, and Cristina Cipriani. The challenges of the nonlinear regime for physics- informed neural networks, 2024. URL http://////riioooraaba/////s//// 03864
James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adamn Paszke, Jake VanderPlas, Skye Wanderman-Milnee, and Qiao Zhang. JAX: compos- able transformations of Python+NumPP programs, 2018. URL. http:://gthhhbb
Shengze Cai, Zhiping Mao, Zhicheng Wang, Minglang Yin, and George Em Karniadakis. Physics-informed neural networks (PINNs) for fluid mechanics: review. Acta Mechanica Sinica, 37(12):1727-1738, 2021.
Felix Dangel, Johannes Müller, and Marius Zeinhofer. Kronecker-factoree apppo:/ curvature for physics- informed neural networks, 2024. URL hntps///aaxxv..og//a///200....
M.W.M.G. Dissanayake and N. Phan-Thien. Neural-nettwork-based approximations for solving partial dif- ferential equations. Communications in Numerical Methods in Engineering, 10(3):195 201, 1994.
Nikita Doikov and Anton Rodomanov. Polynomial preconditioning for gradient methods. In 40th Interna- lional Conference on Machine Learning (ICML), number 40, 2023.
Nikita Doikov, Sebastian U. Stich, and Martin Jaggi. Spectral preconditioning for gradient methods on graded non-convex functions, 2024. URI. htttss://axxii http:/wwwwrggbb///4...ooo
Zachary Frangella, Joel A. Tropp, and Madeleine Udell. Randomized nyström preconditioning, SIAM
Journal on Matrix Analysis and Applications, 41(3):118111212, 2023. doi: 10. 111113721MM1466244
Ehsan Haghighat, Maziar Raissi, Adrian Moure, Hector Gomez, and Ruben Juanes. A physics-informed deep learning framework for inversion and surrogate modeling in solid mechanics. Computer Methods in Applied Mechanics and Engineering, 379:113741, 2021.
H. V. Henderson and R. Searle. On deriving the inverse of a sum of matrices. SIAM Revicu, 23(1):53 60, 1981.
Magnus R. Hestenes and Eduard Stiefel. Methods of conjugate gradients for solving linear systems. Journal of Research of the National Bureau of Standards, 49(6):409-453 436, 1952.
15

---

Published in Transactions on Machine Learning Research (10/2025)
Zheyuan Hu, Khemraj Shukla, George Em Karniadakis, and Kenji Kawaguchi. Tackling the curse of dimen- sionality with physics-informed neural networks. orXiv preprinl arXiv:2307.12 1000007.13000 12306, 2023.
Zheyuan Hu, Zhongqiang Zhang, George Em Karniadakis, and Kenji Kawaguchi. Score-based physics-
informed neural networks for high-dimensional fokker-planck equations, 2024. URL https://arxiw.orrg .org
abs/2402. 07465
Anas Jnini and Flavio Vella. Taylor mode neural operators: Enhancing computational efficiency in physics- informed neural operators.
Anas Jnini, Flavio Vella, and Marius Zeinhofer. Gauss-newton natural gradient descent for physics-informed computational fluid dynamics, 2024. URL httpss//aaxxw..oog///////000 20000000020020000000
Anas Jnini, Lorenzo Breschi, and Flavio Vella. Riemann tensor neural networks: Learning conservative sys- tems with physics-consttrained networks. In Forty-sccond International Conference on Machine Learning, 2025a. URL htttps://openrevieet http::/aww....eeeo//////////ooaa
Anas Jnini, Harshinee Goordoyal, Sujal Dave, Flavio Vclla, Katharine H Fraser, and Artem Korobenko. Phystee-constraanee deeponet for surrogate cfd models: curved backward-facine step case. orXiv preprint arXiv:2503.11 11 1169, 2025b.
George Karniadakis, Zhicheng Wang, and Xuhui Meng. Solution multiplicity and effects of data and eddy viscosity navier-stokes solutions inferred by physics-informed neural networks. Bulletin of the American Physical Society, 2023.
George Em Karniadakis, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. Physics- informed machine learning. Nature Reviews Physics, 3(6):422-44D) 2021.
Diederík P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. 2017.
Elham Kiyani, Khemraj Shukla, Jorge F. Urbán, Jérôme Darbon, and George Em Karniadakis. Which optimizer works best for physics-informed neural networks and kolmogorov-aroold networks?, 2025. URL
http:://aaxxi..og/ba//////////// 1637
L. G. Kovasznay. Laminar flow behind a two-dimensional grid. Mathematical Proceedings of the Cambridge Philosophical Society, 44(1):58-62, 1948. doi: 10. 10:1017/5030504400033399
K. Levenberg. method for the solution of certain non-lincar problems in least squares. Quarterly of Applied Mathematics, 2(2):164 168, 1944.
Zongyi Li, Hongkai Zheng, Nikola Kovachki, David Jin, Haoxuan Chen, Boren Liu, Kamyar Avizzaenesheeli, and Anima Anandkumar. Physics-informer neural operator for learning partial differential equations. arXiv preprind xxX:::2111.08794, 2021.
D. W. Marquardt. An algorithm for least-squares estimation of nonlinear parameters. SIAM Journal on Applied Mathematics, 11(2):431-444 1963.
James Martens. Deep learning via hessian-free optimization. pp. 735 742, 08 2010.
Per-Gumnar Martinsson and Joel A. Tropp. Randomized numerical linear algebra: Foundations and algo- rithms. Acta Numerica, 29:403-572, 372, 2020. doi: 10.1017/500622492200000
Maricela Best Mckay, Avleen Kaur, Chen Greif, and Brian Wetton. Near-optimal sketchy natural gradients for physics-informed neural networks. In Forty- International Conference on Machine Learning, 2025. URL htttps://peenrviies feret/fooummidee http//rrrmmmdd-bkkkommmmmnn
Johannes Müller and Marius Zeinhofer. Achieving high accuracy with PINNs via energy natural gradient descent. In Andreas Krause, Emma Brunskill, Kyunghyun Cho, Barbara Engelhardt, Sivan Sabato, and Jonathan Scarlett (ods.), Proceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research, pp. 25471-25485. PMLR, 23-29 Jul 2023. URL mlr. press//202/mullr223 html
16

---

Published in Transactions on Machine Learning Research (10/2025)
Johannes Müller and Marius Zeinhofer. Position: Optimization in sciml should employ the function space geometry, 2024. URL
07318
Maziar Raissi, Paris Perdikaris, and George Karniadakis. Physics-informed neural networks: deep learn- ing framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics, 378:686-7017 2019.
Anton Rodomanov and Dmitry Kropotov. Volume sampling in optimization. SIAM Journal on Optimization, 30(3):1878 1904, 2020.
Tim De Ryck, Florent Bonnet, Siddhartha Mishra, and Emmanuel de Bézenac. An operator preconditioning perspective on training in physics-informed informed machine learning, 2024. URL httph: .rg. arg/abbs/310. 05801
Yousef Saad. Iterative Methods for Sparse Linear Systems. STAM, 2 edition, 2003.
Nicol N Schraudolph. Fast curvature matrix-vector products for second-order gradient descent. Newral computation, 14(7):1723 1738,,177,81880 2002.
Nilo Schwencke and Cyril Furtlehner. Anagram: natural gradient relative to adapted model for efficient pinns learning, Xiv preprint or arXiv:2412.1076.. 2024.
Jonathan Richard Shewchuk. An introduction the conjugate gradient method without the agonizing pain, 1994. Technical Report, Carnegie Mellon University.
Zekun Shi, Zheyuan Hu, Min Lin, and Kenji Kawaguchi. Stochastic taylor derivative estimator: Efficient amortization for arbitrary differential operators, 2025. URL httpps//aaaxxv..rr////////2221
Leslie N. Smith and Nicholay Topin. Super-coonvergence:: Very fast training of neural networks using large learning rates, 2018. URL httpp::/aaaxvv.orr//b////// 07120
Yang Song, Jiaming Song, and Stefano Ermon. Accelerating natural gradient with higher-order invariance, 2018. URL
01273
Jingtong Sun, Julius Berner, Lorenz Richter, Marius Zeinhofer, Johannes Müller, Kamyar Avizzadenesheli, and Anima Anandkumar. Dynamical measure transport and neural pde solvers for sampling, 2024. URL
/ttps://arzile www.aaiwwwww.aaaaaaeee borg/abs/200 2407 2t/.//4077007777 07873
Sifan Wang, Yujun Teng, and Paris Perdikaris. Understanding and mitigating gradient pathologies in physics- informed neural networks, 2020. URT https::/aarxww.oor//bb/////////s 04536
Sifan Wang, Han Wang, and Paris Perdikaris. Learning the solution operator of parametric partial differential equations with daysics-informer deeponets. SIAM Journal on Scientific Computing, (335):33055 A3081, 2021.
Sifan Wang, Shyam Sankaran, Hanwen Wang, and Paris Perdikaris. An expert's guide to training physics- informed neural networks, 2023.
Sifan Wang, Ananyae Kumar Bhartari, Bowen Li, and Paris Perdikaris. Gradient alignment in physics- informed neural networks: A second-order optimization perspective, 2025. URL https://arxiv.org org/ abs/2502.00604
Stefano Zampini, Umberto Zerbinati, George Turkiyyah, and David Keyes. Petseml: Second-order solvers for
training regression problems in scientific machine learning, 2024. URL https://arxdv.org/ abs//bb3.
12188
Qi Zeng, Spencer H. Bryngelson, and Florian Tobias Schaefer. Competitive physics informed networks. ICLR 2022 Workshop on Gamification and Multiagent Solutions, 2022. URL /tttt::////// http:///pppeeeeeeree scJ61c
17

---

Published in Transactions on Machine Learning Research (10/2025)
A Appendix
A.1 Proofs of Propositions and Theorems
A.1.1 Proof of Proposition 3.1 (Primal normal equations)
The objective function is
PA(A0)) = (+(()))+++))))) + J(0,) (r())()(())))))))))))))))))))) + (Op,) A02) + NN A0/ A0n-
To find the minimizer A we compute the gradient of Fx (A(r)) with respect to Ad, and set to zero.
Let's expand the first term:
((r(0))+ A0.) (aax:(((r(()))))))))))))( (r(0x) + ((o,) A0&) = 2 2 + (0,) J(OR) A0%)
22(2)) (r((0)) +r.0)) J(Ok) A0k + J(0,) J(Ok)) AOK-
The gradient of this term with respect to A0g is:
Vas. (fr(Ow) r(O.) + rOw)" ((0z) A0 JOz) J(Os) A0&) J(0w) + J(0x)" J(Ox) J(Ok) A0x- The gradient of the regularization term
AO A0 with respect to A0 is:
Ve ( A0] = A
Combining these, the gradient of (x)) is:
Vao, Fx(40)) J(0k) (r()))r(++)))))))))))))) J(Ox)) J(0b) A + +AO
Setting the gradient zero for the optimal AO;:
J(Ox) + [(O))] J(0s) AØ AI, A6%- = 0.
Rearranging the terms, get:
(J(B)) J(Ox) + AI,) A = ((()))r((())
Given that Vol(Ok) is defined as J(O)) ((Ox))((O)) in this context (representing the gradient of the unweighted least-squares loss 1 (0x)),)))),, we have:
(J(O2) J(0e) An) A VoL((0.)
The matrix J(Ox)T J(Ox) is positive semi-definite For A > O, the matrix J(Ox)T J(Bb) + AI, is positive definite, ensuring that (A(0)) is strietly convex and thus has a unique minimizer. This completes the proof of Proposition 3.1
A.1.2 Proof of Theorem 3.1 (Dual Normal Equations)
The Lagrangian is given by:
L(A0,, Uk. uk) = (11)(())))))))))))))))))))))([ + 11 + 1100 + Us 1yx- - J(O2) A0x)
Setting the partial derivatives of L(AOk, Jk, yk., uk) with respect to Adk, yks and us to zero yields the KKT system:
Aa, VA. 8: XA0, J(0k) 1/ 0,
(19)
V L: r(0x) Ux 1 = 0,
(20)
L: y4 ((x)) 20; 0.
(21)
18

---

Published in Transactions on Machine Learning Research (10/2025)
From equation 20 Vk -r(0k) y/-- Substituting into equation 19 gives
A A -(0%) J(Ok) (r(0x) Mk).
Using (gl((0)) J(0x)] Tr(0x) (gradient of the unweighted least-squares loss (11r(0+)1B), this becomes À AB& - VoI((k) - =[0x) yk, which rearranges to Equation equation 7
A0 - + VoL(Oe))
Substituting this expression for A0t into equation
21 yields for the optimal yt=
v ((k) A - ((Aa) (J(OL) vi + voL(Ox))
Multiplying by X and expanding gives:
xy ((())))((O) vt J(OL)VoL(O)))
Using the definition K₁ = (Ok) (Definition 3.1 we have:
Ayt KkvV ((x)) VoL(
which Equation equation 6 The equivalence between the primal system equation 3.1 and this dual formulation be confirmed by substitution. Becauss the objective is strictly convex and the equality constraints are affine with non-empty feasible set, the KKT conditions are necessary and sufficient and strong duality holds. This completes the proof of Theorem 3.1
Rearranging gives (K& + =
A.1.3 Proof of Proposition 3.2 (Primal and Dual GA Characterization)
First, the primal condition follows by differentiating (Ox)a + fvull + á1aa2 w.r.t.
(11J((0)) 11 }llal²) = J(0k)
[(0x))( (J(Os)a + foo) + Xa,
and setting this zero yields equation
(J(b)) J(O+) AI)) a = J(bb fee.
To derive the dual formulation, we introduce the auxiliary variable y J(0g) a and enforce the constraint y J(x)aa via Lagrange multiplier m The Lagrangian for the minimization problem in equation S is
L(a, y. v) 11 foell2+ 2| + Ty J(O))a).
The KKT conditions are obtained by setting the partial derivatives to zero:
VaC: Xa - J(0k) Tu = 0,
VyC: : y fvv + V 0,
V.L: y J(Ox)a 0.
From VyC = O, we get v = = - + fwo). Substituting this into VaL = gives Xa = -((x)) Toy))(y + fee).
which, for the optimal as and corresponding ya matches the expression for ag in equation 10
a = - ((x) (ya,k + fee).
Meanwhile, from V.L = 0, we have y = J(x))a Substituting the expression for a:
y = )(w)) (-)0)(=+ J(O2) (y + (oo) = ((J(O))((0)))) (y foo).
Using Ks (Oxa) this becomes
== - KR (y (vw)-
Rearranging gives Xy - Kafees so (K₂- + Amm) y= Kr fov. For the optimal Yo.k: this the first
part of equation TO This completes the proof of Proposition
3.2
19

---

Published in Transactions on Machine Learning Research (10/2025)
A.2 Algorithmic Implementations
This section provides the detailed algorithms referenced in the main text.
Algorithm KernelEntry via double VJP
1: Input: collocation points 24, Zj. parameters 0, residual fn. Tia
2: Define fi(0) = Tiu (2, 0), fy( = (z)
(Mis vjpy) + jax. vip(sj,0) 0)
(ys, vip;) -- jax.vip(fs.0
Initialize B Qixd {d dim of each residual}
for k == to d do
7: +- (ex) (backprop seed e through fi}
8: e (u) {backprop u through fi)
9: B + fcolumn of B}
10: end for
11: return B JJ}
A.2.1 Algorithm: Assemble Residual Gramian IC
Algorithm Assemble residual Gramian
1: Input: collocation points (r,)m parameters 0, residual fn. rin
2: mi + Nqda, 122 t- Nandaq (lype split]
3: x 0mxm
for i
m do for + to m do
if (i m Aj (i> A>>mAj>N MA <) then continue
8:
end if
9:
B 4 KERNELENNTRY(E,,233 rtm)
10:
Kij BB
11:
if fitjon or (is n (ssnAj>m) or (1 ni Aj S n1) then
12:
Ky K B
13: end if
14: end for
15: end for
16: return K.
20

---

Published in Transactions on Machine Learning Research (10/2025)
A.2.2 Algorithm:Dense Dual Solve
Algorithm Dense dual solver step (DINS:DUALSSOIVE)
1: Input: 8, VoL(O), rin, collocation sets (Xa, Xan), damping X, flag use_geodesic acceleration. oot
2: K + ASSPMBLeGGAAMIN points, 0, rm.
K+ K + Amm
4: Lehol + chol(K)
5: b+ -(() 90
Solve Lchol y/ = then Le = y
v - (JOO] y² + go)
use Bece__egodesc___eeee acceleration then
9. fvo - dra Mtn (poiniiss + tv) t=
&
10: ba + -( fvv
11: Solve Lehoi Va = ba then Lehotila = Yo
12: -- - (J(OT (y + fwo))
13:
AO +
14:
if ||ul|2 Cuorm then
15:
< 2 ||aa|2/||ul|]
16:
if r S 0.5 then
17:
AØ < A6 a.5
18: end if
19: end
20: else
21: AO
22: end if
23: return A0
A.3 Algoritlm:Preconddtioneed Conjugate Gradient
Algorithm Nyström construction of U,A
1: Input: 0, residual pts. {I,) Pin, landmarks
2: Choose landmark indox sct I, [I = E; let C = (1,...,,m) m3 /
Form Krr and Ker via KERNELENNRR
Eigendecompose K11 QA fr positive eigenpairs)
U, Q; UC KoIQA="
SVD UA1/2 VEET
U+ V; Å A+D E 8: return (U
21

---

Published in Transactions on Machine Learning Research (10/2025)
Algorithm 6 Pre-conditioned CG step (PCGSTEP)
1: Input: 0, g V ((), Tin, { X, rank f, tol. &, Mmax (U, + BOULD.VYSTTOAAPPPOXXAAAA Build preconditioner p-1 using (U. A)
(Ee),Tin 2)
b+ -J(0)g y 0 r<b
= p- p p (T)
6: for j + 0 to mmax - do
7: Ap - + Ap {two JVP/W VJP calls}
8: if (p, Ap) = then
break
10: end if
11:
a p/ Ap)
12:
y y ap r aAp
13:
if ||rll S Selbbl22 then
14:
break
15:
end
16: Z + p-tr; Pnew + (r,
17: if p 0 then
18: break
19: end if
20: f- Pnew/p: /p: p 4- 3p; x- Pnew
21: end for
22: AO + -----1110 ++>-~1(J(0)) ++
23: return A6
A.3.1 Algorithm: Dual Natural-Gradien Descent Workflow
A.4 Hyperparameters for Allen-Cahn Reaction-Diffusion Experiment
Table 2: + 1-D Allen Cahn, (t,z) E 0, x [1,1]] diffusion 10-4 Setting
Category
PDE
44 10 Unn + 5u 5u = 0; periodic BCs; u00, = x² (TT)
Network
MLP, tanh; layers [3, 100, 100, 100, 100, 1]; R 801 params
Training
Nn 500, Nan 900; budget 000 s; 10 seeds
Dense D-NGD
Dense Cholesky; d = min (loss, 10 5) 31-pt line search; with /without GA
Adam
n = 10 10-3; 3; (81, (2) (0.9, 0.999);
10-3
SGD
One-cycle (peak 10-3 final 10 10-4); momentum 0.9
L-BFGS
Jaxopt; history 300; strong- Wolfe; tol 10 6
Wolle:
A.5 Hyperparameters for 10+1-Dimensional Heat Equation Experiment
Table 3: Heat on [0. 1| [0. 1), == 1 MLP, tanh; [11, 256, 256, 128, 128, 1]; 118 401 params
Network
Training
No = 10 000, Nog 1 000; budget 3 000 s; 10 seeds
Dense DNGD
À min(loss, 10 10-3); 31-pt line search; GA variant identical
Baselines
Adam / SGD / L-BFGS as in Table 2
22

---

Published in Transactions on Machine Learning Research (10/2025)
A.6 Hyperparameters for Logarithmic Fokker-Planck Experiment
Table 4: Eq. equation 13
E [-5, [0,
Network
MLP, tanh; [10, 256, 256, 128, 128, 1 118 145 params
Training
Interior residuals only, No 3 000; budget 3 000 s; 10 seeds
Dense DNGD
Àk min(loss, 10 9); 31-pt line search; GA variant identical
Baselines
Adam SGD / L- LFGS as in Table 2
A.7
Hyperparameters for Kovásznay Flow Experiment
Table 5: Steady 2-D Kovásznay benchmark (Re 40).
Network
MLP, tanh; [2, 50, 50, 50, 50, 3): 7953 params
Training
No = 400, Nao = 400; budget 3 000 s; 10 seeds
À min(loss, 10 10-5); 31-pt line scarch; GA variant identical
Baselines
Adam SGD / L-BFGS as in Table 2
A.8
Hyperparameters for Lid-Driven Cavity Experiment
Table 6: Steady lid-driven cavity (Re 3000; budget 9 000 s)-
Network
MLP, tanh; [2, 128, 128, 128, 128, 128, 3): = 66 000 params
Training
Nn = 10000, Nao 2 000; 10 seeds
Curriculum (DNGD)
Re 100, 400, 1000 (50 it. each) Re 3000
Curriculum (baselines)
50 000 warm-up iterations at each Re 100, 400, 1000 (per Wang et al. (2023
Re 3000
Iterative PCGD-NGD
Nyström rank 2500; CG tol 10-10; max 500; A min 10 10-9); 31-pt line search
Baselines
Adam (exp-decay LR), SGD (one-cycle), L-BFGS (history 300)
A.9 Hyperparameters for 10-Dimensional Poisson Experiment
Table 7: Laplace on [0, 1] with analytic Dirichlet BCs.
Network
MLP, tanh; [10, 100, 100, 100, 100, 1]: 41 501 params
Training
No = 8 000, Nog = 2 000; budget 3 000 s; 10 seeds
Iterative PCGD-NGD
Nyström rank 2500; CG tol 10
10-1; 10; max 500; A = min 10 10-5); 31-pt line
search; no GA
Baselines
Adam SGD
L-BFGS as in Table 2
A.10
Hyperparameters for 103-Dimensional Poisson Experiment
Table 8: Poisson on I0
Network
MLP, tanh; 1100000, 128, 128, 128, 128, 1]; R 12.8M params
Training
Nn = 100 (STDE, re-sample each step); budget 000 s; 10 seeds
Dense DNGD
Dense Cholesky: 2 = min(loss, 10*); 5-pt line search
Dense DNGD +GA
Same with geodesic acceleration
Baselines
Adam SGD L-BFGS as in Table
23

---

Published in Transactions on Machine Learning Research (10/2025)
A.11 Navier-Stokes (m n) Sweep: Primal GN Dual D-NGD (500 iters)
We report the wall-time required to complete 500 optimizer iterations for the steady Sevier-Stokes setup across a grid of residual counts and parameter counts Two solvers are compared: (i) Primal Gauss- Newton/Levenberg
in paramcter space; (ii) Dual D-NGD, which computes the same step in residual space. All experiments in this section were run on an NVIDIA A100 GPU.
Experimental note. Each cell in the following two tables shows the average time per iteration in seconds, calculated from total run of 500 iterations for the given n).
Table 9: Primal (parameter-space) Gauss Newton: time per iteration (seconds) across (m,
m
n 303
n 2853
n 5403
n 7953
n 10503
n 13053
n 15603
n 18153
m
= 300
0.0012
0.0054
0.0146
0.0294
0.0480
0.0824
0.1226
0.1681
1800
0.0039
0.0210
0.0419
0.0676
0.1062
0.1554
0.2138
0.2868
3300
0.0045
0.0556
0.0969
0.1496
0.2121
0.2898
0.3857
0.4917
ra
4800
0.0060
0.1031
0.1805
0,2705
0.3734
0.4913
0.6278
0.7791
6300
0.0084
0.1732
0.2985
0.4389
0.5950
0.7675
0.9606
Table 10:
Dual D-NGD: time per iteration (seconds).
303
2853
n 5403
n 7953
т 10503
13053
15603
18153
300
0.0015
0.0025
0.0035
0.0043
0.0053
0.0061
0.0069
0.0078
T7l
1800
0.0059
0.0072
0.0098
0.0127
0.0138
0.0164
0.0188
0.0213
3300
0.0069
0.0120
0.0179
0.0232
0.0266
0.0307
0.0325
0.0418
m = 4800
0.0099
0.0194
0.0290
0.0344
0.0440
0.0466
0.0511
0.0657
6300
0.0161
0.0308
0.0416
0.0591
0.0639
0.0726
0.0807
Table 11: Primal vs. Dual regime map (empirical): winner at each (m, n) based on measured time per
iteration. Green Dual faster; Orange Primal faster.
303
2853
5403
7953
10503 n 13053
15603
18153
300
Primal
Dual
Dual
Dual
Dual
Dual
Dual
Dual
m = 1800
Primal
Dual
Dual
Dual
Dual
Dual
Dual
Dual
m 3300
Primal
Dual
Dual
Dual
Dual
Dual
Dual
Dual
m 4800
Primal
Dual
Dual
Dual
Dual
Dual
Dual
Dual
m 6300
Primal
Dual
Dual
Dual
Dual
Dual
Dual
A.12 Effect of the Number of Landmarks on the Nyström Preconditioner
The efficiency of the Preconditioned Conjugate Gradient (PCG) solver hinges on the quality of the precon- ditioner. Our choice of Nyström-based spectral proconditioner motivated by the empirical observation that the PINN Gramian matrix, K,, often exhibits a rapidly decaying eigenspectrum. This property implies that low-rank approximation can effectively capture the dominant spectral information responsible for il-conditioning.
To substantiate this, we conducted an analysis on the Lid-driven cavily benchmark at Re=3000, where the number of residuals is large (m = 50, 000). Figure 5 plots the eigenvalues the Gramian matrix at different stages of training, clearly illustrating the rapid spectral decay where a small fraction of eigenvalues contains most of the spectral energy. The analysis in Table 12 further demonstrates the practical benefit of this property, showing sharp decrease in the number of PCG iterations required for convergence as the
24

---

Published in Transactions on Machine Learning Research (10/2025)
number of landmarks (1) in the Nyström approximation increases. A relatively small < m is sufficient to significantly accelerate the solver, confirming the effectiveness and efficiency of our preconditioning strategy.
Eigenealuc Spoctrum Evolution (First 3000 Eigeavalues)
100
107
Magninal
03"4
1008
2006
20.00
(090
3000
1000
2000
1000
1000
5000
Zigtarallr Luiles
Esgausalo:: Indlo
[tecation 100
Etceation 250
([o soslo)
10n-
00-1
Maguinda
22
2000
[tecation
[teration 430
150
20-1
33
Maquicale
201
1006
2066
2000
4003
930
1000
2000
12:0
4000
Bagenooliie Lndex
Dipersloo lncke
Figure 5: Evolution of the first 5000 eigenvalues of the Gramian matrix at different training iterations for the Lid-driven cavity problem. The plots show that the spectrum decays rapidly in the earlier stages of training.
Table 12: Convergence of the Preconditioned Conjugate Gradient (PCG) solver a function of training iteration and the number of landmarks (I) used in the Nyström approximation. Each cell shows the number of iterations required to reach the convergence tolerance. A value of 500 indicates that the solver did not converge within the maximum allowed iterations, highlighting the necessity of a sufficiently large rank for the preconditioner, particularly in later, more challenging stages of training.
Epoch
CG Iterations for Number Landmarks (1)
500
1000
2500
4000
5000
10
10
10
10
10
50
500
172
11
10
10
100
500
337
20
11
11
250
500
500
230
38
17
350
500
500
338
60
23
450
500
500
457
80
31
25


Reviews:
[{"review_id": "VG2wydUbzF", "reviewer": "Reviewer_ttkY", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths\n\n- The manuscript is clearly written, and the proposed advantage of using the dual problem when the batch-output dimension is much smaller than the NN parameter dimension appears sound.\n- The benchmark problems are quite comprehensive, highlighting an overall improvement of the method over standard gradient descent methods.\n\nWeaknesses\n\n- A comparison with a recent work claiming similar level of relative accuracy would be illuminating: Elham Kiyani, Khemraj Shukla, Jorge F. Urbán, Jérôme Darbon, George Em Karniadakis, Which Optimizer Works Best for Physics-Informed Neural Networks and Kolmogorov-Arnold Networks? (https://arxiv.org/abs/2501.16371). They propose a method called self-similar BFGS and Broyden methods. Given the good performance of the methods in the manuscript, comparison showing similar performance would bolster the results in the paper.\n\n- The authors should discuss why they assume the number of collocation points (or the output dimension $m$) should be much less than the number of parameters. Although this is a familiar setting in generic deep learning settings, in the case of PINNs, I am uncertain if the number of parameters should be much more than the number of collocation points.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper proposes a gradient descent method for minimizing physics-informed losses, called the dual natural gradient descent. The method builds on a previously introduced Gauss-Newton Natural Gradient (GNNG), but proposes to solve a dual problem. The dual problem's Gramian has matrix size that scales with the number of collocation points (spatiotemporal collocation points) rather than the number of parameters in the neural network model. Assuming the former is significantly smaller than the latter, the cost of computing the dual problem is more efficient. The new method is named D-NGD and is able to achieve L2 relative error below 1e-4 for benchmark problems.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- A performance comparison with self-similar BFGS and Broyden methods for the set of experiments in the manuscript would be illuminating.\n\n- A discussion regarding $n$ vs. $m$ and when the scaling of this method becomes favorable in terms of the wall-time. E.g. repeat of the experiments with different values of $n$ and $m$ would be helpful.", "broader_impact_concerns": "None.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "cOaCKOlwho", "reviewer": "Reviewer_U1wM", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The submission demonstrates several strong aspects, particularly in its mathematical rigor and clarity. It presents well-defined algorithms and proofs, showcasing a robust theoretical foundation. The introduction of the Dual Natural Gradient PINN is noteworthy, supported by extensive experimental results across various scenarios that effectively illustrate the generalizability of the research. Additionally, the paper integrates a diverse range of ideas into the context of PINNs, paving the way for potential future research directions.\n\nHowever, there are some weaknesses that require attention. While the paper emphasizes the advantages of DNGD as an improved version of Gauss-Newton, it primarily compares this method with first-order techniques like Adam and SGD. A comparison with other second-order methods is essential, especially given the claim of superior time complexity. It would be beneficial to evaluate how quickly and effectively DNGD converges to a low Relative L2 error relative to these alternatives.\n\nFurthermore, in the Iterative Dual Solver process, the mention of the number of landmarks for the Nyström approximation is vague, only stating that it is less than a certain value. A more detailed explanation of what constitutes an appropriate value for this parameter is needed, along with specifics on its usage in the experiments. This applies to the variable in the Low-Rank Nyström Spectral Preconditioner for the Dual System as well. Without clear information on these values, it is challenging to support claims of effective approximation.\n\nLastly, the variable Cdense_thresh in the algorithm section lacks context. Given the paper's focus on reducing time complexity, criteria for this variable are necessary to assess whether the claimed advantages in time complexity are valid. Providing such details would strengthen the overall submission and enhance its credibility.\n\n\nHere is the additional questions:\n\n### Questions\n\n- How does D-NGD compare to the Gauss-Newton method in terms of results?\n\n- While I believe that time complexity alone is sufficient since it is higher than memory complexity, could the paper also address the advantage in memory complexity as a metric?\n\n- The degree of approximation of \\ell in the Low-Rank Nyström Spectral Preconditioner for the Dual System likely varies significantly depending on the value of the variable. Could you explain how \\ell is determined or how results vary based on this value?\n\n- The algorithm applies different methods based on the Cdense_thresh criterion, but the experimental results do not distinguish between them. It would be beneficial to present experimental results for both D-NGD-DDS and D-NGD-PCCS, categorized by Cdense_thresh, to justify this distinction.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The submission addresses the challenges faced in solving Physics-Informed Neural Networks (PINNs) by comparing various optimization methods. It highlights that traditional approaches, such as Adam and SGD, do not yield satisfactory accuracy, making Gauss-Newton a more favorable choice due to its ability to utilize geometric information from the functional space. However, Gauss-Newton's time and memory complexities remain significant, particularly concerning the parameter and residual space sizes.\n\nTo overcome these limitations, the submission introduces the Dual Formulation via KKT in the residual space, which allows for a reduced time complexity compared to the conventional Gauss-Newton method. Additionally, it integrates geodesic acceleration to enhance accuracy further. Recognizing the computational difficulties associated with the Hessian matrix, the authors propose a Hessian-free solution for the dual system, employing the Conjugate Gradient method for efficient approximation. This innovative approach provides a pathway to improve the performance of PINNs while managing computational resources effectively.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "See the above questions and weaknesses.", "broader_impact_concerns": "The work presents no significant ethical concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "Y0ChZnATX8", "reviewer": "Reviewer_Rjch", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\n1. The proposed approaches are theoretically principled,\n2. Empirical evidence demonstrates improvements upon the chosen baselines while taking the uncertainties into consideration.\n\nWeaknesses:\n\n1. Some potential related works are not discussed. See Requested Changes.\n2. Relatedly, Song et al. 2018. noted that using the geodesic acceleration is not an essential component of NGD itself, but rather an approximation of the Riemannian Euler method. It might be beneficial to discuss this perspective further in the paper.\n\nSong et al., Accelerating Natural Gradient with Higher-Order Invariance, ICML 2018", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "Physics-Informed Neural Networks (PINNs) are proposed as a neural approach to solving partial differential equations more efficiently. However, they often suffer from training instabilities, which has been hypothesized to be due to difficult optimizations. Natural Gradient Descent (NGD) has been proposed to improve optimizations, yet previous attempts often do not scale to large models. The paper proposes to perform NGD in a dual space given by the residuals, thus improving the scalabilities. Using geodesic acceleration as a correction term was shown to improve the performances. Nyström spectral preconditioner helps to improve the computation speed for large batch sizes. The benefits of the proposed approach are demonstrated on several PINN benchmarks, showing improvements upon the baselines.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "I am not an expert in PINN, but I am reasonably familiar with natural gradient descent. I find the explorations on improving the training of PINNs using NGD interesting, while not being able to determine their significance due to lack of familarities. However, it would be beneficial if the authors can discuss the relation of the proposed approaches to the following works:\n\n1. While not its core contribution, similar idea on efficiently solving for natural gradient update was presented in Appendix I of Benzing 2022, where it was demonstrated that the exact NGD update can be obtained in a lower dimensional space by invoking the matrix inversion lemma. This appears related to the proposed primal-dual view.\n\n2. Song et al. 2018 proposed to accelerate NGD by accounting for the  geodesic acceleration. This appears rather related to the geodesic acceleration method proposed in the paper.\n\nBenzing, Gradient Descent on Neurons and its Link to Approximate Second-Order Optimization, ICML 2022\n\nSong et al., Accelerating Natural Gradient with Higher-Order Invariance, ICML 2018", "broader_impact_concerns": "The work is largely theoretical, and I do not find major concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)