Energy Manifold Natural Gradient Descent: Riemannian Optimization for Neural PDE Solvers

# Energy Manifold Natural Gradient Descent: Riemannian Optimization for Neural PDE Solvers
Source: [https://arxiv.org/html/2607.22004](https://arxiv.org/html/2607.22004)
\\nameZhangyong Liang\\emailzyliang1994@tju\.edu\.cn \\addrNational Center for Applied Mathematics Tianjin University Tianjin, 300072, China\\nameHuanhuan Gao\\emailgao\_huanhuan@jlu\.edu\.cn \\addrSchool of Mechanical and Aerospace Engineering Jilin University Changchun, 130025, China

###### Abstract

Energy natural gradient descent \(ENGD\) aligns parameter updates with the curvature of an underlying function\-space energy, but existing formulations assume an unconstrained Euclidean parameter domain\. We introduceEnergyManifoldNaturalGradientDescent \(EMNGD\), a manifold optimization framework for physics\-informed and variational neural PDE solvers whose parameters lie on a Riemannian manifold\. EMNGD restricts the energy\-induced quadratic model to feasible tangent directions and uses retractions to preserve parameter constraints throughout optimization\. Under coercivity, we prove that the push\-forward of the undamped EMNGD direction is the best feasible approximation to the function\-space Newton vector in the energy metric\. We establish coordinate invariance, exact reduction to ENGD in Euclidean space, global first\-order convergence with Armijo backtracking, and robustness to inexact tangent solves\. For quadratic residual energies and generalized Gauss–Newton pullbacks, the Woodbury identity transfers the tangent system to sample space without changing the direction\. Nyström approximation provides scalable sample\-space solves with controlled direction error and recovers the exact direction after iterative convergence\. On the evaluated neural PDE benchmarks, EMNGD achieves higher accuracy and faster convergence than the compared state\-of\-the\-art baselines\. Woodbury preserves the EMNGD direction, while scalable\-solver diagnostics quantify the accuracy–cost trade\-off of preconditioning and residual subsampling\.

Keywords:energy natural gradient descent, manifold optimization, neural PDE solvers, woodbury identity, nyström approximation

## 1Introduction

Neural PDE solvers parameterize the unknown solution and minimize a PDE\-based loss\. Early work used residual minimization with neural networks\(Dissanayake and Phan\-Thien,[1994](https://arxiv.org/html/2607.22004#bib.bib243); Lagariset al\.,[1998](https://arxiv.org/html/2607.22004#bib.bib260)\)\. PINNs minimize strong\-form residuals\(Raissiet al\.,[2019](https://arxiv.org/html/2607.22004#bib.bib317)\), while the deep Galerkin method uses a related residual formulation\(Sirignano and Spiliopoulos,[2018](https://arxiv.org/html/2607.22004#bib.bib284)\)\. The deep Ritz method minimizes a variational energy\(E and Yu,[2018](https://arxiv.org/html/2607.22004#bib.bib229)\)\. Other neural PDE solvers include deep BSDE methods, deep splitting methods, and Fourier neural operators\(Hanet al\.,[2018](https://arxiv.org/html/2607.22004#bib.bib310); Eet al\.,[2017](https://arxiv.org/html/2607.22004#bib.bib228); Liet al\.,[2021](https://arxiv.org/html/2607.22004#bib.bib338)\)\. Recent reviews survey the broader field\(Becket al\.,[2020](https://arxiv.org/html/2607.22004#bib.bib297); Weinanet al\.,[2021](https://arxiv.org/html/2607.22004#bib.bib337)\)\.

Training neural PDE solvers to high accuracy remains difficult\. Stiff residual losses and poor conditioning can slow first\-order optimization\(Wanget al\.,[2021](https://arxiv.org/html/2607.22004#bib.bib322); Krishnapriyanet al\.,[2021](https://arxiv.org/html/2607.22004#bib.bib324)\)\. Loss weighting and adaptive residual sampling address part of the problem\(Wanget al\.,[2021](https://arxiv.org/html/2607.22004#bib.bib322); van der Meeret al\.,[2022](https://arxiv.org/html/2607.22004#bib.bib334); Wanget al\.,[2022b](https://arxiv.org/html/2607.22004#bib.bib323)\)\. Curricula and related training strategies provide further controls\(Luet al\.,[2021](https://arxiv.org/html/2607.22004#bib.bib20); Nabianet al\.,[2021](https://arxiv.org/html/2607.22004#bib.bib328); Dawet al\.,[2022](https://arxiv.org/html/2607.22004#bib.bib23)\)\. Other studies examine residual imbalance and training failure modes\(Zapfet al\.,[2022](https://arxiv.org/html/2607.22004#bib.bib22); Wanget al\.,[2022a](https://arxiv.org/html/2607.22004#bib.bib21); Wuet al\.,[2023](https://arxiv.org/html/2607.22004#bib.bib24)\)\. Greedy methods, saddle\-point formulations, and particle\-swarm methods offer alternatives to direct gradient optimization\(Haoet al\.,[2021](https://arxiv.org/html/2607.22004#bib.bib18); Zenget al\.,[2022](https://arxiv.org/html/2607.22004#bib.bib19); Davi and Braga\-Neto,[2022](https://arxiv.org/html/2607.22004#bib.bib327)\)\.

Second\-order methods instead change the geometry of the update\. Energy natural gradient descent \(ENGD\) pulls function\-space energy curvature back to parameter space\(Müller and Zeinhofer,[2023](https://arxiv.org/html/2607.22004#bib.bib4),[2024](https://arxiv.org/html/2607.22004#bib.bib3)\)\. Related PDE\-constrained methods use mass or stiffness matrices as function\-space Gramians\(Schwedeset al\.,[2016](https://arxiv.org/html/2607.22004#bib.bib27),[2017](https://arxiv.org/html/2607.22004#bib.bib26)\)\. Sobolev, Fisher–Rao, and Wasserstein natural gradients have also been studied for PINNs\(Nurbekyanet al\.,[2022](https://arxiv.org/html/2607.22004#bib.bib224)\)\. Gauss–Newton natural gradients and Kronecker\-factored curvature provide further approximations\(Jniniet al\.,[2024](https://arxiv.org/html/2607.22004#bib.bib5); Dangelet al\.,[2024](https://arxiv.org/html/2607.22004#bib.bib6)\)\.

For quadratic residual energies, the Woodbury identity moves the linear solve from parameter space to sample space\. MinSR uses an analogous sample\-space construction in variational Monte Carlo\(Chen and Heyl,[2023](https://arxiv.org/html/2607.22004#bib.bib7); Rendeet al\.,[2024](https://arxiv.org/html/2607.22004#bib.bib12)\)\. SPRING momentum and randomized Nyström sketches have been adapted to PINNs\(Goldshlageret al\.,[2024](https://arxiv.org/html/2607.22004#bib.bib13); Frangellaet al\.,[2023](https://arxiv.org/html/2607.22004#bib.bib14)\)\. Classical Nyström methods construct low\-rank positive\-semidefinite kernel approximations\(Gittens and Mahoney,[2016](https://arxiv.org/html/2607.22004#bib.bib15)\)\. Recent work extends Nyström constructions to Riemannian manifolds\(Nieet al\.,[2026](https://arxiv.org/html/2607.22004#bib.bib16)\)\.

Existing ENGD formulations assume an unconstrained Euclidean parameter domain\. Some neural PDE models impose parameter constraints whose feasible values form a Riemannian manifoldℳ\\mathcal\{M\}\. Atx∈ℳx\\in\\mathcal\{M\}, the realization map sends allowable tangent directions into function space:

Tx​ℳ→d​Pxd​Px​\(Tx​ℳ\)⊂X\.T\_\{x\}\\mathcal\{M\}\\xrightarrow\{\\ dP\_\{x\}\\ \}dP\_\{x\}\(T\_\{x\}\\mathcal\{M\}\)\\subset X\.
The energy Hessian defines the curvature of the resulting function\-space changes\. An ambient ENGD step followed by projection does not generally minimize the constrained quadratic model\. For an ambient curvature operatorAxA\_\{x\}, gradientgxg\_\{x\}, and tangent projectorΠx\\Pi\_\{x\}, one generally has

Πx​Ax−1​gx≠\(Πx​Ax​Πx\|Tx​ℳ\)−1​Πx​gx\.\\Pi\_\{x\}A\_\{x\}^\{\-1\}g\_\{x\}\\neq\\left\(\\left\.\\Pi\_\{x\}A\_\{x\}\\Pi\_\{x\}\\right\|\_\{T\_\{x\}\\mathcal\{M\}\}\\right\)^\{\-1\}\\Pi\_\{x\}g\_\{x\}\.\(1\)
An external projection restores feasibility but can change the minimizer of the tangent quadratic model\. Penalty formulations also change the PDE energy and energy curvature\.

The mismatch raises a natural question:How can an energy natural\-gradient method respect a parameter manifold without changing the PDE energy?

To answer the question, we proposeEnergyManifoldNaturalGradientDescent \(EMNGD\)\. Figure[1](https://arxiv.org/html/2607.22004#S1.F1)depicts EMNGD on the energy landscape overℳ\\mathcal\{M\}\. The left view traces feasible iterates fromθ0\\theta\_\{0\}toward the low\-energy solutionθ∗\\theta^\{\\ast\}, while the inset shows the local update atxkx\_\{k\}\. The energy\-metric solve producesηxk∈Txk​ℳ\\eta\_\{x\_\{k\}\}\\in T\_\{x\_\{k\}\}\\mathcal\{M\}, andRxkR\_\{x\_\{k\}\}maps the tangent pointx~k\+1=xk−αk​ηxk\\widetilde\{x\}\_\{k\+1\}=x\_\{k\}\-\\alpha\_\{k\}\\eta\_\{x\_\{k\}\}to the feasible iteratexk\+1x\_\{k\+1\}\. The construction separates parameter constraints from function\-space energy geometry\.

![Refer to caption](https://arxiv.org/html/2607.22004v1/x1.png)Figure 1:Schematic of an EMNGD update on a parameter manifold\. A tangent step−αk​ηxk\-\\alpha\_\{k\}\\eta\_\{x\_\{k\}\}atxkx\_\{k\}is retracted to the feasible iteratexk\+1x\_\{k\+1\}\.##### Contributions\.

The contributions are as follows:

- •Intrinsic energy manifold geometry\.EMNGD extends ENGD from an unconstrained Euclidean parameter domain to a constrained Riemannian parameter manifold, aligning the parameter geometry with PDE residual constraints\. The energy\-induced quadratic model is defined directly over feasible tangent directions, while retraction\-based updates preserve the parameter constraints\. The intrinsic construction retains the original residual energy and differs from post\-hoc projection of an ambient ENGD step\.
- •Best\-admissible Newton correction\.The main theoretical result characterizes EMNGD as the best admissible approximation to the function\-space Newton correction under the energy metric\. The admissible correction is restricted jointly by the neural realization map and the tangent space of the parameter manifold\. For quadratic energies, the natural\-gradient vector represents the projected current solution error, while the negative retracted step moves toward the corresponding projected solution correction\.
- •Geometric consistency and convergence\.Positive damping makes the tangent energy metric positive definite and yields a unique EMNGD direction\. Coordinate invariance and exact reduction to ENGD in Euclidean space establish consistency across parameter representations\. Under the stated metric\-equivalence and retraction\-smoothness assumptions, Armijo backtracking yields global first\-order convergence\. Controlled inexact tangent solves also preserve the descent property\.
- •Tangent\-space scalable solvers\.For quadratic residual energies and generalized Gauss–Newton pullbacks, exact Woodbury duality transfers the EMNGD tangent system to sample space without changing the damped direction\. Nyström sketch\-and\-solve provides a low\-rank approximate direction, while Nyström preconditioning recovers the exact direction after iterative convergence\. An explicit error bound connects kernel approximation quality and damping with the accuracy of the computed tangent direction\.
- •Scalability with direction control\.Numerical studies verify the Euclidean reduction and the primal–dual agreement of the Woodbury implementation\. Large\-sample diagnostics quantify the effects of sketch rank, damping, and residual subsampling on direction error, convergence, memory consumption, and runtime\. The results identify the regimes in which scalable solvers retain EMNGD accuracy and the regimes in which approximation or sampling error becomes dominant\.

##### Notation\.

We denote the space ofpp\-integrable functions onΩ⊆ℝd\\Omega\\subseteq\\mathbb\{R\}^\{d\}byLp​\(Ω\)L^\{p\}\(\\Omega\)and use the canonical norm ofLp​\(Ω\)L^\{p\}\(\\Omega\)\. For a sufficiently smooth functionuu, let∂iu=∂u/∂xi\\partial\_\{i\}u=\\partial u/\\partial x\_\{i\}\. Let\(Dl​u\)i1,…,il≔∂i1…​∂ilu\(D^\{l\}u\)\_\{i\_\{1\},\\dots,i\_\{l\}\}\\coloneqq\\partial\_\{i\_\{1\}\}\\dots\\partial\_\{i\_\{l\}\}udenote thellth\-derivative tensor\. Let∇u=\(∂1u,…,∂du\)⊤\\nabla u=\(\\partial\_\{1\}u,\\dots,\\partial\_\{d\}u\)^\{\\top\}denote the gradient\. Define the Laplace operator byΔ​u≔∑i=1d∂i2u\\Delta u\\coloneqq\\sum\_\{i=1\}^\{d\}\\partial\_\{i\}^\{2\}u\. We denote the*Sobolev space*of functions with weak derivatives up to orderkkinLp​\(Ω\)L^\{p\}\(\\Omega\)byWk,p​\(Ω\)W^\{k,p\}\(\\Omega\), which is a Banach space with the norm

∥u∥Wk,p​\(Ω\)p≔∑l=0k∥Dl​u∥Lp​\(Ω\)p,\\lVert u\\rVert\_\{W^\{k,p\}\(\\Omega\)\}^\{p\}\\coloneqq\\sum\_\{l=0\}^\{k\}\\lVert D^\{l\}u\\rVert\_\{L^\{p\}\(\\Omega\)\}^\{p\},in the following, we mostly work with the casep=2p=2and writeHk​\(Ω\)H^\{k\}\(\\Omega\)instead ofWk,2​\(Ω\)W^\{k,2\}\(\\Omega\)\.

Letd,m,L,N0,…,NLd,m,L,N\_\{0\},\\dots,N\_\{L\}be natural numbers\. Letθ=\(\(A1,b1\),…,\(AL,bL\)\)\\theta=\\left\(\(A\_\{1\},b\_\{1\}\),\\dots,\(A\_\{L\},b\_\{L\}\)\\right\), whereAl∈ℝNl×Nl−1A\_\{l\}\\in\\mathbb\{R\}^\{N\_\{l\}\\times N\_\{l\-1\}\},bl∈ℝNlb\_\{l\}\\in\\mathbb\{R\}^\{N\_\{l\}\},N0=dN\_\{0\}=d, andNL=mN\_\{L\}=m\. Each pair\(Al,bl\)\(A\_\{l\},b\_\{l\}\)defines an affine mapTl:ℝNl−1→ℝNlT\_\{l\}\\colon\\mathbb\{R\}^\{N\_\{l\-1\}\}\\to\\mathbb\{R\}^\{N\_\{l\}\}\. Given an activation functionρ:ℝ→ℝ\\rho\\colon\\mathbb\{R\}\\to\\mathbb\{R\}, the*neural network function with parameters*θ\\thetais

uθ:ℝd→ℝm,x↦TL​\(ρ​\(TL−1​\(ρ​\(⋯​ρ​\(T1​\(x\)\)\)\)\)\)\.u\_\{\\theta\}\\colon\\mathbb\{R\}^\{d\}\\to\\mathbb\{R\}^\{m\},\\quad x\\mapsto T\_\{L\}\(\\rho\(T\_\{L\-1\}\(\\rho\(\\cdots\\rho\(T\_\{1\}\(x\)\)\)\)\)\)\.
The*number of trainable parameters*of such a network is∑l=0L−1\(nl\+1\)​nl\+1\\sum\_\{l=0\}^\{L\-1\}\(n\_\{l\}\+1\)n\_\{l\+1\}\. We call a network with depth22*shallow*and a deeper network*deep*\. In the remainder, we restrict ourselves to the casem=1m=1since we only consider real\-valued functions\. Our experiments usetanh\\tanhactivations, which are required for the smoothness ofuθu\_\{\\theta\}and the mapθ↦uθ\\theta\\mapsto u\_\{\\theta\}\. ForA∈ℝn×mA\\in\\mathbb\{R\}^\{n\\times m\}, we denote any pseudo inverse ofAAbyA\+A^\{\+\}\.

## 2Preliminaries

Various neural solvers for the approximate solution of PDEs have been suggested\(Becket al\.,[2020](https://arxiv.org/html/2607.22004#bib.bib297); Weinanet al\.,[2021](https://arxiv.org/html/2607.22004#bib.bib337); Kovachkiet al\.,[2021](https://arxiv.org/html/2607.22004#bib.bib344)\)\. Neural PDE solvers parameterize an approximate solution and minimize either a residual energy or a variational energy\. The preliminary discussion introduces both formulations, fixes a common function\-space setup, and summarizes the optimization motivation for natural gradients\.

##### Residual\-form neural PDE solvers\.

Residual\-form neural PDE solvers minimize the PDE residual and boundary mismatch\. Consider a general partial differential equation of the form

ℒ​u=fin​Ωℬ​u=gon​∂Ω,\\displaystyle\\begin\{split\}\\mathcal\{L\}u&=f\\quad\\text\{in \}\\Omega\\\\ \\mathcal\{B\}u&=g\\quad\\text\{on \}\\partial\\Omega,\\end\{split\}\(2\)whereΩ⊆ℝd\\Omega\\subseteq\\mathbb\{R\}^\{d\}is open,ℒ\\mathcal\{L\}is a possibly nonlinear partial differential operator, andℬ\\mathcal\{B\}is a boundary\-value operator\. We seekuuin a Hilbert spaceXX\. Assume thatffis square integrable onΩ\\Omegaandggis square integrable on∂Ω\\partial\\Omega\. Equation \([2](https://arxiv.org/html/2607.22004#S2.E2)\) then admits the minimization formulation

E​\(u\)=∫Ω\(ℒ​u−f\)2​dx\+τ​∫∂Ω\(ℬ​u−g\)2​ds,E\(u\)=\\int\_\{\\Omega\}\(\\mathcal\{L\}u\-f\)^\{2\}\\mathrm\{d\}x\+\\tau\\int\_\{\\partial\\Omega\}\(\\mathcal\{B\}u\-g\)^\{2\}\\mathrm\{d\}s,\(3\)for a penalization parameterτ\>0\\tau\>0\. A functionu∈Xu\\in Xsolves \([2](https://arxiv.org/html/2607.22004#S2.E2)\) exactly whenE​\(u\)=0E\(u\)=0\. For an approximate solution, parameterizeuθu\_\{\\theta\}by a neural network and minimize the parametersθ∈ℝp\\theta\\in\\mathbb\{R\}^\{p\}using

L​\(θ\)≔∫Ω\(ℒ​uθ−f\)2​dx\+τ​∫∂Ω\(ℬ​uθ−g\)2​ds\.L\(\\theta\)\\coloneqq\\int\_\{\\Omega\}\(\\mathcal\{L\}u\_\{\\theta\}\-f\)^\{2\}\\mathrm\{d\}x\+\\tau\\int\_\{\\partial\\Omega\}\\mathcal\{\(\}\\mathcal\{B\}u\_\{\\theta\}\-g\)^\{2\}\\mathrm\{d\}s\.\(4\)
Residual minimization for neural PDE solvers traces back to\(Dissanayake and Phan\-Thien,[1994](https://arxiv.org/html/2607.22004#bib.bib243); Lagariset al\.,[1998](https://arxiv.org/html/2607.22004#bib.bib260)\)\. The deep Galerkin method and physics\-informed neural networks use related residual objectives\(Sirignano and Spiliopoulos,[2018](https://arxiv.org/html/2607.22004#bib.bib284); Raissiet al\.,[2019](https://arxiv.org/html/2607.22004#bib.bib317)\)\. Data terms can be added to the loss\. Numerical implementations discretize the integrals with interior and boundary samples\.

##### Variational neural PDE solvers\.

Weak PDE formulations often use an energy functional whose Euler–Lagrange equations recover the weak form\.Ritz \([1909](https://arxiv.org/html/2607.22004#bib.bib247)\)used the idea to compute polynomial approximation coefficients\.E and Yu \([2018](https://arxiv.org/html/2607.22004#bib.bib229)\)introduced the name*deep Ritz method*for neural networks\. Given a variational energyE:X→ℝE\\colon X\\to\\mathbb\{R\}on a Hilbert spaceXX, parameterize the ansatz byuθu\_\{\\theta\}\. The loss isL​\(θ\)≔E​\(uθ\)L\(\\theta\)\\coloneqq E\(u\_\{\\theta\}\)\. For−Δ​u=f\-\\Delta u=f, the residual energy isu↦∥Δ​u\+f∥L2​\(Ω\)2u\\mapsto\\lVert\\Delta u\+f\\rVert\_\{L^\{2\}\(\\Omega\)\}^\{2\}\. The variational energy isu↦12​∥∇u∥L2​\(Ω\)2−∫Ωf​u​dxu\\mapsto\\frac\{1\}\{2\}\\lVert\\nabla u\\rVert\_\{L^\{2\}\(\\Omega\)\}^\{2\}\-\\int\_\{\\Omega\}fu\\mathrm\{d\}x\. The two energies require different smoothness and belong to different Sobolev spaces\.

Essential boundary values enter the deep Ritz method differently from PINNs\. For PINNs, the unique minimizer is the PDE solution for everyτ\>0\\tau\>0\. In the deep Ritz method, the penalized minimizer solves a perturbed Robin problem\. Accurate approximation of the original problem requires large penalty parameters, which cause ill\-conditioning\(Müller and Zeinhofer,[2022a](https://arxiv.org/html/2607.22004#bib.bib25); Courte and Zeinhofer,[2023](https://arxiv.org/html/2607.22004#bib.bib112)\)\.

##### Function\-space setup\.

Residual and variational formulations minimize an energyE:X→ℝE\\colon X\\to\\mathbb\{R\}\. The parameter loss isL​\(θ\)≔E​\(uθ\)L\(\\theta\)\\coloneqq E\(u\_\{\\theta\}\)\. Assume thatXXis a Hilbert space,uθ∈Xu\_\{\\theta\}\\in X, andEEhas a unique minimizeru∗∈Xu^\{\*\}\\in X\. And assume thatP:ℝp→XP\\colon\\mathbb\{R\}^\{p\}\\to X,θ↦uθ\\theta\\mapsto u\_\{\\theta\}, is differentiable\. DefineℱΘ=\{uθ:θ∈ℝp\}\\mathcal\{F\}\_\{\\Theta\}=\\\{u\_\{\\theta\}:\\theta\\in\\mathbb\{R\}^\{p\}\\\}\. The generalized tangent space ofℱΘ\\mathcal\{F\}\_\{\\Theta\}is

Tθ​ℱΘ≔span⁡\{∂θiuθ:i=1,…,p\}\.T\_\{\\theta\}\\mathcal\{F\}\_\{\\Theta\}\\coloneqq\\operatorname\{span\}\\left\\\{\\partial\_\{\\theta\_\{i\}\}u\_\{\\theta\}:i=1,\\dots,p\\right\\\}\.\(5\)

##### Optimization challenge\.

First\-order optimization can stagnate on residual losses, even for simple PDEs\. Residual stiffness contributes to poor conditioning\(Wanget al\.,[2021](https://arxiv.org/html/2607.22004#bib.bib322)\)\. Squared residuals can further increase the condition number\(Zenget al\.,[2022](https://arxiv.org/html/2607.22004#bib.bib19)\)\. Poor conditioning slows iterative solvers such as gradient descent\. Figure[2](https://arxiv.org/html/2607.22004#S2.F2)illustrates the challenge in one dimension\. Across Poisson, heat, and nonlinear equations, SGD, Adam, BFGS, L\-BFGS, and Adam–L\-BFGS either stagnate at large relativeL2L^\{2\}errors or need many iterations\.

![Refer to caption](https://arxiv.org/html/2607.22004v1/x2.png)Figure 2:RelativeL2L^\{2\}errors of standard optimizers on one\-dimensional PDE benchmarks\.
##### Natural Gradient Descent\.

Amari \([1996](https://arxiv.org/html/2607.22004#bib.bib348)\)originally proposed natural gradient descent \(NGD\) for Euclidean parameter optimization\. Given a positive\-definite metricG​\(θ\)G\(\\theta\), NGD solvesG​\(θ\)​dθ=∇L​\(θ\)G\(\\theta\)d\_\{\\theta\}=\\nabla L\(\\theta\)and takes−dθ\-d\_\{\\theta\}as the descent direction\. In statistical models,G​\(θ\)G\(\\theta\)is usually the Fisher information matrix\. The metric accounts for local model geometry and yields a coordinate\-invariant direction whenG​\(θ\)G\(\\theta\)transforms as a pullback metric\.Müller and Zeinhofer \([2023](https://arxiv.org/html/2607.22004#bib.bib4)\)introduced ENGD for neural PDE solvers by replacing the Fisher metric with an energy\-induced metric\. Foruθ=P​\(θ\)u\_\{\\theta\}=P\(\\theta\), the energy Gram matrix has entriesGE​\(θ\)i​j=D2​E​\(uθ\)​\(∂θiuθ,∂θjuθ\)G\_\{E\}\(\\theta\)\_\{ij\}=D^\{2\}E\(u\_\{\\theta\}\)\(\\partial\_\{\\theta\_\{i\}\}u\_\{\\theta\},\\partial\_\{\\theta\_\{j\}\}u\_\{\\theta\}\)\. The damped system\(GE​\(θ\)\+λ​I\)​d=∇L​\(θ\)\(G\_\{E\}\(\\theta\)\+\\lambda I\)d=\\nabla L\(\\theta\)defines the ENGD direction\. The energy metric gives the direction a direct function\-space interpretation\. The following sections extend the construction to constrained parameter manifolds\.

##### Manifold Optimization\.

Let\(ℳ,g0\)\(\\mathcal\{M\},g^\{0\}\)be a smooth Riemannian manifold and consider the minimization of a differentiable objectiveF:ℳ→ℝF\\colon\\mathcal\{M\}\\to\\mathbb\{R\}\. Atx∈ℳx\\in\\mathcal\{M\}, the tangent spaceTx​ℳT\_\{x\}\\mathcal\{M\}contains the feasible local directions\. The metricgx0g\_\{x\}^\{0\}defines an inner product on that tangent space\. The Riemannian gradientgrad0⁡F​\(x\)∈Tx​ℳ\\operatorname\{grad\}^\{0\}F\(x\)\\in T\_\{x\}\\mathcal\{M\}is defined by

gx0​\(grad0⁡F​\(x\),ξ\)=d​Fx​\[ξ\]for all​ξ∈Tx​ℳ\.g\_\{x\}^\{0\}\(\\operatorname\{grad\}^\{0\}F\(x\),\\xi\)=dF\_\{x\}\[\\xi\]\\qquad\\text\{for all \}\\xi\\in T\_\{x\}\\mathcal\{M\}\.\(6\)
A retractionRx:Tx​ℳ→ℳR\_\{x\}\\colon T\_\{x\}\\mathcal\{M\}\\to\\mathcal\{M\}maps a tangent vector back to the manifold\. A Riemannian optimization step first computesηx∈Tx​ℳ\\eta\_\{x\}\\in T\_\{x\}\\mathcal\{M\}and then sets

xk\+1=Rxk​\(−αk​ηxk\)\.x\_\{k\+1\}=R\_\{x\_\{k\}\}\(\-\\alpha\_\{k\}\\eta\_\{x\_\{k\}\}\)\.\(7\)
Manifold optimization is useful when parameters satisfy hard constraints\. The realization map sends a tangent direction to a first\-order change in function space\.

## 3Energy manifold natural gradient descent \(EMNGD\)

We next develop the geometric formulation of energy natural gradient descent\. Classical ENGD uses function\-space energy curvature on the tangent space of the current neural model\. The Euclidean formulation has several limitations\. The curvature system can be expensive and ill\-conditioned\. The Newton interpretation is local to the current tangent space\. Damping, pseudoinverses, or least\-squares solves are often needed for numerical stability\. Existing scalable ENGD variants only change how that system is solved or approximated\. Examples include kernel, dual, randomized, and low\-rank linear algebra\. Such variants do not encode hard constraints, quotient symmetries, or other feasible\-set geometries\. The manifold formulation begins with the tangent spaceTx​ℳT\_\{x\}\\mathcal\{M\}\. The pullback energy Hessian acts on that space, and a retraction follows each step\.Amari \([1996](https://arxiv.org/html/2607.22004#bib.bib348)\)popularized natural gradients for parameter estimation in supervised learning and blind source separation\. Natural gradients use a chosen metric to define update directions, including Fisher, product\-Fisher, Wasserstein, and Sobolev geometries\(Kakade,[2001](https://arxiv.org/html/2607.22004#bib.bib96); Li and Montúfar,[2018](https://arxiv.org/html/2607.22004#bib.bib31); Nurbekyanet al\.,[2022](https://arxiv.org/html/2607.22004#bib.bib224)\), and have been applied to reinforcement learning\(Kakade,[2001](https://arxiv.org/html/2607.22004#bib.bib96); Peterset al\.,[2003](https://arxiv.org/html/2607.22004#bib.bib41); Bagnell and Schneider,[2003](https://arxiv.org/html/2607.22004#bib.bib38); Morimuraet al\.,[2008](https://arxiv.org/html/2607.22004#bib.bib121)\), inverse problems\(Nurbekyanet al\.,[2022](https://arxiv.org/html/2607.22004#bib.bib224)\), neural\-network training\(Schraudolph,[2002](https://arxiv.org/html/2607.22004#bib.bib342); Pascanu and Bengio,[2014](https://arxiv.org/html/2607.22004#bib.bib340); Martens,[2020](https://arxiv.org/html/2607.22004#bib.bib159)\), and generative models\(Shenet al\.,[2020](https://arxiv.org/html/2607.22004#bib.bib343); Linet al\.,[2021](https://arxiv.org/html/2607.22004#bib.bib341)\)\. A key issue for natural gradients is the choice of function\-space geometry\. The geometry can be defined axiomatically or through the Hessian of a potential function\(Amari and Cichocki,[2010](https://arxiv.org/html/2607.22004#bib.bib201); Amari,[2016](https://arxiv.org/html/2607.22004#bib.bib208); Wang and Yan,[2022](https://arxiv.org/html/2607.22004#bib.bib206); Müller and Montúfar,[2022](https://arxiv.org/html/2607.22004#bib.bib28)\)\. EMNGD uses the exact function\-space Hessian when that Hessian is positive semidefinite on realized directions\. For residual objectives, the implementation may instead use the generalized Gauss–Newton \(GGN\) curvature\.

ForE​\(u\)=12​‖𝒬​\(u\)‖2E\(u\)=\\tfrac\{1\}\{2\}\\\|\\mathcal\{Q\}\(u\)\\\|^\{2\}, the exact Hessian is

D2​E​\(u\)​\[v,w\]=⟨D​𝒬​\(u\)​v,D​𝒬​\(u\)​w⟩\+⟨𝒬​\(u\),D2​𝒬​\(u\)​\[v,w\]⟩\.D^\{2\}E\(u\)\[v,w\]=\\langle D\\mathcal\{Q\}\(u\)v,D\\mathcal\{Q\}\(u\)w\\rangle\+\\langle\\mathcal\{Q\}\(u\),D^\{2\}\\mathcal\{Q\}\(u\)\[v,w\]\\rangle\.
The GGN retains the first term\. The approximation equals the exact Hessian when𝒬\\mathcal\{Q\}is affine\. For nonlinear𝒬\\mathcal\{Q\}, the approximation discards the residual\-weighted second derivative\. Related curvature methods have been proposed for supervised neural\-network training\(Ren and Goldfarb,[2019](https://arxiv.org/html/2607.22004#bib.bib346); Caiet al\.,[2019](https://arxiv.org/html/2607.22004#bib.bib345); Gargianiet al\.,[2020](https://arxiv.org/html/2607.22004#bib.bib347); Martens,[2020](https://arxiv.org/html/2607.22004#bib.bib159)\)\. Our applications may involve infinite\-dimensional or non\-strongly\-convex objectives\.

###### Assumption 1\(Geometric and analytic setting\)\.

Let\(ℳ,g0\)\(\\mathcal\{M\},g^\{0\}\)be a finite\-dimensional smooth Riemannian manifold\. And letRRbe a retraction, i\.e\.,

Rx​\(0x\)=x,D​Rx​\(0x\)=idTx​ℳ\.R\_\{x\}\(0\_\{x\}\)=x,\\qquad DR\_\{x\}\(0\_\{x\}\)=\\mathrm\{id\}\_\{T\_\{x\}\\mathcal\{M\}\}\.
LetXXbe a real Hilbert space,E:X→ℝE:X\\to\\mathbb\{R\}be twice Fréchet differentiable, andP:ℳ→XP:\\mathcal\{M\}\\to Xbe twice differentiable\. We minimizeF=E∘PF=E\\circ Ponℳ\\mathcal\{M\}\. The differentialJx=d​Px:Tx​ℳ→XJ\_\{x\}=dP\_\{x\}:T\_\{x\}\\mathcal\{M\}\\to Xmaps a tangent direction to function space\. The baseline Riemannian gradient is defined by

gx0​\(grad0⁡F​\(x\),ζ\)=d​Fx​\[ζ\],ζ∈Tx​ℳ\.g\_\{x\}^\{0\}\(\\operatorname\{grad\}^\{0\}F\(x\),\\zeta\)=dF\_\{x\}\[\\zeta\],\\qquad\\zeta\\in T\_\{x\}\\mathcal\{M\}\.

The parameter manifoldℳ\\mathcal\{M\}and the imageJx​\(Tx​ℳ\)⊂XJ\_\{x\}\(T\_\{x\}\\mathcal\{M\}\)\\subset Xhave different roles\. The parameter manifold defines the allowed directions\. The image contains the corresponding first\-order changes in function space\. The energyE:X→ℝE\\colon X\\to\\mathbb\{R\}is twice differentiable\. The setting covers both PINNs and the deep Ritz method\. The energy Hessian induces on each tangent space the pullback bilinear form

gxE​\(ξ,ζ\)≔D2​E​\(P​\(x\)\)​\(Jx​ξ,Jx​ζ\),ξ,ζ∈Tx​ℳ\.g\_\{x\}^\{E\}\(\\xi,\\zeta\)\\coloneqq D^\{2\}E\(P\(x\)\)\(J\_\{x\}\\xi,J\_\{x\}\\zeta\),\\qquad\\xi,\\zeta\\in T\_\{x\}\\mathcal\{M\}\.\(8\)
IfgxEg\_\{x\}^\{E\}is positive definite, the bilinear form defines a Riemannian metric onℳ\\mathcal\{M\}\. IfgxEg\_\{x\}^\{E\}is not positive definite, we use the damped form

gxE,λ​\(ξ,ζ\)≔gxE​\(ξ,ζ\)\+λ​gx0​\(ξ,ζ\),λ≥0,g\_\{x\}^\{E,\\lambda\}\(\\xi,\\zeta\)\\coloneqq g\_\{x\}^\{E\}\(\\xi,\\zeta\)\+\\lambda g\_\{x\}^\{0\}\(\\xi,\\zeta\),\\qquad\\lambda\\geq 0,\(9\)forλ\>0\\lambda\>0, damping gives a regularized approximation to the minimum\-norm pseudoinverse solution\. The least\-squares implementations use the same regularized tangent system\. The EMNGD direction is the tangent vectorηx∈Tx​ℳ\\eta\_\{x\}\\in T\_\{x\}\\mathcal\{M\}satisfying

gxE,λ​\(ηx,ζ\)=d​Fx​\[ζ\]=gx0​\(grad0⁡F​\(x\),ζ\)for all​ζ∈Tx​ℳ\.g\_\{x\}^\{E,\\lambda\}\(\\eta\_\{x\},\\zeta\)=dF\_\{x\}\[\\zeta\]=g\_\{x\}^\{0\}\(\\operatorname\{grad\}^\{0\}F\(x\),\\zeta\)\\qquad\\text\{for all \}\\zeta\\in T\_\{x\}\\mathcal\{M\}\.\(10\)
The equation is a linear system on the tangent space\. EMNGD then updates by

xk\+1=Rxk​\(−αk​ηxk\)\.x\_\{k\+1\}=R\_\{x\_\{k\}\}\(\-\\alpha\_\{k\}\\eta\_\{x\_\{k\}\}\)\.\(11\)
In local coordinatesx=ϕ​\(ξ\)x=\\phi\(\\xi\), the Euclidean energy Gram matrix becomes the pullbackG~E​\(ξ\)=Jϕ​\(ξ\)⊤​GE​\(ϕ​\(ξ\)\)​Jϕ​\(ξ\)\\widetilde\{G\}\_\{E\}\(\\xi\)=J\_\{\\phi\}\(\\xi\)^\{\\top\}G\_\{E\}\(\\phi\(\\xi\)\)J\_\{\\phi\}\(\\xi\)\. The unconstrained parameter case used in our experiments corresponds toℳ=ℝp\\mathcal\{M\}=\\mathbb\{R\}^\{p\},g0g^\{0\}equal to the Euclidean metric, and the retractionRθ​\(v\)=θ\+vR\_\{\\theta\}\(v\)=\\theta\+v\.

###### Proposition 2\(Operator form and variational characterization\)\.

Assume thatD2​E​\(P​\(x\)\)D^\{2\}E\(P\(x\)\)is symmetric positive semidefinite as a bilinear form onXX\. ThengxEg\_\{x\}^\{E\}is symmetric positive semidefinite onTx​ℳT\_\{x\}\\mathcal\{M\}\. Ifλ\>0\\lambda\>0, thengxE,λg\_\{x\}^\{E,\\lambda\}is positive definite and there exists a uniquegx0g\_\{x\}^\{0\}\-self\-adjoint positive definite operator

Axλ:Tx​ℳ→Tx​ℳA\_\{x\}^\{\\lambda\}:T\_\{x\}\\mathcal\{M\}\\to T\_\{x\}\\mathcal\{M\}such that

gx0​\(Axλ​ξ,ζ\)=gxE,λ​\(ξ,ζ\)for all​ξ,ζ∈Tx​ℳ\.g\_\{x\}^\{0\}\(A\_\{x\}^\{\\lambda\}\\xi,\\zeta\)=g\_\{x\}^\{E,\\lambda\}\(\\xi,\\zeta\)\\qquad\\text\{for all \}\\xi,\\zeta\\in T\_\{x\}\\mathcal\{M\}\.\(12\)
The EMNGD direction is uniquely

ηx=\(Axλ\)−1​grad0⁡F​\(x\),\\eta\_\{x\}=\(A\_\{x\}^\{\\lambda\}\)^\{\-1\}\\operatorname\{grad\}^\{0\}F\(x\),\(13\)which is the unique minimizer of

Qx​\(ξ\)=12​gxE,λ​\(ξ,ξ\)−d​Fx​\[ξ\],ξ∈Tx​ℳ\.Q\_\{x\}\(\\xi\)=\\frac\{1\}\{2\}g\_\{x\}^\{E,\\lambda\}\(\\xi,\\xi\)\-dF\_\{x\}\[\\xi\],\\qquad\\xi\\in T\_\{x\}\\mathcal\{M\}\.\(14\)

ProofFor anyξ∈Tx​ℳ\\xi\\in T\_\{x\}\\mathcal\{M\},gxE​\(ξ,ξ\)=D2​E​\(P​\(x\)\)​\[Jx​ξ,Jx​ξ\]≥0g\_\{x\}^\{E\}\(\\xi,\\xi\)=D^\{2\}E\(P\(x\)\)\[J\_\{x\}\\xi,J\_\{x\}\\xi\]\\geq 0\. Ifλ\>0\\lambda\>0andξ≠0\\xi\\neq 0, thengxE,λ​\(ξ,ξ\)≥λ​gx0​\(ξ,ξ\)\>0g\_\{x\}^\{E,\\lambda\}\(\\xi,\\xi\)\\geq\\lambda g\_\{x\}^\{0\}\(\\xi,\\xi\)\>0\. The Riesz representation theorem on the finite\-dimensional inner\-product space\(Tx​ℳ,gx0\)\(T\_\{x\}\\mathcal\{M\},g\_\{x\}^\{0\}\)givesAxλA\_\{x\}^\{\\lambda\}, and symmetry ofgxE,λg\_\{x\}^\{E,\\lambda\}makesAxλA\_\{x\}^\{\\lambda\}self\-adjoint\. Substituting \([12](https://arxiv.org/html/2607.22004#S3.E12)\) into \([10](https://arxiv.org/html/2607.22004#S3.E10)\) givesAxλ​ηx=grad0⁡F​\(x\)A\_\{x\}^\{\\lambda\}\\eta\_\{x\}=\\operatorname\{grad\}^\{0\}F\(x\)\. The first\-order optimality condition forQxQ\_\{x\}is exactly \([10](https://arxiv.org/html/2607.22004#S3.E10)\), and strict convexity follows from positive definiteness\.

###### Theorem 3\(Coordinate form and coordinate invariance\)\.

Letϕ:U⊂ℝq→ℳ\\phi:U\\subset\\mathbb\{R\}^\{q\}\\to\\mathcal\{M\}be a local chart withx=ϕ​\(y\)x=\\phi\(y\)andei=∂iϕ​\(y\)e\_\{i\}=\\partial\_\{i\}\\phi\(y\)\. Writingηx=∑i=1qvi​ei\\eta\_\{x\}=\\sum\_\{i=1\}^\{q\}v\_\{i\}e\_\{i\}, define

\(GEϕ\)i​j=gxE​\(ei,ej\),\(G0ϕ\)i​j=gx0​\(ei,ej\),bi=d​Fx​\[ei\]\.\(G\_\{E\}^\{\\phi\}\)\_\{ij\}=g\_\{x\}^\{E\}\(e\_\{i\},e\_\{j\}\),\\qquad\(G\_\{0\}^\{\\phi\}\)\_\{ij\}=g\_\{x\}^\{0\}\(e\_\{i\},e\_\{j\}\),\\qquad b\_\{i\}=dF\_\{x\}\[e\_\{i\}\]\.
Then the coordinate vectorvvsatisfies

\(GEϕ\+λ​G0ϕ\)​v=b\.\(G\_\{E\}^\{\\phi\}\+\\lambda G\_\{0\}^\{\\phi\}\)v=b\.\(15\)
Forλ\>0\\lambda\>0, the tangent system has a unique solution and the tangent vectorηx\\eta\_\{x\}is independent of the chosen chart\. Forℳ=ℝp\\mathcal\{M\}=\\mathbb\{R\}^\{p\}with the Euclidean metric,P​\(θ\)=uθP\(\\theta\)=u\_\{\\theta\}, andRθ​\(v\)=θ\+vR\_\{\\theta\}\(v\)=\\theta\+v, \([15](https://arxiv.org/html/2607.22004#S3.E15)\) becomes

\(GE​\(θ\)\+λ​I\)​d=∇L​\(θ\)\.\(G\_\{E\}\(\\theta\)\+\\lambda I\)d=\\nabla L\(\\theta\)\.\(16\)
Forλ=0\\lambda=0, the Moore–Penrose convention gives the minimum\-norm pseudoinverse directiond=GE​\(θ\)\+​∇L​\(θ\)d=G\_\{E\}\(\\theta\)^\{\+\}\\nabla L\(\\theta\)when the undamped system is singular\.

ProofTesting \([10](https://arxiv.org/html/2607.22004#S3.E10)\) with each basis vectoreje\_\{j\}gives \([15](https://arxiv.org/html/2607.22004#S3.E15)\)\. Forλ\>0\\lambda\>0, the coefficient matrix is positive definite because

v⊤​\(GEϕ\+λ​G0ϕ\)​v=gxE,λ​\(∑ivi​ei,∑ivi​ei\)\>0,v^\{\\top\}\(G\_\{E\}^\{\\phi\}\+\\lambda G\_\{0\}^\{\\phi\}\)v=g\_\{x\}^\{E,\\lambda\}\\Big\(\\sum\_\{i\}v\_\{i\}e\_\{i\},\\sum\_\{i\}v\_\{i\}e\_\{i\}\\Big\)\>0,\(17\)for every nonzerovv\. The vector reconstructed from the coordinate solution satisfies the weak equation \([10](https://arxiv.org/html/2607.22004#S3.E10)\); uniqueness of that equation implies chart independence\. The Euclidean reduction follows from the identity chart, for whichG0=IG\_\{0\}=Iandb=∇L​\(θ\)b=\\nabla L\(\\theta\)\.

![Refer to caption](https://arxiv.org/html/2607.22004v1/figures/emngd.png)Figure 3:EMNGD update on a parameter manifold\. From left to right, EMNGD solves the energy\-metric equation inTxk​ℳT\_\{x\_\{k\}\}\\mathcal\{M\}, retracts−αk​ηxk\-\\alpha\_\{k\}\\eta\_\{x\_\{k\}\}toxk\+1∈ℳx\_\{k\+1\}\\in\\mathcal\{M\}, and repeats the update towardx∗x^\{\\ast\}\.We define the*Hilbert*and*energy Gram matrices*by

GH​\(θ\)i​j≔⟨∂θiuθ,∂θjuθ⟩X,\\displaystyle G\_\{H\}\(\\theta\)\_\{ij\}\\coloneqq\\langle\\partial\_\{\\theta\_\{i\}\}u\_\{\\theta\},\\partial\_\{\\theta\_\{j\}\}u\_\{\\theta\}\\rangle\_\{X\},\\quad\(18\)and

GE​\(θ\)i​j≔D2​E​\(uθ\)​\(∂θiuθ,∂θjuθ\)\.G\_\{E\}\(\\theta\)\_\{ij\}\\coloneqq D^\{2\}E\(u\_\{\\theta\}\)\(\\partial\_\{\\theta\_\{i\}\}u\_\{\\theta\},\\partial\_\{\\theta\_\{j\}\}u\_\{\\theta\}\)\.\(19\)
The*Hilbert natural\-gradient direction*∇HL​\(θ\)=GH​\(θ\)\+​∇L​\(θ\)\\nabla^\{H\}L\(\\theta\)=G\_\{H\}\(\\theta\)^\{\+\}\\nabla L\(\\theta\)uses the Sobolev inner product⟨⋅,⋅⟩X\\langle\\cdot,\\cdot\\rangle\_\{X\}for neural\-network training\(Nurbekyanet al\.,[2022](https://arxiv.org/html/2607.22004#bib.bib224)\)\. For a Sobolev spaceXX, the direction is also called the*Sobolev natural gradient*; H\-NG denotes the*Hilbert natural gradient*\. Natural\-gradient theory establishes111For regular and singular Gram matrices and finite\-dimensional spaces, see\(Amari,[2016](https://arxiv.org/html/2607.22004#bib.bib208); van Oostrumet al\.,[2022](https://arxiv.org/html/2607.22004#bib.bib182)\)\. The appendix gives an argument for infinite\-dimensional spaces\.that222Here, the Hilbert space gradient∇E​\(u\)∈X\\nabla E\(u\)\\in Xis the unique element satisfying⟨∇E​\(u\),v⟩X=D​E​\(u\)​v\\langle\\nabla E\(u\),v\\rangle\_\{X\}=DE\(u\)v, whereD​EDEdenotes the Fréchet derivative\.

D​Pθ​∇HL​\(θ\)=ΠTθ​ℱΘ​\(∇E​\(uθ\)\)\.DP\_\{\\theta\}\\nabla^\{H\}L\(\\theta\)=\\Pi\_\{T\_\{\\theta\}\\mathcal\{F\}\_\{\\Theta\}\}\(\\nabla E\(u\_\{\\theta\}\)\)\.\(20\)
In words, following the natural gradient amounts to moving along the projection of the Hilbert space gradient onto the model’s tangent space in function space\. The observation that identifying the function space gradient via the Hessian leads to a Newton update motivates the concept of energy natural gradients that we now introduce\.

###### Definition 4\(Energy Manifold Natural Gradient\)\.

Under Assumption[1](https://arxiv.org/html/2607.22004#Thmtheorem1), the*energy manifold natural gradient*atx∈ℳx\\in\\mathcal\{M\}is the tangent vectorηx\\eta\_\{x\}solving \([10](https://arxiv.org/html/2607.22004#S3.E10)\)\. The associated descent direction is−ηx\-\\eta\_\{x\}, and the algorithmic update is the retracted step \([11](https://arxiv.org/html/2607.22004#S3.E11)\)\. In the Euclidean parameter caseℳ=ℝp\\mathcal\{M\}=\\mathbb\{R\}^\{p\},Rθ​\(v\)=θ\+vR\_\{\\theta\}\(v\)=\\theta\+v, andλ=0\\lambda=0, Definition[4](https://arxiv.org/html/2607.22004#Thmtheorem4)reduces to

∇EL​\(θ\)≔GE​\(θ\)\+​∇L​\(θ\),\\nabla^\{E\}L\(\\theta\)\\coloneqq G\_\{E\}\(\\theta\)^\{\+\}\\nabla L\(\\theta\),\(21\)the usual energy natural gradient direction\.

### 3\.1Scalable solvers for the EMNGD tangent system

For a linear PDE operatorℒ\\mathcal\{L\}, the residual yields a quadratic energy, and the energy Gram matrix takes the form

GE​\(θ\)i​j=∫Ωℒ​\(∂θiuθ\)​ℒ​\(∂θjuθ\)​dx\+τ​∫∂Ωℬ​\(∂θiuθ\)​ℬ​\(∂θjuθ\)​ds\\displaystyle\\begin\{split\}G\_\{E\}\(\\theta\)\_\{ij\}&=\\int\_\{\\Omega\}\\mathcal\{L\}\(\\partial\_\{\\theta\_\{i\}\}u\_\{\\theta\}\)\\mathcal\{L\}\(\\partial\_\{\\theta\_\{j\}\}u\_\{\\theta\}\)\\mathrm\{d\}x\\\\ &\+\\tau\\int\_\{\\partial\\Omega\}\\mathcal\{B\}\(\\partial\_\{\\theta\_\{i\}\}u\_\{\\theta\}\)\\mathcal\{B\}\(\\partial\_\{\\theta\_\{j\}\}u\_\{\\theta\}\)\\mathrm\{d\}s\\end\{split\}\(22\)
The residual\-energy expression also exposes the low\-rank structure used by scalable implementations\. After quadrature, the residual loss can be written as

L​\(θ\)=12​‖r​\(θ\)‖22,r​\(θ\)∈ℝN,L\(\\theta\)=\\frac\{1\}\{2\}\\\|r\(\\theta\)\\\|\_\{2\}^\{2\},\\qquad r\(\\theta\)\\in\\mathbb\{R\}^\{N\},\(23\)where the entries ofrrcollect the weighted interior and boundary residuals\. LetJ​\(θ\)=Dθ​r​\(θ\)∈ℝN×pJ\(\\theta\)=D\_\{\\theta\}r\(\\theta\)\\in\\mathbb\{R\}^\{N\\times p\}be the residual Jacobian\. For a linear PDE operator,J⊤​JJ^\{\\top\}Jis the exact pullback of the quadratic function\-space Hessian\. For a nonlinear residual map,J⊤​JJ^\{\\top\}Jis the GGN pullback and omits residual\-weighted second\-derivative terms from the parameter Hessian\. Then

GE​\(θ\)=J​\(θ\)⊤​J​\(θ\),∇L​\(θ\)=J​\(θ\)⊤​r​\(θ\)\.G\_\{E\}\(\\theta\)=J\(\\theta\)^\{\\top\}J\(\\theta\),\\qquad\\nabla L\(\\theta\)=J\(\\theta\)^\{\\top\}r\(\\theta\)\.\(24\)
The damped EMNGD direction satisfies

∇λEL​\(θ\)=\(J⊤​J\+λ​I\)−1​J⊤​r,\\nabla^\{E\}\_\{\\lambda\}L\(\\theta\)=\\big\(J^\{\\top\}J\+\\lambda I\\big\)^\{\-1\}J^\{\\top\}r,\(25\)whereJ=J​\(θ\)J=J\(\\theta\)andr=r​\(θ\)r=r\(\\theta\)\. Applying the push\-through identity, equivalently the Woodbury matrix identity, gives the sample\-space form

\(J⊤​J\+λ​I\)−1​J⊤​r=J⊤​\(J​J⊤\+λ​I\)−1​r\.\\big\(J^\{\\top\}J\+\\lambda I\\big\)^\{\-1\}J^\{\\top\}r=J^\{\\top\}\\big\(JJ^\{\\top\}\+\\lambda I\\big\)^\{\-1\}r\.\(26\)
For an embedded parameter manifold, let𝐉x:ℝp→ℝN\\mathbf\{J\}\_\{x\}:\\mathbb\{R\}^\{p\}\\to\\mathbb\{R\}^\{N\}be the ambient residual Jacobian\. LetΠx:ℝp→Tx​ℳ\\Pi\_\{x\}:\\mathbb\{R\}^\{p\}\\to T\_\{x\}\\mathcal\{M\}be the orthogonal projector onto the tangent space\. The corresponding intrinsic direction is

ηx=Πx​𝐉x⊤​\(𝐉x​Πx​𝐉x⊤\+λ​I\)−1​r​\(x\)\.\\eta\_\{x\}=\\Pi\_\{x\}\\mathbf\{J\}\_\{x\}^\{\\top\}\\big\(\\mathbf\{J\}\_\{x\}\\Pi\_\{x\}\\mathbf\{J\}\_\{x\}^\{\\top\}\+\\lambda I\\big\)^\{\-1\}r\(x\)\.\(27\)
The Woodbury identity computes the same EMNGD tangent direction from anN×NN\\times Nsystem rather than ap×pp\\times psystem\. The reduction is useful when quadrature or collocation samples are far fewer than trainable parameters\. The dominant solve then occurs in sample space\. The matrixJ​J⊤JJ^\{\\top\}is the energy analogue of the empirical neural tangent kernel\(Jacotet al\.,[2018](https://arxiv.org/html/2607.22004#bib.bib272)\)\. Efficient kernel construction techniques can be used without changing the EMNGD geometry\(Novaket al\.,[2022](https://arxiv.org/html/2607.22004#bib.bib17)\)\.

##### Exact Woodbury duality and Nyström sketches\.

Let𝒥x=D​rx:Tx​ℳ→ℝN\\mathcal\{J\}\_\{x\}=Dr\_\{x\}:T\_\{x\}\\mathcal\{M\}\\to\\mathbb\{R\}^\{N\}denote the residual differential, and let𝒥x∗\\mathcal\{J\}\_\{x\}^\{\*\}be the adjoint induced bygx0g\_\{x\}^\{0\}and the Euclidean inner product\. The residual gradient isgrad0⁡F​\(x\)=𝒥x∗​r​\(x\)\\operatorname\{grad\}^\{0\}F\(x\)=\\mathcal\{J\}\_\{x\}^\{\*\}r\(x\)\. For the quadratic residual energy or the corresponding GGN metric, the intrinsic damped direction and sample\-space kernel are

ηx=\(𝒥x∗​𝒥x\+λ​Ix\)−1​𝒥x∗​r​\(x\),Kx=𝒥x​𝒥x∗\.\\eta\_\{x\}=\(\\mathcal\{J\}\_\{x\}^\{\*\}\\mathcal\{J\}\_\{x\}\+\\lambda I\_\{x\}\)^\{\-1\}\\mathcal\{J\}\_\{x\}^\{\*\}r\(x\),\\qquad K\_\{x\}=\\mathcal\{J\}\_\{x\}\\mathcal\{J\}\_\{x\}^\{\*\}\.\(28\)
The push\-through identity givesηx=𝒥x∗​\(Kx\+λ​I\)−1​r​\(x\)\\eta\_\{x\}=\\mathcal\{J\}\_\{x\}^\{\*\}\(K\_\{x\}\+\\lambda I\)^\{\-1\}r\(x\)\. Woodbury is an exact dual representation of the damped tangent direction\. In embedded coordinates,𝒥x=𝐉x​Πx\\mathcal\{J\}\_\{x\}=\\mathbf\{J\}\_\{x\}\\Pi\_\{x\}, which recovers \([27](https://arxiv.org/html/2607.22004#S3.E27)\)\.

###### Proposition 5\(Nyström direction error\)\.

LetK^x⪰0\\widehat\{K\}\_\{x\}\\succeq 0be a rank\-ℓ\\ellNyström approximation ofKxK\_\{x\}, and define

η~x=𝒥x∗​\(K^x\+λ​I\)−1​r​\(x\),λ\>0\.\\widetilde\{\\eta\}\_\{x\}=\\mathcal\{J\}\_\{x\}^\{\*\}\(\\widehat\{K\}\_\{x\}\+\\lambda I\)^\{\-1\}r\(x\),\\qquad\\lambda\>0\.
Then

‖η~x−ηx‖0≤‖𝒥x‖0→2​‖Kx−K^x‖2λ2​‖r​\(x\)‖2\.\\\|\\widetilde\{\\eta\}\_\{x\}\-\\eta\_\{x\}\\\|\_\{0\}\\leq\\frac\{\\\|\\mathcal\{J\}\_\{x\}\\\|\_\{0\\to 2\}\\,\\\|K\_\{x\}\-\\widehat\{K\}\_\{x\}\\\|\_\{2\}\}\{\\lambda^\{2\}\}\\\|r\(x\)\\\|\_\{2\}\.\(29\)

ProofThe resolvent identity gives

\(K^x\+λ​I\)−1−\(Kx\+λ​I\)−1=\(K^x\+λ​I\)−1​\(Kx−K^x\)​\(Kx\+λ​I\)−1\.\(\\widehat\{K\}\_\{x\}\+\\lambda I\)^\{\-1\}\-\(K\_\{x\}\+\\lambda I\)^\{\-1\}=\(\\widehat\{K\}\_\{x\}\+\\lambda I\)^\{\-1\}\(K\_\{x\}\-\\widehat\{K\}\_\{x\}\)\(K\_\{x\}\+\\lambda I\)^\{\-1\}\.
The positive damping is essential: positive semidefiniteness bounds both inverse norms byλ−1\\lambda^\{\-1\}\. Withoutλ\>0\\lambda\>0, the inverse norm can diverge whenKxK\_\{x\}orK^x\\widehat\{K\}\_\{x\}is singular\. Applying𝒥x∗\\mathcal\{J\}\_\{x\}^\{\*\}and the operator\-norm bound gives \([29](https://arxiv.org/html/2607.22004#S3.E29)\)\.

A Riemannian Nyström construction can also approximate the tangent Gram operator𝒢x=𝒥x∗​𝒥x\\mathcal\{G\}\_\{x\}=\\mathcal\{J\}\_\{x\}^\{\*\}\\mathcal\{J\}\_\{x\}directly\. For a rank\-ℓ\\elltangent sketch𝒫x\\mathcal\{P\}\_\{x\}, the intrinsic approximation is

𝒢^x=𝒢x​𝒫x​\(𝒫x∗​𝒢x​𝒫x\)†​𝒫x∗​𝒢x\.\\widehat\{\\mathcal\{G\}\}\_\{x\}=\\mathcal\{G\}\_\{x\}\\mathcal\{P\}\_\{x\}\(\\mathcal\{P\}\_\{x\}^\{\*\}\\mathcal\{G\}\_\{x\}\\mathcal\{P\}\_\{x\}\)^\{\\dagger\}\\mathcal\{P\}\_\{x\}^\{\*\}\\mathcal\{G\}\_\{x\}\.
The intrinsic approximation isgx0g\_\{x\}^\{0\}\-self\-adjoint, positive semidefinite, and has rank at mostℓ\\ell\(Nieet al\.,[2026](https://arxiv.org/html/2607.22004#bib.bib16)\)\. EMNGD instead sketches the dual kernelKxK\_\{x\}after the exact Woodbury transformation\. A sketch\-and\-solve update usesη~x\\widetilde\{\\eta\}\_\{x\}and is approximate\. Equation \([29](https://arxiv.org/html/2607.22004#S3.E29)\) bounds the resulting direction error\. A Nyström preconditioner changes only the conditioning of an iterative solve of\(Kx\+λ​I\)​a=r​\(x\)\(K\_\{x\}\+\\lambda I\)a=r\(x\)\. After convergence,ηx=𝒥x∗​a\\eta\_\{x\}=\\mathcal\{J\}\_\{x\}^\{\*\}ais the exact Woodbury direction\.

Equivalently, \([25](https://arxiv.org/html/2607.22004#S3.E25)\) is the solution of the Tikhonov\-regularized least\-squares problem

∇λEL​\(θ\)=arg⁡minψ∈ℝp⁡\{12​‖J​ψ−r‖22\+λ2​‖ψ‖22\}\.\\nabla^\{E\}\_\{\\lambda\}L\(\\theta\)=\\arg\\min\_\{\\psi\\in\\mathbb\{R\}^\{p\}\}\\left\\\{\\frac\{1\}\{2\}\\\|J\\psi\-r\\\|\_\{2\}^\{2\}\+\\frac\{\\lambda\}\{2\}\\\|\\psi\\\|\_\{2\}^\{2\}\\right\\\}\.\(30\)
The first\-order condition of \([30](https://arxiv.org/html/2607.22004#S3.E30)\) is

\(J⊤​J\+λ​I\)​ψ=J⊤​r,\(J^\{\\top\}J\+\\lambda I\)\\psi=J^\{\\top\}r,which gives \([25](https://arxiv.org/html/2607.22004#S3.E25)\); the identity

\(J⊤​J\+λ​I\)​J⊤=J⊤​\(J​J⊤\+λ​I\)\(J^\{\\top\}J\+\\lambda I\)J^\{\\top\}=J^\{\\top\}\(JJ^\{\\top\}\+\\lambda I\)then gives \([26](https://arxiv.org/html/2607.22004#S3.E26)\)\. Ifr∈range⁡\(J\)r\\in\\operatorname\{range\}\(J\), singular\-value decomposition shows that∇λEL​\(θ\)→J\+​r\\nabla^\{E\}\_\{\\lambda\}L\(\\theta\)\\to J^\{\+\}rasλ↓0\\lambda\\downarrow 0\. The damped direction then converges to the minimum\-norm residual\-matching direction\. SPRING\-style momentum replaces the regularization center with a previous direction\(Goldshlageret al\.,[2024](https://arxiv.org/html/2607.22004#bib.bib13)\)\. The tangent space, energy metric, and retraction remain unchanged\.

For a quadratic energyE​\(u\)=12​a​\(u,u\)−f​\(u\)E\(u\)=\\frac\{1\}\{2\}a\(u,u\)\-f\(u\), the deep Ritz method uses a symmetric and coercive bilinear formaaandf∈X∗f\\in X^\{\*\}, which gives

GE​\(θ\)i​j=a​\(∂θiuθ,∂θjuθ\)\.G\_\{E\}\(\\theta\)\_\{ij\}=a\(\\partial\_\{\\theta\_\{i\}\}u\_\{\\theta\},\\partial\_\{\\theta\_\{j\}\}u\_\{\\theta\}\)\.\(31\)

## 4Theoretical Analysis

The following theorem is the central structural statement of EMNGD\. The theorem states that the EMNGD vector is the best feasible approximation to the function\-space Newton correction under the energy metric\.

###### Theorem 6\(Main theorem: energy\-manifold Newton projection\)\.

Assume thatD2​E​\(P​\(x\)\)D^\{2\}E\(P\(x\)\)is symmetric, bounded, and coercive onXX\. LetHx=D2​E​\(P​\(x\)\)H\_\{x\}=D^\{2\}E\(P\(x\)\)andNx∈XN\_\{x\}\\in Xbe the function\-space Newton vector defined by

Hx​\[Nx,v\]=D​E​\(P​\(x\)\)​\[v\]for all​v∈X\.H\_\{x\}\[N\_\{x\},v\]=DE\(P\(x\)\)\[v\]\\qquad\\text\{for all \}v\\in X\.\(32\)
LetSx=Jx​\(Tx​ℳ\)⊂XS\_\{x\}=J\_\{x\}\(T\_\{x\}\\mathcal\{M\}\)\\subset X\. If the undamped EMNGD equation has a solution, the push\-forward satisfies

Jx​ηx=ΠSxHx​Nx,J\_\{x\}\\eta\_\{x\}=\\Pi\_\{S\_\{x\}\}^\{H\_\{x\}\}N\_\{x\},\(33\)whereΠSxHx\\Pi\_\{S\_\{x\}\}^\{H\_\{x\}\}is theHxH\_\{x\}\-orthogonal projection ontoSxS\_\{x\}\. In Euclidean coordinates, the projection identity becomes

D​Pθ​∇EL​\(θ\)=ΠTθ​ℱΘD2​E​\(uθ\)​\(D2​E​\(uθ\)−1​∇E​\(uθ\)\)\.DP\_\{\\theta\}\\nabla^\{E\}L\(\\theta\)=\\Pi\_\{T\_\{\\theta\}\\mathcal\{F\}\_\{\\Theta\}\}^\{D^\{2\}E\(u\_\{\\theta\}\)\}\\big\(D^\{2\}E\(u\_\{\\theta\}\)^\{\-1\}\\nabla E\(u\_\{\\theta\}\)\\big\)\.

ProofThe chain rule givesd​Fx​\[ζ\]=D​E​\(P​\(x\)\)​\[Jx​ζ\]dF\_\{x\}\[\\zeta\]=DE\(P\(x\)\)\[J\_\{x\}\\zeta\]\. The definition ofNxN\_\{x\}givesHx​\[Nx,Jx​ζ\]H\_\{x\}\[N\_\{x\},J\_\{x\}\\zeta\]\. The undamped EMNGD equation gives

Hx​\[Jx​ηx,Jx​ζ\]=Hx​\[Nx,Jx​ζ\]for all​ζ∈Tx​ℳ\.H\_\{x\}\[J\_\{x\}\\eta\_\{x\},J\_\{x\}\\zeta\]=H\_\{x\}\[N\_\{x\},J\_\{x\}\\zeta\]\\qquad\\text\{for all \}\\zeta\\in T\_\{x\}\\mathcal\{M\}\.
Equivalently,Nx−Jx​ηxN\_\{x\}\-J\_\{x\}\\eta\_\{x\}isHxH\_\{x\}\-orthogonal to every vector inSxS\_\{x\}, whileJx​ηx∈SxJ\_\{x\}\\eta\_\{x\}\\in S\_\{x\}\. The relation is precisely the Hilbert\-space characterization of the orthogonal projection\.

###### Corollary 7\(Quadratic energies\)\.

LetE​\(u\)=12​a​\(u,u\)−ℓ​\(u\)\+cE\(u\)=\\frac\{1\}\{2\}a\(u,u\)\-\\ell\(u\)\+c, wherea:X×X→ℝa:X\\times X\\to\\mathbb\{R\}is symmetric, bounded, and coercive, andℓ∈X∗\\ell\\in X^\{\*\}\. Ifu∗u^\{\*\}is the unique minimizer, equivalentlya​\(u∗,v\)=ℓ​\(v\)a\(u^\{\*\},v\)=\\ell\(v\)for allv∈Xv\\in X, then the undamped EMNGD direction satisfies

Jx​ηx=ΠJx​\(Tx​ℳ\)a​\(P​\(x\)−u∗\)\.J\_\{x\}\\eta\_\{x\}=\\Pi^\{a\}\_\{J\_\{x\}\(T\_\{x\}\\mathcal\{M\}\)\}\(P\(x\)\-u^\{\*\}\)\.\(34\)
The natural\-gradient vector is the projected errorP​\(x\)−u∗P\(x\)\-u^\{\*\}, and the descent updateRx​\(−α​ηx\)R\_\{x\}\(\-\\alpha\\eta\_\{x\}\)moves toward the projected correctionu∗−P​\(x\)u^\{\*\}\-P\(x\)\.

ProofFor a quadratic energy,D​E​\(u\)​\[v\]=a​\(u,v\)−ℓ​\(v\)=a​\(u−u∗,v\)DE\(u\)\[v\]=a\(u,v\)\-\\ell\(v\)=a\(u\-u^\{\*\},v\)\. The Newton vector in theaa\-inner product isu−u∗u\-u^\{\*\}\. The claim follows from Theorem[6](https://arxiv.org/html/2607.22004#Thmtheorem6)withu=P​\(x\)u=P\(x\)andHx=aH\_\{x\}=a\.

###### Proposition 8\(Damping as regularized projection\)\.

Under the assumptions of Theorem[6](https://arxiv.org/html/2607.22004#Thmtheorem6), letNxN\_\{x\}be the function\-space Newton vector\. Forλ\>0\\lambda\>0, the damped EMNGD direction is the unique minimizer of

minξ∈Tx​ℳ⁡\{12​‖Nx−Jx​ξ‖Hx2\+λ2​‖ξ‖02\}\.\\min\_\{\\xi\\in T\_\{x\}\\mathcal\{M\}\}\\left\\\{\\frac\{1\}\{2\}\\\|N\_\{x\}\-J\_\{x\}\\xi\\\|\_\{H\_\{x\}\}^\{2\}\+\\frac\{\\lambda\}\{2\}\\\|\\xi\\\|\_\{0\}^\{2\}\\right\\\}\.\(35\)
Damping turns the exact tangent\-space projection of the Newton vector into a Tikhonov\-regularized tangent\-space projection\.

ProofDifferentiating the objective in \([35](https://arxiv.org/html/2607.22004#S4.E35)\) in the directionζ\\zetagives the stationarity condition

Hx​\[Jx​ξ,Jx​ζ\]\+λ​gx0​\(ξ,ζ\)=Hx​\[Nx,Jx​ζ\]=d​Fx​\[ζ\],H\_\{x\}\[J\_\{x\}\\xi,J\_\{x\}\\zeta\]\+\\lambda g\_\{x\}^\{0\}\(\\xi,\\zeta\)=H\_\{x\}\[N\_\{x\},J\_\{x\}\\zeta\]=dF\_\{x\}\[\\zeta\],which is exactly \([10](https://arxiv.org/html/2607.22004#S3.E10)\)\. Strict convexity follows fromλ\>0\\lambda\>0\.

Equations \([33](https://arxiv.org/html/2607.22004#S4.E33)\) and \([34](https://arxiv.org/html/2607.22004#S4.E34)\) link parameter\-space energy NG to a function\-space Newton update\. For quadratic energies, the function\-space natural\-gradient vector is the current errorP​\(x\)−u∗P\(x\)\-u^\{\*\}\.

###### Proposition 9\(Descent direction\)\.

Assumeλ\>0\\lambda\>0andgrad0⁡F​\(x\)≠0\\operatorname\{grad\}^\{0\}F\(x\)\\neq 0\. Letηx\\eta\_\{x\}be the EMNGD direction\. Then the retraction curveγ​\(α\)=Rx​\(−α​ηx\)\\gamma\(\\alpha\)=R\_\{x\}\(\-\\alpha\\eta\_\{x\}\)satisfies

dd​α​F​\(γ​\(α\)\)\|α=0=−gxE,λ​\(ηx,ηx\)<0\.\\left\.\\frac\{d\}\{d\\alpha\}F\(\\gamma\(\\alpha\)\)\\right\|\_\{\\alpha=0\}=\-g\_\{x\}^\{E,\\lambda\}\(\\eta\_\{x\},\\eta\_\{x\}\)<0\.\(36\)

ProofSinceRRis a retraction,γ′​\(0\)=D​Rx​\(0x\)​\[−ηx\]=−ηx\\gamma^\{\\prime\}\(0\)=DR\_\{x\}\(0\_\{x\}\)\[\-\\eta\_\{x\}\]=\-\\eta\_\{x\}\. The chain rule and \([10](https://arxiv.org/html/2607.22004#S3.E10)\) give

dd​α​F​\(γ​\(α\)\)\|α=0=d​Fx​\[−ηx\]=−gxE,λ​\(ηx,ηx\)\.\\left\.\\frac\{d\}\{d\\alpha\}F\(\\gamma\(\\alpha\)\)\\right\|\_\{\\alpha=0\}=dF\_\{x\}\[\-\\eta\_\{x\}\]=\-g\_\{x\}^\{E,\\lambda\}\(\\eta\_\{x\},\\eta\_\{x\}\)\.
Positive definiteness forλ\>0\\lambda\>0andgrad0⁡F​\(x\)≠0\\operatorname\{grad\}^\{0\}F\(x\)\\neq 0implyηx≠0\\eta\_\{x\}\\neq 0, so the derivative is strictly negative\.

###### Assumption 10\(Uniform metric equivalence and retraction smoothness\)\.

LetΩ=\{x∈ℳ:F​\(x\)≤F​\(x0\)\}\\Omega=\\\{x\\in\\mathcal\{M\}:F\(x\)\\leq F\(x\_\{0\}\)\\\}\. Assume thatFFis bounded below onΩ\\Omegaand that there exist constants0<m≤M<∞0<m\\leq M<\\inftysuch that

m​‖ξ‖02≤gxE,λ​\(ξ,ξ\)≤M​‖ξ‖02\(x∈Ω,ξ∈Tx​ℳ\)\.m\\\|\\xi\\\|\_\{0\}^\{2\}\\leq g\_\{x\}^\{E,\\lambda\}\(\\xi,\\xi\)\\leq M\\\|\\xi\\\|\_\{0\}^\{2\}\\qquad\(x\\in\\Omega,\\ \\xi\\in T\_\{x\}\\mathcal\{M\}\)\.\(37\)
Assume also that there existsLR\>0L\_\{R\}\>0such that every trial step considered by the line search satisfies

F​\(Rx​\(s\)\)≤F​\(x\)\+d​Fx​\[s\]\+LR2​‖s‖02\.F\(R\_\{x\}\(s\)\)\\leq F\(x\)\+dF\_\{x\}\[s\]\+\\frac\{L\_\{R\}\}\{2\}\\\|s\\\|\_\{0\}^\{2\}\.\(38\)

###### Theorem 11\(Global first\-order convergence with Armijo line search\)\.

Suppose Assumptions[1](https://arxiv.org/html/2607.22004#Thmtheorem1)and[10](https://arxiv.org/html/2607.22004#Thmtheorem10)hold withλ\>0\\lambda\>0\. At iterationkk, compute the exact EMNGD directionηk=\(Axkλ\)−1​grad0⁡F​\(xk\)\\eta\_\{k\}=\(A\_\{x\_\{k\}\}^\{\\lambda\}\)^\{\-1\}\\operatorname\{grad\}^\{0\}F\(x\_\{k\}\)\. Chooseαk\\alpha\_\{k\}by backtracking fromα0\>0\\alpha\_\{0\}\>0with contraction factorβ∈\(0,1\)\\beta\\in\(0,1\)until

F​\(Rxk​\(−αk​ηk\)\)≤F​\(xk\)−c​αk​d​Fxk​\[ηk\],F\(R\_\{x\_\{k\}\}\(\-\\alpha\_\{k\}\\eta\_\{k\}\)\)\\leq F\(x\_\{k\}\)\-c\\alpha\_\{k\}\\,dF\_\{x\_\{k\}\}\[\\eta\_\{k\}\],\(39\)holds for somec∈\(0,1\)c\\in\(0,1\)\. Then the line search terminates, the iterates remain inΩ\\Omega, and

‖grad0⁡F​\(xk\)‖0→0\.\\\|\\operatorname\{grad\}^\{0\}F\(x\_\{k\}\)\\\|\_\{0\}\\to 0\.\(40\)
Every accumulation point is a first\-order stationary point\.

ProofLetgk=grad0⁡F​\(xk\)g\_\{k\}=\\operatorname\{grad\}^\{0\}F\(x\_\{k\}\)andAk=AxkλA\_\{k\}=A\_\{x\_\{k\}\}^\{\\lambda\}\. The metric bounds imply

d​Fxk​\[ηk\]=gxkE,λ​\(ηk,ηk\)≥m​‖ηk‖02,d​Fxk​\[ηk\]=gxk0​\(gk,Ak−1​gk\)≥1M​‖gk‖02\.dF\_\{x\_\{k\}\}\[\\eta\_\{k\}\]=g\_\{x\_\{k\}\}^\{E,\\lambda\}\(\\eta\_\{k\},\\eta\_\{k\}\)\\geq m\\\|\\eta\_\{k\}\\\|\_\{0\}^\{2\},\\qquad dF\_\{x\_\{k\}\}\[\\eta\_\{k\}\]=g\_\{x\_\{k\}\}^\{0\}\(g\_\{k\},A\_\{k\}^\{\-1\}g\_\{k\}\)\\geq\\frac\{1\}\{M\}\\\|g\_\{k\}\\\|\_\{0\}^\{2\}\.
Using \([38](https://arxiv.org/html/2607.22004#S4.E38)\) withs=−α​ηks=\-\\alpha\\eta\_\{k\}gives

F​\(Rxk​\(−α​ηk\)\)≤F​\(xk\)−α​\(1−LR​α2​m\)​d​Fxk​\[ηk\]\.F\(R\_\{x\_\{k\}\}\(\-\\alpha\\eta\_\{k\}\)\)\\leq F\(x\_\{k\}\)\-\\alpha\\left\(1\-\\frac\{L\_\{R\}\\alpha\}\{2m\}\\right\)dF\_\{x\_\{k\}\}\[\\eta\_\{k\}\]\.
Every sufficiently smallα\\alphasatisfies \([39](https://arxiv.org/html/2607.22004#S4.E39)\)\. Backtracking terminates and returns a step bounded below by a positive constant\. The accepted steps yieldF​\(xk\+1\)≤F​\(xk\)−C​‖gk‖02F\(x\_\{k\+1\}\)\\leq F\(x\_\{k\}\)\-C\\\|g\_\{k\}\\\|\_\{0\}^\{2\}for someC\>0C\>0independent ofkk\. Summing and using thatFFis bounded below proves∑k‖gk‖02<∞\\sum\_\{k\}\\\|g\_\{k\}\\\|\_\{0\}^\{2\}<\\infty, so \([40](https://arxiv.org/html/2607.22004#S4.E40)\) holds\. Continuity of the Riemannian gradient gives stationarity of any accumulation point\.

###### Proposition 12\(Inexact tangent solves\)\.

Letηx∗=\(Axλ\)−1​grad0⁡F​\(x\)\\eta\_\{x\}^\{\*\}=\(A\_\{x\}^\{\\lambda\}\)^\{\-1\}\\operatorname\{grad\}^\{0\}F\(x\)be the exact direction\. Suppose an approximate directionη~x\\widetilde\{\\eta\}\_\{x\}satisfies

‖η~x−ηx∗‖Axλ≤q​‖ηx∗‖Axλ,q∈\[0,1\),\\\|\\widetilde\{\\eta\}\_\{x\}\-\\eta\_\{x\}^\{\*\}\\\|\_\{A\_\{x\}^\{\\lambda\}\}\\leq q\\\|\\eta\_\{x\}^\{\*\}\\\|\_\{A\_\{x\}^\{\\lambda\}\},\\qquad q\\in\[0,1\),\(41\)where‖ξ‖Axλ2=gx0​\(Axλ​ξ,ξ\)\\\|\\xi\\\|\_\{A\_\{x\}^\{\\lambda\}\}^\{2\}=g\_\{x\}^\{0\}\(A\_\{x\}^\{\\lambda\}\\xi,\\xi\)\. Then

d​Fx​\[η~x\]≥\(1−q\)​‖ηx∗‖Axλ2\>0,dF\_\{x\}\[\\widetilde\{\\eta\}\_\{x\}\]\\geq\(1\-q\)\\\|\\eta\_\{x\}^\{\*\}\\\|\_\{A\_\{x\}^\{\\lambda\}\}^\{2\}\>0,\(42\)whenevergrad0⁡F​\(x\)≠0\\operatorname\{grad\}^\{0\}F\(x\)\\neq 0\. As a result,−η~x\-\\widetilde\{\\eta\}\_\{x\}remains a descent direction\.

ProofWritee=η~x−ηx∗e=\\widetilde\{\\eta\}\_\{x\}\-\\eta\_\{x\}^\{\*\}\. SinceAxλ​ηx∗=grad0⁡F​\(x\)A\_\{x\}^\{\\lambda\}\\eta\_\{x\}^\{\*\}=\\operatorname\{grad\}^\{0\}F\(x\),

d​Fx​\[η~x\]=‖ηx∗‖Axλ2\+⟨ηx∗,e⟩Axλ\.dF\_\{x\}\[\\widetilde\{\\eta\}\_\{x\}\]=\\\|\\eta\_\{x\}^\{\*\}\\\|\_\{A\_\{x\}^\{\\lambda\}\}^\{2\}\+\\langle\\eta\_\{x\}^\{\*\},e\\rangle\_\{A\_\{x\}^\{\\lambda\}\}\.
Cauchy–Schwarz and \([41](https://arxiv.org/html/2607.22004#S4.E41)\) give the lower bound \([42](https://arxiv.org/html/2607.22004#S4.E42)\)\.

##### Computational Complexity and Scalability\.

An EMNGD iteration comprises energy\-operator construction or application, solution of a tangent linear system, and manifold operations\.

We useppfor the ambient parameter dimension,q=dimℳq=\\dim\\mathcal\{M\}for the intrinsic manifold dimension,NNfor the number of residual samples, andq=pq=pfor euclidean parameters\. Letℓ\\ellbe the Nyström rank andmKrylovm\_\{\\mathrm\{Krylov\}\}the number of Krylov iterations\. The costs of one tangent\-operator and one sample\-space kernel application are denoted byCAC\_\{A\}andCKC\_\{K\}, respectively\. The following costs cover additional linear algebra after residual and derivative evaluation\. Residual and derivative costs depend on the PDE operator, network architecture, and automatic\-differentiation implementation\.

Table 1:Dominant linear\-algebra costs of EMNGD solvers\.- •Direct tangent solve\.For residual and generalized Gauss–Newton models, explicit tangent coordinates give a residual JacobianJx∈ℝN×qJ\_\{x\}\\in\\mathbb\{R\}^\{N\\times q\}\. FormingJx⊤​JxJ\_\{x\}^\{\\top\}J\_\{x\}costsO​\(N​q2\)O\(Nq^\{2\}\), and a dense factorization costsO​\(q3\)O\(q^\{3\}\)\. The stated storage excludesJxJ\_\{x\}; retaining the Jacobian addsO​\(N​q\)O\(Nq\)memory\. Direct tangent solves are practical whenqqis moderate\.
- •Exact Woodbury solve\.For a quadratic residual energy or a generalized Gauss–Newton pullback, the Woodbury identity replaces the tangent solve with a sample\-space solve involvingKx=Jx​Jx⊤K\_\{x\}=J\_\{x\}J\_\{x\}^\{\\top\}\. Explicit construction ofKxK\_\{x\}costsO​\(N2​q\)O\(N^\{2\}q\), dense solution costsO​\(N3\)O\(N^\{3\}\), and the back\-projection costsO​\(N​q\)O\(Nq\)\. The route requiresO​\(N2\)O\(N^\{2\}\)additional storage and is favorable whenN≪qN\\ll q\.
- •Matrix\-free and Nyström solvers\.Matrix\-free Krylov methods apply the damped tangent operator without forming a Gram matrix\. The solve costsO​\(mKrylov​CA\)O\(m\_\{\\mathrm\{Krylov\}\}C\_\{A\}\)and requiresO​\(q\)O\(q\)working storage, apart from automatic\-differentiation buffers\. Nyström preconditioning constructs a rank\-ℓ\\ellapproximation to the sample\-space kernel\. Each preconditioned Krylov iteration applies the exact kernel, so iterative convergence recovers the exact Woodbury direction\. A Nyström sketch\-and\-solve method instead returns an approximate direction\.
- •Geometric overhead and operating regimes\.LetCΠC\_\{\\Pi\}andCRC\_\{R\}denote the costs of tangent projection and retraction\. An accepted update addsCΠ\+CRC\_\{\\Pi\}\+C\_\{R\}to the linear\-algebra cost\. Armijo backtracking withnlsn\_\{\\mathrm\{ls\}\}trial steps addsO​\(nls​\(CF\+CR\)\)O\\\!\\left\(n\_\{\\mathrm\{ls\}\}\(C\_\{F\}\+C\_\{R\}\)\\right\), whereCFC\_\{F\}is the energy\-evaluation cost\. Direct tangent solves suit moderateqq, exact Woodbury solves suitN≪qN\\ll q, and Nyström preconditioning reduces Krylov iterations for kernels with useful low\-rank structure\.

## 5Experiments

The experiments address three questions\. First, does EMNGD recover the expected energy\-metric behavior in the Euclidean specialization? Second, do Woodbury duality and Nyström preconditioning compute reliable tangent directions? Third, how do the resulting solvers behave across PDEs, residual counts, and network sizes? The benchmark tables use the Euclidean controlℳ=ℝp\\mathcal\{M\}=\\mathbb\{R\}^\{p\}with the additive retraction\. Separate residual\-formulation diagnostics assess the Jacobian, Gramian, and sample\-space solves\. The diagnostics test implementation consistency rather than architecture\-matched manifold comparisons\.

### 5\.1Experimental Protocol

##### Manifold parametrization\.

For a layer with weight matrixWℓ∈ℝnℓ×nℓ−1W\_\{\\ell\}\\in\\mathbb\{R\}^\{n\_\{\\ell\}\\times n\_\{\\ell\-1\}\}we use the direction–scale decomposition

Wℓ=Diag⁡\(exp⁡\(ρℓ\)\)​Qℓ⊤,Qℓ∈Ob⁡\(nℓ−1,nℓ\),W\_\{\\ell\}=\\operatorname\{Diag\}\\\!\\big\(\\exp\(\\rho\_\{\\ell\}\)\\big\)\\,Q\_\{\\ell\}^\{\\top\},\\qquad Q\_\{\\ell\}\\in\\operatorname\{Ob\}\(n\_\{\\ell\-1\},n\_\{\\ell\}\),\(43\)whereOb⁡\(m,n\)=\{Q∈ℝm×n:diag⁡\(Q⊤​Q\)=𝟏\}\\operatorname\{Ob\}\(m,n\)=\\\{Q\\in\\mathbb\{R\}^\{m\\times n\}:\\operatorname\{diag\}\(Q^\{\\top\}Q\)=\\mathbf\{1\}\\\}is the oblique manifold of unit\-norm columns\. The log\-scalesρℓ∈ℝnℓ\\rho\_\{\\ell\}\\in\\mathbb\{R\}^\{n\_\{\\ell\}\}and biasesbℓ∈ℝnℓb\_\{\\ell\}\\in\\mathbb\{R\}^\{n\_\{\\ell\}\}remain Euclidean\. Every nonzero weight row is a length times a unit direction\. The direction–scale decomposition preserves the network function class and does not reduce the model\. LetΘ=\(\(Qℓ,ρℓ,bℓ\)\)ℓ=1L\\Theta=\(\(Q\_\{\\ell\},\\rho\_\{\\ell\},b\_\{\\ell\}\)\)\_\{\\ell=1\}^\{L\}\. The parameter manifold is

ℳ=∏ℓ=1L\[Ob⁡\(nℓ−1,nℓ\)×ℝnℓ×ℝnℓ\]\.\\mathcal\{M\}=\\prod\_\{\\ell=1\}^\{L\}\\Big\[\\operatorname\{Ob\}\(n\_\{\\ell\-1\},n\_\{\\ell\}\)\\times\\mathbb\{R\}^\{n\_\{\\ell\}\}\\times\\mathbb\{R\}^\{n\_\{\\ell\}\}\\Big\]\.\(44\)
The baseline metric uses the Frobenius metric on oblique factors and the Euclidean metric on scale and bias factors\. For an oblique factor, the tangent space, orthogonal projection, and normalization retraction are

TQ​Ob⁡\(m,n\)\\displaystyle T\_\{Q\}\\operatorname\{Ob\}\(m,n\)=\{Ξ:diag⁡\(Q⊤​Ξ\)=0\},ΠQ​\(Z\)=Z−Q​Diag⁡\(diag⁡\(Q⊤​Z\)\),\\displaystyle=\\\{\\Xi:\\operatorname\{diag\}\(Q^\{\\top\}\\Xi\)=0\\\},\\qquad\\Pi\_\{Q\}\(Z\)=Z\-Q\\operatorname\{Diag\}\\\!\\big\(\\operatorname\{diag\}\(Q^\{\\top\}Z\)\\big\),\(45\)RQ​\(Ξ\)\\displaystyle R\_\{Q\}\(\\Xi\)=\(Q\+Ξ\)Diag\(diag\(\(Q\+Ξ\)⊤\(Q\+Ξ\)\)\)−1/2,\\displaystyle=\(Q\+\\Xi\)\\operatorname\{Diag\}\\\!\\Big\(\\operatorname\{diag\}\\big\(\(Q\+\\Xi\)^\{\\top\}\(Q\+\\Xi\)\\big\)\\Big\)^\{\-1/2\},where the retractionRQR\_\{Q\}renormalizes the columns ofQ\+ΞQ\+\\Xi\. The Euclidean factors useRρ​\(δ​ρ\)=ρ\+δ​ρR\_\{\\rho\}\(\\delta\\rho\)=\\rho\+\\delta\\rhoandRb​\(δ​b\)=b\+δ​bR\_\{b\}\(\\delta b\)=b\+\\delta b\. The unconstrained Euclidean control drops the oblique constraint, soℳ=ℝp\\mathcal\{M\}=\\mathbb\{R\}^\{p\}andRΘ​\(v\)=Θ\+vR\_\{\\Theta\}\(v\)=\\Theta\+v\.

The PDE benchmarks compare optimizers in the Euclidean control\. The product\-manifold parametrization defines the constrained EMNGD setting\. The residual diagnostics test the tangent\-space and sample\-space computations separately\.

##### Intrinsic EMNGD direction\.

Letr​\(Θ\)r\(\\Theta\)collect the weighted interior, boundary, and initial residuals\. The tangent differential𝒥Θ:TΘ​ℳ→ℝN\\mathcal\{J\}\_\{\\Theta\}:T\_\{\\Theta\}\\mathcal\{M\}\\to\\mathbb\{R\}^\{N\}maps a tangent direction to the residual change\. At iterationkk, the damped EMNGD direction solves

\(𝒥Θk∗​𝒥Θk\+λk​I\)​ηk=𝒥Θk∗​r​\(Θk\)in​TΘk​ℳ,\\big\(\\mathcal\{J\}\_\{\\Theta\_\{k\}\}^\{\*\}\\mathcal\{J\}\_\{\\Theta\_\{k\}\}\+\\lambda\_\{k\}I\\big\)\\eta\_\{k\}=\\mathcal\{J\}\_\{\\Theta\_\{k\}\}^\{\*\}\\,r\(\\Theta\_\{k\}\)\\qquad\\text\{in \}T\_\{\\Theta\_\{k\}\}\\mathcal\{M\},\(46\)The Woodbury form \([27](https://arxiv.org/html/2607.22004#S3.E27)\) gives the same tangent direction with kernelKΘk=𝐉k​ΠΘk​𝐉k⊤K\_\{\\Theta\_\{k\}\}=\\mathbf\{J\}\_\{k\}\\Pi\_\{\\Theta\_\{k\}\}\\mathbf\{J\}\_\{k\}^\{\\top\}\. The line search evaluatesΘk​\(α\)=RΘk​\(−α​ηk\)\\Theta\_\{k\}\(\\alpha\)=R\_\{\\Theta\_\{k\}\}\(\-\\alpha\\eta\_\{k\}\)\.

We globalize EMNGD with the retraction\-based Armijo line search in Algorithm[1](https://arxiv.org/html/2607.22004#alg1)\. The search starts from the full trial stepα=1\\alpha=1\. The projected\-Newton interpretation of the undamped direction motivates that initial value\. Nonlinear realization maps, parameter manifolds, and retractions require a sufficient\-decrease test before accepting the full trial step\.

The experiments evaluate candidate step sizes on the geometric grid

𝒜=\{1,β,β2,…,βmls\}⊂\(0,1\],β∈\(0,1\)\.\\mathcal\{A\}=\\\{1,\\beta,\\beta^\{2\},\\ldots,\\beta^\{m\_\{\\rm ls\}\}\\\}\\subset\(0,1\],\\qquad\\beta\\in\(0,1\)\.
Energy evaluations on the grid can run in parallel\. The implementation selects the largest candidate satisfying the Armijo condition\. Further backtracking extends the grid when no candidate is accepted\.

Algorithm 1Damped EMNGD with retraction\-based Armijo line search1:Input:initial point

x0∈ℳx\_\{0\}\\in\\mathcal\{M\}; positive damping parameters

\{λk\}k≥0\\\{\\lambda\_\{k\}\\\}\_\{k\\geq 0\}; Armijo parameter

c∈\(0,1\)c\\in\(0,1\); backtracking factor

β∈\(0,1\)\\beta\\in\(0,1\); gradient tolerance

εgrad\>0\\varepsilon\_\{\\rm grad\}\>0; maximum iterations

NmaxN\_\{\\max\}
2:for

k=0,…,Nmax−1k=0,\\ldots,N\_\{\\max\}\-1do

3:Compute

gk=grad0⁡F​\(xk\)∈Txk​ℳg\_\{k\}=\\operatorname\{grad\}^\{0\}F\(x\_\{k\}\)\\in T\_\{x\_\{k\}\}\\mathcal\{M\}\.

4:if

‖gk‖gxk0≤εgrad\\\|g\_\{k\}\\\|\_\{g^\{0\}\_\{x\_\{k\}\}\}\\leq\\varepsilon\_\{\\rm grad\}then

5:stop

6:endif

7:Define the

gxk0g^\{0\}\_\{x\_\{k\}\}\-self\-adjoint tangent operator

Akλk:Txk​ℳ→Txk​ℳA\_\{k\}^\{\\lambda\_\{k\}\}:T\_\{x\_\{k\}\}\\mathcal\{M\}\\to T\_\{x\_\{k\}\}\\mathcal\{M\}by

8:

gxk0​\(Akλk​ξ,ζ\)=gxkE​\(ξ,ζ\)\+λk​gxk0​\(ξ,ζ\)g^\{0\}\_\{x\_\{k\}\}\(A\_\{k\}^\{\\lambda\_\{k\}\}\\xi,\\zeta\)=g^\{E\}\_\{x\_\{k\}\}\(\\xi,\\zeta\)\+\\lambda\_\{k\}g^\{0\}\_\{x\_\{k\}\}\(\\xi,\\zeta\)for all

ξ,ζ∈Txk​ℳ\\xi,\\zeta\\in T\_\{x\_\{k\}\}\\mathcal\{M\}\.

9:Solve

Akλk​ηk=gkA\_\{k\}^\{\\lambda\_\{k\}\}\\eta\_\{k\}=g\_\{k\}in

Txk​ℳT\_\{x\_\{k\}\}\\mathcal\{M\}exactly or to a prescribed inner\-solver tolerance\.

10:if

d​Fxk​\[ηk\]≤0dF\_\{x\_\{k\}\}\[\\eta\_\{k\}\]\\leq 0then

11:Increase

λk\\lambda\_\{k\}, or tighten the inner\-solver tolerance, and recompute

ηk\\eta\_\{k\}\.

12:endif

13:Set

αk←1\\alpha\_\{k\}\\leftarrow 1\.

14:while

F​\(Rxk​\(−αk​ηk\)\)\>F​\(xk\)−c​αk​d​Fxk​\[ηk\]F\\\!\\left\(R\_\{x\_\{k\}\}\(\-\\alpha\_\{k\}\\eta\_\{k\}\)\\right\)\>F\(x\_\{k\}\)\-c\\alpha\_\{k\}dF\_\{x\_\{k\}\}\[\\eta\_\{k\}\]do

15:

αk←β​αk\\alpha\_\{k\}\\leftarrow\\beta\\alpha\_\{k\}\.

16:endwhile

17:Update

xk\+1=Rxk​\(−αk​ηk\)x\_\{k\+1\}=R\_\{x\_\{k\}\}\(\-\\alpha\_\{k\}\\eta\_\{k\}\)\.

18:endfor

For quadratic residual energies or generalized Gauss–Newton pullback metrics, Algorithm[2](https://arxiv.org/html/2607.22004#alg2)uses an embedded product manifold with the metric induced by the ambient Euclidean product space\. The induced metric gives𝒥Θ∗=ΠΘ​𝐉Θ⊤\\mathcal\{J\}\_\{\\Theta\}^\{\*\}=\\Pi\_\{\\Theta\}\\mathbf\{J\}\_\{\\Theta\}^\{\\top\}\. General Riemannian metrics require metric\-dependent projectors and adjoints\.

Algorithm 2EMNGD with exact Woodbury and Nyström\-preconditioned solves1:Input:initial point

Θ0∈ℳ\\Theta\_\{0\}\\in\\mathcal\{M\}; positive damping parameters

\{λk\}k≥0\\\{\\lambda\_\{k\}\\\}\_\{k\\geq 0\}; solver mode

𝗆𝗈𝖽𝖾∈\{𝖣𝗂𝗋𝖾𝖼𝗍,𝖭𝗒𝗌𝖯𝖢𝖦\}\\mathsf\{mode\}\\in\\\{\\mathsf\{Direct\},\\mathsf\{NysPCG\}\\\}; Nyström rank

ℓ\\ell; linear\-solver tolerance

εlin\\varepsilon\_\{\\rm lin\}; Armijo parameters

c,β∈\(0,1\)c,\\beta\\in\(0,1\); maximum iterations

NmaxN\_\{\\max\}
2:for

k=0,…,Nmax−1k=0,\\ldots,N\_\{\\max\}\-1do

3:Evaluate

rk=r​\(Θk\)∈ℝNr\_\{k\}=r\(\\Theta\_\{k\}\)\\in\\mathbb\{R\}^\{N\}\.

4:Define the ambient residual Jacobian

𝐉k=D​r​\(Θk\):ℝp→ℝN\\mathbf\{J\}\_\{k\}=Dr\(\\Theta\_\{k\}\):\\mathbb\{R\}^\{p\}\\to\\mathbb\{R\}^\{N\}through explicit assembly or matrix\-free Jacobian products\.

5:Form or apply the

g0g^\{0\}\-orthogonal tangent projector

Πk=ΠΘk:ℝp→TΘk​ℳ\\Pi\_\{k\}=\\Pi\_\{\\Theta\_\{k\}\}:\\mathbb\{R\}^\{p\}\\to T\_\{\\Theta\_\{k\}\}\\mathcal\{M\}\.

6:Define

𝒥k=𝐉k\|TΘk​ℳ\\mathcal\{J\}\_\{k\}=\\left\.\\mathbf\{J\}\_\{k\}\\right\|\_\{T\_\{\\Theta\_\{k\}\}\\mathcal\{M\}\}and

𝒥k∗=Πk​𝐉k⊤\\mathcal\{J\}\_\{k\}^\{\*\}=\\Pi\_\{k\}\\mathbf\{J\}\_\{k\}^\{\\top\}\.

7:Define

Kk=𝒥k​𝒥k∗=𝐉k​Πk​𝐉k⊤K\_\{k\}=\\mathcal\{J\}\_\{k\}\\mathcal\{J\}\_\{k\}^\{\*\}=\\mathbf\{J\}\_\{k\}\\Pi\_\{k\}\\mathbf\{J\}\_\{k\}^\{\\top\}\.

8:if

𝗆𝗈𝖽𝖾=𝖣𝗂𝗋𝖾𝖼𝗍\\mathsf\{mode\}=\\mathsf\{Direct\}then

9:Solve

\(Kk\+λk​IN\)​ak=rk\(K\_\{k\}\+\\lambda\_\{k\}I\_\{N\}\)a\_\{k\}=r\_\{k\}by a direct symmetric positive\-definite solver\.

10:else

11:Construct a rank\-

ℓ\\ellNyström approximation

K^k\\widehat\{K\}\_\{k\}of

KkK\_\{k\}and set

Mk=K^k\+λk​INM\_\{k\}=\\widehat\{K\}\_\{k\}\+\\lambda\_\{k\}I\_\{N\}\.

12:Solve the*exact*system

\(Kk\+λk​IN\)​ak=rk\(K\_\{k\}\+\\lambda\_\{k\}I\_\{N\}\)a\_\{k\}=r\_\{k\}by preconditioned conjugate gradients with

MkM\_\{k\}until

‖\(Kk\+λk​IN\)​ak−rk‖2/‖rk‖2≤εlin\\\|\(K\_\{k\}\+\\lambda\_\{k\}I\_\{N\}\)a\_\{k\}\-r\_\{k\}\\\|\_\{2\}/\\\|r\_\{k\}\\\|\_\{2\}\\leq\\varepsilon\_\{\\rm lin\}\.

13:endif

14:Reconstruct

ηk=𝒥k∗​ak=Πk​𝐉k⊤​ak∈TΘk​ℳ\\eta\_\{k\}=\\mathcal\{J\}\_\{k\}^\{\*\}a\_\{k\}=\\Pi\_\{k\}\\mathbf\{J\}\_\{k\}^\{\\top\}a\_\{k\}\\in T\_\{\\Theta\_\{k\}\}\\mathcal\{M\}\.

15:if

d​FΘk​\[ηk\]≤0dF\_\{\\Theta\_\{k\}\}\[\\eta\_\{k\}\]\\leq 0then

16:Increase

λk\\lambda\_\{k\}, or reduce

εlin\\varepsilon\_\{\\rm lin\}, and recompute

aka\_\{k\}and

ηk\\eta\_\{k\}\.

17:endif

18:Set

αk←1\\alpha\_\{k\}\\leftarrow 1and apply the Armijo backtracking rule from Algorithm[1](https://arxiv.org/html/2607.22004#alg1)\.

19:Update

Θk\+1=RΘk​\(−αk​ηk\)\\Theta\_\{k\+1\}=R\_\{\\Theta\_\{k\}\}\(\-\\alpha\_\{k\}\\eta\_\{k\}\)\.

20:endfor

Algorithm[1](https://arxiv.org/html/2607.22004#alg1)is the intrinsic, coordinate\-free EMNGD method\. The algorithm defines the energy\-metric equation on the current tangent space and uses retraction\-based Armijo backtracking to obtain the next feasible iterate\.

Algorithm[2](https://arxiv.org/html/2607.22004#alg2)specializes EMNGD to quadratic residual energies or generalized Gauss–Newton pullback metrics on embedded product manifolds with the induced Euclidean product metric\. A direct sample\-space solve computes the exact damped EMNGD direction through the Woodbury identity\. Nyström\-preconditioned Krylov iteration solves the same sample\-space system to a prescribed tolerance\. Nyström changes the conditioning of the inner solve but does not change the target EMNGD direction\.

The unconstrained Euclidean specialization follows fromℳ=ℝp\\mathcal\{M\}=\\mathbb\{R\}^\{p\},ΠΘ=Ip\\Pi\_\{\\Theta\}=I\_\{p\}, andRΘ​\(η\)=Θ\+ηR\_\{\\Theta\}\(\\eta\)=\\Theta\+\\eta\. Section[4](https://arxiv.org/html/2607.22004#S4.SS0.SSS0.Px1)compares parameter\-space, sample\-space, and matrix\-free costs\.

In local coordinates atΘk\\Theta\_\{k\}, letGE,kϕG\_\{E,k\}^\{\\phi\},G0,kϕG\_\{0,k\}^\{\\phi\}, andbkb\_\{k\}denote the energy Gramian, baseline Gramian, and coordinate representation ofd​FΘkdF\_\{\\Theta\_\{k\}\}\. The damped EMNGD direction is the unique minimizer of the strictly convex tangent quadratic model

vk=arg⁡minv∈ℝq⁡\{12​v⊤​\(GE,kϕ\+λk​G0,kϕ\)​v−bk⊤​v\},q=dimℳ\.v\_\{k\}=\\arg\\min\_\{v\\in\\mathbb\{R\}^\{q\}\}\\left\\\{\\frac\{1\}\{2\}v^\{\\top\}\\big\(G\_\{E,k\}^\{\\phi\}\+\\lambda\_\{k\}G\_\{0,k\}^\{\\phi\}\\big\)v\-b\_\{k\}^\{\\top\}v\\right\\\},\\qquad q=\\dim\\mathcal\{M\}\.\(47\)
For quadratic residual energies or the generalized Gauss–Newton metric, the same direction solves the Tikhonov\-regularized tangent least\-squares problem

vk=arg⁡minv∈ℝq⁡\{12​‖Jϕ,k​v−rk‖22\+λk2​v⊤​G0,kϕ​v\},v\_\{k\}=\\arg\\min\_\{v\\in\\mathbb\{R\}^\{q\}\}\\left\\\{\\frac\{1\}\{2\}\\\|J\_\{\\phi,k\}v\-r\_\{k\}\\\|\_\{2\}^\{2\}\+\\frac\{\\lambda\_\{k\}\}\{2\}v^\{\\top\}G\_\{0,k\}^\{\\phi\}v\\right\\\},\(48\)whereJϕ,kJ\_\{\\phi,k\}is the residual Jacobian in the selected tangent basis\.

When the undamped tangent system is singular, the Moore–Penrose direction requires an explicit minimum\-g0g^\{0\}\-norm solution convention\. The computational algorithms useλk\>0\\lambda\_\{k\}\>0to ensure uniqueness and improve conditioning\.

##### Initialization, sampling, and budgets\.

For loss and Gram\-matrix integrals, we use fixed regular grids or resampled random points\. We initialize weights and biases from a zero\-mean Gaussian with standard deviation0\.10\.1\. On the product manifold, eachQℓQ\_\{\\ell\}is the columnwise normalization of the corresponding Gaussian matrix\. The log\-scalesρℓ\\rho\_\{\\ell\}and biases remain Euclidean\. The Euclidean control uses the unnormalized Gaussian initialization\. Each PDE subsection states the collocation rule and iteration budget\. Tables list runtime separately from iteration counts\. The studies are not wall\-clock\-matched comparisons\.

#### 5\.1\.1Evaluation Metrics and Baselines

We report relativeL2L^\{2\}error and, where derivative evaluations are available, relativeH1H^\{1\}error\. Evaluation uses denser quadrature than optimization\. The Euclidean\-control studies compare stochastic gradient descent \(SGD\), Adam, BFGS\(Nocedal and Wright,[1999](https://arxiv.org/html/2607.22004#bib.bib1)\), and ENGD\. SGD uses a logarithmic line\-search grid\. Adam starts at10−310^\{\-3\}\. After1\.5×1041\.5\\times 10^\{4\}steps, the learning rate decreases by a factor of10−110^\{\-1\}every10410^\{4\}steps\. The schedule stops at10−710^\{\-7\}or the iteration budget\.

The convergence and mechanism figures display NGD, Hessian\-free, Woodbury, SPRING, and Nyström variants when the corresponding curve is labelled\. The displayed trajectories use the stated configurations and do not replace the 10\-initialization table protocol\. The task\-specific iteration budgets state the computational allocations\.

The preliminary one\-dimensional studies establish the Euclidean reduction and the solver identity\. Figure[4](https://arxiv.org/html/2607.22004#S5.F4)summarizes loss and final relativeL2L^\{2\}error across the displayed PDEs\. The labelled EMNGD curve reaches the lowest displayed final errors\. Figure[5](https://arxiv.org/html/2607.22004#S5.F5)shows that Woodbury ENGD follows the parameter\-space ENGD trajectory\. The sample\-space solve therefore changes the linear algebra, not the direction\.

Figure[6](https://arxiv.org/html/2607.22004#S5.F6)gives a spatial check on one\-dimensional Poisson\. EMNGD overlaps the reference solution and keeps the pointwise error near10−710^\{\-7\}or lower over most of the domain\. ENGD and NGD also track the reference, but retain interior errors near10−510^\{\-5\}\. The diagnostic ordering supports the trajectory results but does not replace matched\-budget comparisons\.

![Refer to caption](https://arxiv.org/html/2607.22004v1/x3.png)Figure 4:Training loss \(top\) and final relativeL2L^\{2\}error \(bottom\) across PDE benchmarks\.![Refer to caption](https://arxiv.org/html/2607.22004v1/x4.png)Figure 5:One\-dimensional PDE benchmark\.![Refer to caption](https://arxiv.org/html/2607.22004v1/x5.png)Figure 6:One\-dimensional Poisson solutions and pointwise errors \(log scale\)\.
#### 5\.1\.2Implementation Details

We implement the solvers in JAX\(Bradburyet al\.,[2018](https://arxiv.org/html/2607.22004#bib.bib2)\)with automatic differentiation\. Least\-squares solves use singular\-value decomposition\. BFGS usesjaxopt\.BFGS\. Unless stated otherwise, experiments run in double precision on one NVIDIA RTX 5090 Laptop GPU\. The implementation is available at[https://github\.com/liangzhangyong/EMNGD](https://github.com/liangzhangyong/EMNGD)\.

### 5\.2Residual\-Formulation Diagnostics

The residual\-formulation diagnostic verifies the hard\-Dirichlet embedding and Woodbury solve before PDE accuracy comparisons\. The trial mapuθ​\(x\)=∏i=12xi​\(1−xi\)​vθ​\(x\)u\_\{\\theta\}\(x\)=\\prod\_\{i=1\}^\{2\}x\_\{i\}\(1\-x\_\{i\}\)v\_\{\\theta\}\(x\)imposes the boundary condition by construction, leavingu∗u^\{\\ast\}outside the training objective\. For 48 fixed interior residual points and a 337\-parameter network, 80 updates reduce the residual loss by more than 14 orders of magnitude, from5\.521×1015\.521\\times 10^\{1\}to2\.941×10−132\.941\\times 10^\{\-13\}, and yield a held\-out relativeL2L^\{2\}error of2\.208×10−42\.208\\times 10^\{\-4\}\. The diagnostic supports the correct interaction of the constraint embedding and Woodbury solver on the stated fixed\-sample problem\.

The primal and Woodbury directions agree to relative error8\.91×10−98\.91\\times 10^\{\-9\}, verifying the dual implementation at numerical precision\. Figures[7](https://arxiv.org/html/2607.22004#S5.F7)and[10](https://arxiv.org/html/2607.22004#S5.F10)reveal evolving residual\-Fisher geometry and a wide spectral range, which motivates damping in the sample\-space solve\. Figure[8](https://arxiv.org/html/2607.22004#S5.F8)displays the48×4848\\times 48kernel system that replaces the337×337337\\times 337parameter\-space system, while Figure[9](https://arxiv.org/html/2607.22004#S5.F9)reconstructs the layerwise Gramian to error4\.36×10−164\.36\\times 10^\{\-16\}\. In contrast, Figure[12](https://arxiv.org/html/2607.22004#S5.F12)has relative Frobenius error7\.08×10−17\.08\\times 10^\{\-1\}, so weight sharing changes the kernel and remains an approximation\.

![Refer to caption](https://arxiv.org/html/2607.22004v1/x6.png)Figure 7:Residual\-Fisher geometry for hard\-Dirichlet EMNGD\.![Refer to caption](https://arxiv.org/html/2607.22004v1/x7.png)Figure 8:Residual Jacobians for Woodbury EMNGD on two\-dimensional Poisson\.![Refer to caption](https://arxiv.org/html/2607.22004v1/x8.png)Figure 9:Layerwise contributions to the residual Gramian\.![Refer to caption](https://arxiv.org/html/2607.22004v1/x9.png)Figure 10:Sample\-space residual\-Gramian spectrum for128128residual samples\.Figure[11](https://arxiv.org/html/2607.22004#S5.F11)compares the exact Woodbury solve with rank\-900 Nyström preconditioning\. Both reduce residual loss to the10−1210^\{\-12\}scale and reach relativeL2L^\{2\}errors near10−810^\{\-8\}\. Woodbury finishes slightly lower \(1\.09×10−81\.09\\times 10^\{\-8\}versus1\.45×10−81\.45\\times 10^\{\-8\}\)\. The Nyström trajectory closely follows Woodbury\.

![Refer to caption](https://arxiv.org/html/2607.22004v1/x10.png)Figure 11:EMNGD with an exact Woodbury solve and rank\-900 Nyström preconditioning\. Left: residual loss\. Right: relativeL2L^\{2\}error\.![Refer to caption](https://arxiv.org/html/2607.22004v1/x11.png)Figure 12:Weight\-sharing approximation for an EMNGD residual\-Jacobian block\.
### 5\.3Poisson Equation

We consider the two\-dimensional Poisson equation

−Δ​u​\(x,y\)=f​\(x,y\)=2​π2​sin⁡\(π​x\)​sin⁡\(π​y\),\-\\Delta u\(x,y\)=f\(x,y\)=2\\pi^\{2\}\\sin\(\\pi x\)\\sin\(\\pi y\),on the unit square\[0,1\]2\[0,1\]^\{2\}with zero boundary values\. The solution is given by

u∗​\(x,y\)=sin⁡\(π​x\)​sin⁡\(π​y\),u^\{\*\}\(x,y\)=\\sin\(\\pi x\)\\sin\(\\pi y\),and the PINNs loss of the problem is

L​\(θ\)=1NΩ​∑i=1NΩ\(Δ​uθ​\(xi,yi\)\+f​\(xi,yi\)\)2\+1N∂Ω​∑i=1N∂Ωuθ​\(xib,yib\)2,\\displaystyle\\begin\{split\}L\(\\theta\)&=\\frac\{1\}\{N\_\{\\Omega\}\}\\sum\_\{i=1\}^\{N\_\{\\Omega\}\}\(\\Delta u\_\{\\theta\}\(x\_\{i\},y\_\{i\}\)\+f\(x\_\{i\},y\_\{i\}\)\)^\{2\}\\\\ &\\qquad\\qquad\\quad\+\\frac\{1\}\{N\_\{\\partial\\Omega\}\}\\sum\_\{i=1\}^\{N\_\{\\partial\\Omega\}\}u\_\{\\theta\}\(x^\{b\}\_\{i\},y^\{b\}\_\{i\}\)^\{2\},\\end\{split\}\(49\)where\{\(xi,yi\)\}i=1,…,NΩ\\\{\(x\_\{i\},y\_\{i\}\)\\\}\_\{i=1,\\dots,N\_\{\\Omega\}\}denote the interior collocation points and\{\(xib,yib\)\}i=1,…,N∂Ω\\\{\(x^\{b\}\_\{i\},y^\{b\}\_\{i\}\)\\\}\_\{i=1,\\dots,N\_\{\\partial\\Omega\}\}denote the collocation points on∂Ω\\partial\\Omega\. For the Poisson problem, the energy inner product onH2​\(Ω\)H^\{2\}\(\\Omega\)is

a​\(u,v\)=∫ΩΔ​u​Δ​v​dx\+∫∂Ωu​v​ds\.a\(u,v\)=\\int\_\{\\Omega\}\\Delta u\\Delta v\\mathrm\{d\}x\+\\int\_\{\\partial\\Omega\}uv\\mathrm\{d\}s\.\(50\)
The energy inner product is not coercive333The inner product is coercive with respect to theH1/2​\(Ω\)H^\{1/2\}\(\\Omega\)norm; see\(Müller and Zeinhofer,[2022b](https://arxiv.org/html/2607.22004#bib.bib339)\)\.onH2​\(Ω\)H^\{2\}\(\\Omega\)and differs from theH2​\(Ω\)H^\{2\}\(\\Omega\)inner product\. We approximate \([50](https://arxiv.org/html/2607.22004#S5.E50)\) with the collocation points from \([49](https://arxiv.org/html/2607.22004#S5.E49)\)\. The reproducibility rerun uses a common22–3232–11network\. SGD and Adam run for the recorded long\-horizon updates\. BFGS, ENGD, and SPRING run for 50 updates\. EMNGD runs for 20 updates\.

Table 2:Single\-seed Poisson2D reproducibility reruns\.Table 3:Poisson2D runtime records under the listed settings\.Table 4:Imported Poisson2D baseline integration results\.Table 5:Poisson2D endpoints of native solver implementations forD=8,577D=8\{,\}577\.The two\-dimensional study separates reproducibility, implementation coverage, and parameter scaling\. Table[2](https://arxiv.org/html/2607.22004#S5.T2)uses a common22–3232–11network\. EMNGD reaches a relativeL2L^\{2\}error of6\.778×10−96\.778\\times 10^\{\-9\}after 20 updates\. The external runners in Table[4](https://arxiv.org/html/2607.22004#S5.T4)use 257 parameters, whereas the native EMNGD run uses 8,577 parameters\. Table[5](https://arxiv.org/html/2607.22004#S5.T5)therefore records endpoint coverage rather than an architecture\-matched ranking\.

Figure[13](https://arxiv.org/html/2607.22004#S5.F13)tests parameter scaling atD=8,577D=8\{,\}577,9,8739\{,\}873, and116,097116\{,\}097\. The labelled EMNGD curve reaches relativeL2L^\{2\}errors near10−710^\{\-7\}within10310^\{3\}iterations in all three settings\. Table[6](https://arxiv.org/html/2607.22004#S5.T6)gives the corresponding terminal values\. Unequal stopping rules prevent a matched wall\-clock or iteration\-budget ranking\.

![Refer to caption](https://arxiv.org/html/2607.22004v1/x12.png)Figure 13:Two\-dimensional Poisson relativeL2L^\{2\}error across three parameter dimensions\.Table 6:Terminal losses and relativeL2L^\{2\}errors for two\-dimensional Poisson\.Table[6](https://arxiv.org/html/2607.22004#S5.T6)lists the terminal values in Figure[13](https://arxiv.org/html/2607.22004#S5.F13)\. The entries are archived trajectory endpoints and recorded hard residual\-manifold EMNGD runs withoutu∗u^\{\\ast\}in training\. TheD=116,097D=116\{,\}097ENGD entry uses a diagonal approximation\. The endpoints do not form a matched\-budget ranking\.

### 5\.4Five\-Dimensional Poisson Equation

We next consider the Poisson equation in five spatial dimensions:

−Δ​u\\displaystyle\-\\Delta u=fin​\[0,1\]5,\\displaystyle=f\\quad\\quad\\ \\ \\ \\quad\\quad\\quad\\text\{in \}\[0,1\]^\{5\},u​\(x\)\\displaystyle u\(x\)=∑k=15sin⁡\(π​xk\)on​∂\[0,1\]5\.\\displaystyle=\\sum\_\{k=1\}^\{5\}\\sin\(\\pi x\_\{k\}\)\\quad\\ \\text\{on \}\\partial\[0,1\]^\{5\}\.
We use the manufactured solution

u∗:ℝ5→ℝ,x↦∑k=15sin⁡\(π​xk\)u^\{\\ast\}\\colon\\mathbb\{R\}^\{5\}\\to\\mathbb\{R\},\\quad x\\mapsto\\sum\_\{k=1\}^\{5\}\\sin\(\\pi x\_\{k\}\)sof=π2​u∗f=\\pi^\{2\}u^\{\\ast\}\. We use the loss and energy inner product from Equations[49](https://arxiv.org/html/2607.22004#S5.E49)and[50](https://arxiv.org/html/2607.22004#S5.E50)\. Each optimization step drawsNΩ=3000N\_\{\\Omega\}=3000interior points andN∂Ω=500N\_\{\\partial\\Omega\}=500boundary points\. The Euclidean\-control network has five inputs, 64 hyperbolic\-tangent hidden units, and one output\.

The five\-dimensional problem tests whether the sample\-space solvers retain the energy\-metric advantage as residual evaluation becomes more expensive\. Figure[14](https://arxiv.org/html/2607.22004#S5.F14)reports relativeL2L^\{2\}error against iterations and wall\-clock time\. Under the Euclidean\-control protocol, energy\-metric curves reach lower errors than the first\-order baselines\. The comparison tests the Euclidean reduction, not a non\-Euclidean manifold advantage\.

![Refer to caption](https://arxiv.org/html/2607.22004v1/x13.png)Figure 14:Five\-dimensional Poisson: relativeL2L^\{2\}error versus iteration and time\.Figure[15](https://arxiv.org/html/2607.22004#S5.F15)separates training loss from evaluation error\. Energy\-based curves reduce both quantities more rapidly than the first\-order curves\. The displayed EMNGD trajectory reaches a low\-error regime in both iteration and time\.

![Refer to caption](https://arxiv.org/html/2607.22004v1/x14.png)Figure 15:Loss and relativeL2L^\{2\}error for five\-dimensional Poisson\.Figure[16](https://arxiv.org/html/2607.22004#S5.F16)then increases the number of collocation residuals fromN=1000N=1000toN=10000N=10000\. Woodbury and Nyström variants reduce loss and relativeL2L^\{2\}error across all three sample sizes\. Nyström follows the Woodbury convergence pattern while avoiding a dense sample\-space solve\. The randomized SPRING curves vary more when the sampled kernel is less stable\.

![Refer to caption](https://arxiv.org/html/2607.22004v1/x15.png)Figure 16:Large\-sample five\-dimensional Poisson benchmark\.Figure[17](https://arxiv.org/html/2607.22004#S5.F17)gives a metric\-focused view of the same benchmark\. The energy\-metric curves have lower loss and relativeL2L^\{2\}error than the first\-order curves\. Hessian\-free and SPRING improve on the first\-order baselines, while the labelled EMNGD curve reaches the lowest displayed error range\.

![Refer to caption](https://arxiv.org/html/2607.22004v1/x16.png)Figure 17:Metric comparison for five\-dimensional Poisson\.
### 5\.5Heat Equation

Let us consider the one\-dimensional heat equation

∂tu​\(t,x\)\\displaystyle\\partial\_\{t\}u\(t,x\)=14​∂x2u​\(t,x\)for​\(t,x\)∈\[0,1\]2\\displaystyle=\\frac\{1\}\{4\}\\partial\_\{x\}^\{2\}u\(t,x\)\\quad\\text\{for \}\(t,x\)\\in\[0,1\]^\{2\}u​\(0,x\)\\displaystyle u\(0,x\)=sin⁡\(π​x\)for​x∈\[0,1\]\\displaystyle=\\sin\(\\pi x\)\\qquad\\;\\,\\text\{for \}x\\in\[0,1\]u​\(t,x\)\\displaystyle u\(t,x\)=0for​\(t,x\)∈\[0,1\]×\{0,1\}\.\\displaystyle=0\\qquad\\qquad\\quad\\text\{for \}\(t,x\)\\in\[0,1\]\\times\\\{0,1\\\}\.
The solution is given by

u∗​\(t,x\)=exp⁡\(−π2​t4\)​sin⁡\(π​x\),u^\{\*\}\(t,x\)=\\exp\\left\(\-\\frac\{\\pi^\{2\}t\}\{4\}\\right\)\\sin\(\\pi x\),and the PINNs loss is

L​\(θ\)\\displaystyle L\(\\theta\)=1NΩT​∑i=1NΩT\(∂tuθ​\(ti,xi\)−14​∂x2uθ​\(ti,xi\)\)2\\displaystyle=\\frac\{1\}\{N\_\{\\Omega\_\{T\}\}\}\\sum\_\{i=1\}^\{N\_\{\\Omega\_\{T\}\}\}\\left\(\\partial\_\{t\}u\_\{\\theta\}\(t\_\{i\},x\_\{i\}\)\-\\frac\{1\}\{4\}\\partial\_\{x\}^\{2\}u\_\{\\theta\}\(t\_\{i\},x\_\{i\}\)\\right\)^\{2\}\+1Nin​∑i=1NΩ\(uθ​\(0,xiin\)−sin⁡\(π​xiin\)\)2\\displaystyle\\quad\+\\frac\{1\}\{N\_\{\\text\{in\}\}\}\\sum\_\{i=1\}^\{N\_\{\\Omega\}\}\\left\(u\_\{\\theta\}\(0,x\_\{i\}^\{\\text\{in\}\}\)\-\\sin\(\\pi x\_\{i\}^\{\\text\{in\}\}\)\\right\)^\{2\}\+1N∂Ω​∑i=1N∂Ωuθ​\(tib,xib\)2,\\displaystyle\\quad\+\\frac\{1\}\{N\_\{\\partial\\Omega\}\}\\sum\_\{i=1\}^\{N\_\{\\partial\\Omega\}\}u\_\{\\theta\}\(t^\{b\}\_\{i\},x^\{b\}\_\{i\}\)^\{2\},where\{\(ti,xi\)\}i=1,…,NΩT\\\{\(t\_\{i\},x\_\{i\}\)\\\}\_\{i=1,\\dots,N\_\{\\Omega\_\{T\}\}\}are interior space\-time collocation points\. The set\{\(tib,xib\)\}i=1,…,N∂Ω\\\{\(t\_\{i\}^\{b\},x\_\{i\}^\{b\}\)\\\}\_\{i=1,\\dots,N\_\{\\partial\\Omega\}\}contains spatial\-boundary points\. The set\{\(xiin\)\}i=1,…,Nin\\\{\(x\_\{i\}^\{\\text\{in\}\}\)\\\}\_\{i=1,\\dots,N\_\{\\text\{in\}\}\}contains initial\-condition points\. The energy inner product is defined on

a:\(H1​\(I,L2​\(Ω\)\)∩L2​\(I,H2​\(Ω\)\)\)2→ℝ,a\\colon\\left\(H^\{1\}\(I,L^\{2\}\(\\Omega\)\)\\cap L^\{2\}\(I,H^\{2\}\(\\Omega\)\)\\right\)^\{2\}\\to\\mathbb\{R\},and given by

a​\(u,v\)\\displaystyle a\(u,v\)=∫01∫Ω\(∂tu−14​∂x2u\)​\(∂tv−14​∂x2v\)​dx​dt\\displaystyle=\\int\_\{0\}^\{1\}\\int\_\{\\Omega\}\\left\(\\partial\_\{t\}u\-\\frac\{1\}\{4\}\\partial\_\{x\}^\{2\}u\\right\)\\left\(\\partial\_\{t\}v\-\\frac\{1\}\{4\}\\partial\_\{x\}^\{2\}v\\right\)\\,\\mathrm\{d\}x\\mathrm\{d\}t\+∫Ωu​\(0,x\)​v​\(0,x\)​dx\+∫I×∂Ωu​v​ds​dt\.\\displaystyle\\quad\+\\int\_\{\\Omega\}u\(0,x\)v\(0,x\)\\,\\mathrm\{d\}x\+\\int\_\{I\\times\\partial\\Omega\}uv\\,\\mathrm\{d\}s\\mathrm\{d\}t\.
The heat problem tests convergence under a fixed recorded\-time budget\. We discretize the energy inner product with the loss quadrature points\. The Euclidean\-control experiment uses a width\-64 hyperbolic\-tangent network\. Table[7](https://arxiv.org/html/2607.22004#S5.T7)reports updates, runtime, and final relativeL2L^\{2\}error\. EMNGD reaches the lowest recorded error in 199\.03 seconds\.

Table 7:Audited Heat\-1D runtime and relative\-error records across methods\.Figure[18](https://arxiv.org/html/2607.22004#S5.F18)tests the same behavior atD=4,417D=4\{,\}417,5,4415\{,\}441, and99,58599\{,\}585\. The labelled EMNGD curve reaches relativeL2L^\{2\}errors near10−810^\{\-8\}within a few thousand iterations in all three settings\. The other displayed optimizers remain at higher errors over the plotted trajectories\. The runtimes are configuration\-specific and do not represent hardware\-independent complexity estimates\.

![Refer to caption](https://arxiv.org/html/2607.22004v1/x17.png)Figure 18:One\-dimensional heat relativeL2L^\{2\}error across three parameter dimensions\.

## 6Discussion

##### Feasible and energy geometry\.

EMNGD combines two distinct structures\. The parameter manifold defines admissible local variations, and retractions preserve feasibility of finite updates\. The pullback energy metric ranks the admissible variations by their function\-space effects\. EMNGD therefore minimizes the energy\-induced quadratic model directly on the tangent space\. The construction keeps the residual energy unchanged while incorporating parameter constraints throughTx​ℳT\_\{x\}\\mathcal\{M\}andRxR\_\{x\}\. Post\-hoc projection of an ambient ENGD direction generally solves a different local problem when projection and inversion do not commute\. The distinction matters whenever the manifold geometry restricts directions that have strong energy curvature in the ambient space\.

##### Geometry and linear algebra\.

The geometric definition does not depend on a particular tangent solver\. For quadratic residual energies and generalized Gauss–Newton pullbacks, Woodbury gives an exact sample\-space form of the damped tangent direction\. The identity replaces the primal system with the kernelKx=𝐉x​Πx​𝐉x⊤K\_\{x\}=\\mathbf\{J\}\_\{x\}\\Pi\_\{x\}\\mathbf\{J\}\_\{x\}^\{\\top\}and reconstructs the tangent vector throughΠx​𝐉x⊤\\Pi\_\{x\}\\mathbf\{J\}\_\{x\}^\{\\top\}\. Sample\-space natural\-gradient solves also arise in MinSR methods for variational Monte Carlo\(Chen and Heyl,[2023](https://arxiv.org/html/2607.22004#bib.bib7); Rendeet al\.,[2024](https://arxiv.org/html/2607.22004#bib.bib12)\)\. For EMNGD, Woodbury changes the linear algebra but not the tangent metric, feasible space, or retraction\.

Nyström methods have two different roles\. Sketch\-and\-solve uses a low\-rank kernel approximation and produces an approximate tangent direction\. Nyström\-preconditioned Krylov iteration uses the approximation only to accelerate the exact Woodbury system\. Iterative convergence then recovers the same damped direction as the direct sample\-space solve\. The distinction is important when interpreting accuracy and computational cost\.

Figure[19](https://arxiv.org/html/2607.22004#S6.F19)isolates the solver behavior in the overparameterized regimeN≪pN\\ll p\. Panel \(a\) compares the scaling of the primal and sample\-space formulations as the ambient dimension grows\. The sample\-space cost remains controlled by the residual count\. Panel \(b\) records primal–dual direction agreement at numerical precision and shows the dependence of a Nyström approximation on sketch rank\. Panel \(c\) compares the corresponding linearized residual trajectories\. The diagnostic separates exact Woodbury duality from acceleration strategies that approximate the kernel\.

![Refer to caption](https://arxiv.org/html/2607.22004v1/x18.png)Figure 19:Woodbury and Nyström diagnostics for tangent\-solver scaling and agreement\.
##### Operating regimes and limitations\.

Woodbury is most useful when the residual count is below the intrinsic parameter dimension\. The computational bottleneck then moves from a parameter\-space system to anN×NN\\times Nsample kernel\. Large residual sets still require quadratic kernel storage and can require expensive dense solves\. The formulation therefore does not remove the cost of residual evaluation, Jacobian products, or sample\-space conditioning\. Matrix\-free products and iterative solves become necessary when dense kernels no longer fit in memory\.

Figure[20](https://arxiv.org/html/2607.22004#S6.F20)summarizes three practical limits\. Sample\-kernel memory grows quadratically withNN\. Reducing the damping parameter approaches the undamped natural\-gradient direction but worsens conditioning and amplifies residual perturbations\. Residual subsampling replaces the full kernel by a sampled system, so the computed direction can differ from the full\-residual direction\. Nyström preconditioning is effective when the regularized kernel has a low effective dimension\. Excessive rank reduction can instead degrade the direction\. Large damping improves conditioning but moves the update toward the baseline Riemannian gradient\.

![Refer to caption](https://arxiv.org/html/2607.22004v1/x19.png)Figure 20:Sample\-kernel growth, damping sensitivity, and residual\-subsampling effects\.The residual\-Jacobian diagnostic in Figure[21](https://arxiv.org/html/2607.22004#S6.F21)provides the same perspective for the ENGD–Woodbury implementation\. The kernel spectrum spans several orders of magnitude, which explains the sensitivity of the undamped solve to small perturbations\. The result also quantifies the change in direction caused by residual batches and insufficient Nyström rank\. Such effects are numerical properties of the sampled tangent system rather than changes in the EMNGD geometry\. Damping, validation batches, and controlled Krylov tolerances provide practical safeguards, but no setting removes the underlying sample\-space trade\-off\.

![Refer to caption](https://arxiv.org/html/2607.22004v1/x20.png)Figure 21:Residual\-Jacobian sensitivity to damping, sampling, and Nyström rank\.
##### Evidence scope\.

The reported benchmarks verify Euclidean reduction, primal–dual Woodbury agreement, and high\-accuracy optimization on the stated neural PDE problems\. The residual\-formulation tests also confirm that hard boundary embeddings can be handled without using the exact PDE solution during training\. Several tables use different architectures, stopping rules, or hardware settings\. Such records provide solver coverage and endpoint evidence rather than a uniform ranking\. Most accuracy studies evaluate the Euclidean specialization\. Architecture\-matched tests on genuinely constrained neural PDE models will clarify the empirical value of the manifold component beyond the geometric guarantees\.

## 7Conclusion

We introducedEnergyManifoldNaturalGradientDescent \(EMNGD\), an intrinsic extension of ENGD from an unconstrained Euclidean parameter domain to a Riemannian parameter manifold\. EMNGD combines the feasible geometry of the parameter manifold with the energy geometry induced by the neural PDE objective\. The method solves the energy\-metric quadratic model over feasible tangent directions and uses a retraction to preserve the parameter constraints\.

Under coercivity, the push\-forward of the undamped EMNGD direction is the best feasible approximation to the function\-space Newton vector in the energy metric\. The analysis also establishes well\-posedness under damping, coordinate invariance, exact reduction to ENGD in Euclidean space, and global first\-order convergence with Armijo backtracking\. Controlled inexact tangent solves retain the descent property\. For quadratic residual energies and generalized Gauss–Newton pullbacks, the Woodbury identity gives an exact sample\-space representation of the damped tangent direction\. Nyström sketching produces an approximate direction, whereas Nyström\-preconditioned Krylov solves recover the exact direction after convergence\. The resulting solvers are most useful when the residual count is smaller than the intrinsic parameter dimension\.

Experiments verify Euclidean reduction, primal–dual agreement, and accurate neural PDE optimization under the reported settings\. The framework provides a geometric basis for constrained neural PDE solvers without altering the residual energy\. Future work can evaluate architecture\-matched constrained manifolds and develop matrix\-free solvers for larger residual systems\.

Acknowledgments and Disclosure of Funding

The study was supported by the National Natural Science Foundation of China \(12572138\)\. The authors declare no competing interests\.

## References

- Information geometry of divergence functions\.Bulletin of the Polish academy of sciences\. Technical sciences58\(1\),pp\. 183–195\.Cited by:[§3](https://arxiv.org/html/2607.22004#S3.p1.1)\.
- S\. Amari \(1996\)Neural learning in structured parameter spaces\-natural riemannian gradient\.Advances in neural information processing systems9\.Cited by:[§2](https://arxiv.org/html/2607.22004#S2.SS0.SSS0.Px5.p1.8),[§3](https://arxiv.org/html/2607.22004#S3.p1.1)\.
- S\. Amari \(2016\)Information geometry and its applications\.Vol\.194,Springer,Japan\.Cited by:[§3](https://arxiv.org/html/2607.22004#S3.p1.1),[footnote 1](https://arxiv.org/html/2607.22004#footnote1)\.
- J\. A\. Bagnell and J\. G\. Schneider \(2003\)Covariant policy search\.InIJCAI,pp\. 1019–1024\.Cited by:[§3](https://arxiv.org/html/2607.22004#S3.p1.1)\.
- C\. Beck, M\. Hutzenthaler, A\. Jentzen, and B\. Kuckuck \(2020\)An overview on deep learning\-based approximation methods for partial differential equations\.arXiv preprint arXiv:2012\.12348\.Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p1.1),[§2](https://arxiv.org/html/2607.22004#S2.p1.1)\.
- J\. Bradbury, R\. Frostig, P\. Hawkins, M\. J\. Johnson, C\. Leary, D\. Maclaurin, G\. Necula, A\. Paszke, J\. VanderPlas, S\. Wanderman\-Milne, and Q\. Zhang \(2018\)JAX: composable transformations of Python\+NumPy programs\.External Links:[Link](http://github.com/google/jax)Cited by:[§5\.1\.2](https://arxiv.org/html/2607.22004#S5.SS1.SSS2.p1.1)\.
- T\. Cai, R\. Gao, J\. Hou, S\. Chen, D\. Wang, D\. He, Z\. Zhang, and L\. Wang \(2019\)Gram\-gauss\-newton method: learning overparameterized neural networks for regression problems\.arXiv preprint arXiv:1905\.11675\.Cited by:[§3](https://arxiv.org/html/2607.22004#S3.p3.2)\.
- A\. Chen and M\. Heyl \(2023\)Efficient optimization of deep neural quantum states toward machine precision\.arXiv preprint arXiv:2302\.01941\.Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p4.1),[§6](https://arxiv.org/html/2607.22004#S6.SS0.SSS0.Px2.p1.2)\.
- L\. Courte and M\. Zeinhofer \(2023\)Robin Pre\-Training for the Deep Ritz Method\.Northern Lights Deep Learning Conference\.Cited by:[§2](https://arxiv.org/html/2607.22004#S2.SS0.SSS0.Px2.p2.1)\.
- F\. Dangel, J\. Müller, and M\. Zeinhofer \(2024\)Kronecker\-factored approximate curvature for physics\-informed neural networks\.arXiv preprint arXiv:2405\.15603\.Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p3.1)\.
- C\. Davi and U\. Braga\-Neto \(2022\)PSO\-pinn: physics\-informed neural networks trained with particle swarm optimization\.arXiv preprint arXiv:2202\.01943\.Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p2.1)\.
- A\. Daw, J\. Bu, S\. Wang, P\. Perdikaris, and A\. Karpatne \(2022\)Rethinking the importance of sampling in physics\-informed neural networks\.arXiv preprint arXiv:2207\.02338\.Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p2.1)\.
- M\. Dissanayake and N\. Phan\-Thien \(1994\)Neural\-network\-based approximations for solving partial differential equations\.communications in Numerical Methods in Engineering10\(3\),pp\. 195–201\.Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p1.1),[§2](https://arxiv.org/html/2607.22004#S2.SS0.SSS0.Px1.p2.1)\.
- W\. E, J\. Han, and A\. Jentzen \(2017\)Deep learning\-based numerical methods for high\-dimensional parabolic partial differential equations and backward stochastic differential equations\.Communications in Mathematics and Statistics5\(4\),pp\. 349–380\.Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p1.1)\.
- W\. E and B\. Yu \(2018\)The Deep Ritz Method: A Deep Learning\-Based Numerical Algorithm for Solving Variational Problems\.Communications in Mathematics and Statistics6\(1\),pp\. 1–12\.Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p1.1),[§2](https://arxiv.org/html/2607.22004#S2.SS0.SSS0.Px2.p1.7)\.
- Z\. Frangella, J\. A\. Tropp, and M\. Udell \(2023\)Randomized Nyström preconditioning\.SIAM Journal on Matrix Analysis and Applications44\(2\),pp\. 718–752\.Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p4.1)\.
- M\. Gargiani, A\. Zanelli, M\. Diehl, and F\. Hutter \(2020\)On the promise of the stochastic generalized gauss\-newton method for training dnns\.arXiv preprint arXiv:2006\.02409\.Cited by:[§3](https://arxiv.org/html/2607.22004#S3.p3.2)\.
- A\. Gittens and M\. W\. Mahoney \(2016\)Revisiting the Nyström method for improved large\-scale machine learning\.Journal of Machine Learning Research17\(117\),pp\. 1–65\.Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p4.1)\.
- G\. Goldshlager, N\. Abrahamsen, and L\. Lin \(2024\)A Kaczmarz\-inspired approach to accelerate the optimization of neural network wavefunctions\.Journal of Computational Physics516,pp\. 113351\.Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p4.1),[§3\.1](https://arxiv.org/html/2607.22004#S3.SS1.SSS0.Px1.p8.3)\.
- J\. Han, A\. Jentzen, and E\. Weinan \(2018\)Solving high\-dimensional partial differential equations using deep learning\.Proceedings of the National Academy of Sciences115\(34\),pp\. 8505–8510\.Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p1.1)\.
- W\. Hao, X\. Jin, J\. W\. Siegel, and J\. Xu \(2021\)An efficient greedy training algorithm for neural networks and applications in PDEs\.arXiv preprint arXiv:2107\.04466\.Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p2.1)\.
- A\. Jacot, F\. Gabriel, and C\. Hongler \(2018\)Neural tangent kernel: Convergence and generalization in neural networks\.InAdvances in neural information processing systems,pp\. 8571–8580\.Cited by:[§3\.1](https://arxiv.org/html/2607.22004#S3.SS1.p6.3)\.
- A\. Jnini, F\. Vella, and M\. Zeinhofer \(2024\)Gauss\-Newton natural gradient descent for physics\-informed computational fluid dynamics\.arXiv preprint arXiv:2402\.10680\.Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p3.1)\.
- S\. M\. Kakade \(2001\)A natural policy gradient\.Advances in Neural Information Processing Systems14\.Cited by:[§3](https://arxiv.org/html/2607.22004#S3.p1.1)\.
- N\. Kovachki, Z\. Li, B\. Liu, K\. Azizzadenesheli, K\. Bhattacharya, A\. Stuart, and A\. Anandkumar \(2021\)Neural operator: learning maps between function spaces\.arXiv preprint arXiv:2108\.08481\.Cited by:[§2](https://arxiv.org/html/2607.22004#S2.p1.1)\.
- A\. Krishnapriyan, A\. Gholami, S\. Zhe, R\. Kirby, and M\. W\. Mahoney \(2021\)Characterizing possible failure modes in physics\-informed neural networks\.Advances in Neural Information Processing Systems34,pp\. 26548–26560\.Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p2.1)\.
- I\. E\. Lagaris, A\. Likas, and D\. I\. Fotiadis \(1998\)Artificial neural networks for solving ordinary and partial differential equations\.IEEE transactions on neural networks9\(5\),pp\. 987–1000\.Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p1.1),[§2](https://arxiv.org/html/2607.22004#S2.SS0.SSS0.Px1.p2.1)\.
- W\. Li and G\. Montúfar \(2018\)Natural gradient via optimal transport\.Information Geometry1\(2\),pp\. 181–214\.External Links:ISBN 2511\-249X,[Link](https://doi.org/10.1007/s41884-018-0015-3)Cited by:[§3](https://arxiv.org/html/2607.22004#S3.p1.1)\.
- Z\. Li, N\. B\. Kovachki, K\. Azizzadenesheli, B\. liu, K\. Bhattacharya, A\. Stuart, and A\. Anandkumar \(2021\)Fourier neural operator for parametric partial differential equations\.InInternational Conference on Learning Representations,External Links:[Link](https://openreview.net/forum?id=c8P9NQVtmnO)Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p1.1)\.
- A\. T\. Lin, W\. Li, S\. Osher, and G\. Montúfar \(2021\)Wasserstein proximal of gans\.InInternational Conference on Geometric Science of Information,pp\. 524–533\.Cited by:[§3](https://arxiv.org/html/2607.22004#S3.p1.1)\.
- L\. Lu, X\. Meng, Z\. Mao, and G\. E\. Karniadakis \(2021\)DeepXDE: a deep learning library for solving differential equations\.SIAM Review63\(1\),pp\. 208–228\.Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p2.1)\.
- J\. Martens \(2020\)New insights and perspectives on the natural gradient method\.The Journal of Machine Learning Research21\(1\),pp\. 5776–5851\.Cited by:[§3](https://arxiv.org/html/2607.22004#S3.p1.1),[§3](https://arxiv.org/html/2607.22004#S3.p3.2)\.
- T\. Morimura, E\. Uchibe, J\. Yoshimoto, and K\. Doya \(2008\)A new natural policy gradient by stationary distribution metric\.InJoint European Conference on Machine Learning and Knowledge Discovery in Databases,pp\. 82–97\.Cited by:[§3](https://arxiv.org/html/2607.22004#S3.p1.1)\.
- J\. Müller and G\. Montúfar \(2022\)Geometry and convergence of natural policy gradients\.MPI MiS Preprint 31/2022\.External Links:[Link](https://www.mis.mpg.de/publications/preprints/2022/prepr2022-31.html)Cited by:[§3](https://arxiv.org/html/2607.22004#S3.p1.1)\.
- J\. Müller and M\. Zeinhofer \(2022a\)Error estimates for the deep ritz method with boundary penalty\.InMathematical and Scientific Machine Learning,pp\. 215–230\.Cited by:[§2](https://arxiv.org/html/2607.22004#S2.SS0.SSS0.Px2.p2.1)\.
- J\. Müller and M\. Zeinhofer \(2022b\)Notes on exact boundary values in residual minimisation\.InMathematical and Scientific Machine Learning,pp\. 231–240\.Cited by:[footnote 3](https://arxiv.org/html/2607.22004#footnote3)\.
- J\. Müller and M\. Zeinhofer \(2023\)Achieving high accuracy with PINNs via energy natural gradient descent\.InInternational Conference on Machine Learning,pp\. 25471–25485\.Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p3.1),[§2](https://arxiv.org/html/2607.22004#S2.SS0.SSS0.Px5.p1.8)\.
- J\. Müller and M\. Zeinhofer \(2024\)Position: optimization in SciML should employ the function space geometry\.InForty\-first International Conference on Machine Learning,Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p3.1)\.
- M\. A\. Nabian, R\. J\. Gladstone, and H\. Meidani \(2021\)Efficient training of physics\-informed neural networks via importance sampling\.Computer\-Aided Civil and Infrastructure Engineering36\(8\),pp\. 962–977\.Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p2.1)\.
- H\. Nie, B\. Gao, A\. Han, P\. Jawanpuria, B\. Mishra, and Z\. Wen \(2026\)Nyström approximation on manifolds\.arXiv preprint arXiv:2605\.14933\.Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p4.1),[§3\.1](https://arxiv.org/html/2607.22004#S3.SS1.SSS0.Px1.p6.6)\.
- J\. Nocedal and S\. J\. Wright \(1999\)Numerical optimization\.Springer\.Cited by:[§5\.1\.1](https://arxiv.org/html/2607.22004#S5.SS1.SSS1.p1.7)\.
- R\. Novak, J\. Sohl\-Dickstein, and S\. S\. Schoenholz \(2022\)Fast finite width neural tangent kernel\.InInternational Conference on Machine Learning,pp\. 17018–17044\.Cited by:[§3\.1](https://arxiv.org/html/2607.22004#S3.SS1.p6.3)\.
- L\. Nurbekyan, W\. Lei, and Y\. Yang \(2022\)Efficient natural gradient descent methods for large\-scale optimization problems\.arXiv:2202\.06236\.Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p3.1),[§3](https://arxiv.org/html/2607.22004#S3.p1.1),[§3](https://arxiv.org/html/2607.22004#S3.p11.3.2)\.
- R\. Pascanu and Y\. Bengio \(2014\)Revisiting natural gradient for deep networks\.InInternational Conference on Learning Representations,External Links:[Link](https://openreview.net/forum?id=vz8AumxkAfz5U)Cited by:[§3](https://arxiv.org/html/2607.22004#S3.p1.1)\.
- J\. Peters, S\. Vijayakumar, and S\. Schaal \(2003\)Reinforcement learning for humanoid robotics\.InProceedings of the third IEEE\-RAS international conference on humanoid robots,pp\. 1–20\.Cited by:[§3](https://arxiv.org/html/2607.22004#S3.p1.1)\.
- M\. Raissi, P\. Perdikaris, and G\. E\. Karniadakis \(2019\)Physics\-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations\.Journal of Computational physics378,pp\. 686–707\.Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p1.1),[§2](https://arxiv.org/html/2607.22004#S2.SS0.SSS0.Px1.p2.1)\.
- Y\. Ren and D\. Goldfarb \(2019\)Efficient subsampled gauss\-newton and natural gradient methods for training neural networks\.arXiv preprint arXiv:1906\.02353\.Cited by:[§3](https://arxiv.org/html/2607.22004#S3.p3.2)\.
- R\. Rende, L\. L\. Viteritti, L\. Bardone, F\. Becca, and S\. Goldt \(2024\)A simple linear algebra identity to optimize large\-scale neural network quantum states\.Communications Physics7\(1\),pp\. 260\.Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p4.1),[§6](https://arxiv.org/html/2607.22004#S6.SS0.SSS0.Px2.p1.2)\.
- W\. Ritz \(1909\)Über eine neue Methode zur Lösung gewisser Variationsprobleme der mathematischen Physik\.\.Journal für die reine und angewandte Mathematik \(Crelles Journal\)1909\(135\),pp\. 1–61\.Cited by:[§2](https://arxiv.org/html/2607.22004#S2.SS0.SSS0.Px2.p1.7)\.
- N\. N\. Schraudolph \(2002\)Fast curvature matrix\-vector products for second\-order gradient descent\.Neural computation14\(7\),pp\. 1723–1738\.Cited by:[§3](https://arxiv.org/html/2607.22004#S3.p1.1)\.
- T\. Schwedes, S\. W\. Funke, and D\. A\. Ham \(2016\)An iteration count estimate for a mesh\-dependent steepest descent method based on finite elements and Riesz inner product representation\.arXiv preprint arXiv:1606\.08069\.Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p3.1)\.
- T\. Schwedes, D\. A\. Ham, S\. W\. Funke, and M\. D\. Piggott \(2017\)Mesh dependence in PDE\-constrained optimisation\.InMesh Dependence in PDE\-Constrained Optimisation,pp\. 53–78\.Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p3.1)\.
- Z\. Shen, Z\. Wang, A\. Ribeiro, and H\. Hassani \(2020\)Sinkhorn natural gradient for generative models\.Advances in Neural Information Processing Systems33,pp\. 1646–1656\.Cited by:[§3](https://arxiv.org/html/2607.22004#S3.p1.1)\.
- J\. Sirignano and K\. Spiliopoulos \(2018\)DGM: A deep learning algorithm for solving partial differential equations\.Journal of computational physics375,pp\. 1339–1364\.Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p1.1),[§2](https://arxiv.org/html/2607.22004#S2.SS0.SSS0.Px1.p2.1)\.
- R\. van der Meer, C\. W\. Oosterlee, and A\. Borovykh \(2022\)Optimally weighted loss functions for solving pdes with neural networks\.Journal of Computational and Applied Mathematics405,pp\. 113887\.Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p2.1)\.
- J\. van Oostrum, J\. Müller, and N\. Ay \(2022\)Invariance properties of the natural gradient in overparametrised systems\.Information Geometry,pp\. 1–17\.Cited by:[footnote 1](https://arxiv.org/html/2607.22004#footnote1)\.
- L\. Wang and M\. Yan \(2022\)Hessian informed mirror descent\.Journal of Scientific Computing92\(3\),pp\. 1–22\.Cited by:[§3](https://arxiv.org/html/2607.22004#S3.p1.1)\.
- S\. Wang, S\. Sankaran, and P\. Perdikaris \(2022a\)Respecting causality is all you need for training physics\-informed neural networks\.arXiv preprint arXiv:2203\.07404\.Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p2.1)\.
- S\. Wang, Y\. Teng, and P\. Perdikaris \(2021\)Understanding and mitigating gradient flow pathologies in physics\-informed neural networks\.SIAM Journal on Scientific Computing43\(5\),pp\. A3055–A3081\.Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p2.1),[§2](https://arxiv.org/html/2607.22004#S2.SS0.SSS0.Px4.p1.1)\.
- S\. Wang, X\. Yu, and P\. Perdikaris \(2022b\)When and why PINNs fail to train: a neural tangent kernel perspective\.Journal of Computational Physics449,pp\. 110768\.Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p2.1)\.
- E\. Weinan, J\. Han, and A\. Jentzen \(2021\)Algorithms for solving high dimensional PDEs: from nonlinear monte carlo to machine learning\.Nonlinearity35\(1\),pp\. 278\.Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p1.1),[§2](https://arxiv.org/html/2607.22004#S2.p1.1)\.
- C\. Wu, M\. Zhu, Q\. Tan, Y\. Kartha, and L\. Lu \(2023\)A comprehensive study of non\-adaptive and residual\-based adaptive sampling for physics\-informed neural networks\.Computer Methods in Applied Mechanics and Engineering403,pp\. 115671\.Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p2.1)\.
- B\. Zapf, J\. Haubner, M\. Kuchta, G\. Ringstad, P\. K\. Eide, and K\. Mardal \(2022\)Investigating molecular transport in the human brain from mri with physics\-informed neural networks\.Scientific Reports12\(1\),pp\. 1–12\.Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p2.1)\.
- Q\. Zeng, S\. H\. Bryngelson, and F\. T\. Schaefer \(2022\)Competitive physics informed networks\.InICLR 2022 Workshop on Gamification and Multiagent Solutions,External Links:[Link](https://openreview.net/forum?id=rMz_scJ6lc)Cited by:[§1](https://arxiv.org/html/2607.22004#S1.p2.1),[§2](https://arxiv.org/html/2607.22004#S2.SS0.SSS0.Px4.p1.1)\.