math.NAMay 16, 2025

Accelerating Natural Gradient Descent for PINNs with Randomized Numerical Linear Algebra

Authors: Ivan Bioli, Carlo Marcati, Giancarlo Sangalli

Organizations: Department of Mathematics, University of Pavia, Via A. Ferrata 5, 27100 Pavia, Italy · Department of Civil Engineering and Architecture, University of Pavia, Via A. Ferrata 3, 27100 Pavia, Italy · Istituto di Matematica Applicata e Tecnologie Informatiche “E. Magenes”, CNR, Via A. Ferrata 1, 27100 Pavia, Italy

Read PDF
Ask AI
Cite

Copy citation

DOINot available for this paper

Abstract

Natural Gradient Descent (NGD) has emerged as a promising optimization algorithm for training neural network-based solvers for partial differential equations (PDEs), such as Physics-Informed Neural Networks (PINNs). However, its practical use is often limited by the high computational cost of solving linear systems involving the Gramian matrix. While matrix-free NGD methods based on the conjugate gradient (CG) method avoid explicit matrix inversion, the ill-conditioning of the Gramian significantly slows the convergence of CG. In this work, we extend matrix-free NGD to broader classes of problems than previously considered and propose the use of Randomized Numerical Linear Algebra (RandNLA) techniques for efficient preconditioning of the inner CG solver. The resulting algorithms demonstrate substantial performance improvements over existing NGD-based methods on a range of PDE problems discretized using neural networks, and offer competitive results compared to other state-of-the-art optimizers.

Figures & tables

StrategySetup CostMMP CostMemory Footprint
Precompute GramianO(q[FP]+p2q)\mathcal{O}(q\texttt{[FP]}+p^{2}q)O(p2ℓ)\mathcal{O}(p^{2}\ell)O(p2+pB)\mathcal{O}(p^{2}+pB)
Precompute Jacobian onlyO(q[FP])\mathcal{O}(q\texttt{[FP]})O(pqℓ)\mathcal{O}(pq\ell)O(pq)\mathcal{O}(pq)
Recompute batch Jacobian   
on the fly–O(q[FP]+pqℓ)\mathcal{O}(q\texttt{[FP]}+pq\ell)O(pB)\mathcal{O}(pB)
Table 1 . Summary of computational and memory costs for the three matrix-vector product strategies. Here pp is the number of parameters, [FP] is the cost of a forward pass, ℓ\ell is the number of simultaneous matvecs, qq is the number of quadrature points, and BB is the batch size for processing quadrature points. We assume that p≤qp\leq q , i.e., we focus on the regime where a highly accurate quadrature rule is employed and the number of quadrature points qq exceeds the number of neural network parameters pp .
OptimizerPoisson 3DHeat 3+1DKovasznayDeep-Ritz PoissonFEINNs 2D
NGD full1.56×10−61.56\times 10^{-6}    
(1.37−1.77 ×10−6)(1.37-1.77\,\times 10^{-6})     
317.7 s1.75×10−61.75\times 10^{-6}    
(1.49−1.86 ×10−6)(1.49-1.86\,\times 10^{-6})     
279.2 s1.08×10−61.08\times 10^{-6}    
(0.965−1.25 ×10−6)(0.965-1.25\,\times 10^{-6})     
Table 2 . Relative H1(Ω)\mathrm{H}^{1}(\Omega) error at the plateau iterate for each optimizer and test problem. The first line reports the median error, the second line reports in smaller font and in parentheses the 2525 th and 7575 th percentiles in compact form as (a−b ×10k)(a-b\,\times 10^{k}) , and the third line reports the corresponding cumulative wall-clock time in seconds. The plateau criterion and the dot ( ∙\bullet ) marking the corresponding iterate in the convergence curves are described in Section 5.2 ; see Figure 10 .
Optimizer3×64×64×64×1\mathbf{3\times 64\times 64\times 64\times 1}
q=104\mathbf{q=10^{4}}3×64×64×64×1\mathbf{3\times 64\times 64\times 64\times 1}
q=105\mathbf{q=10^{5}}3×128×128×128×1\mathbf{3\times 128\times 128\times 128\times 1}
q=105\mathbf{q=10^{5}} 
NGD full1.56×10−61.56\times 10^{-6}
(1.37−1.77 ×10−6)(1.37-1.77\,\times 10^{-6}) 
317.7 s 
Table 3 . Scalability experiments on the 3D Poisson problem discretized via PINNs; see Section 5.3 . Relative H1(Ω)\mathrm{H}^{1}(\Omega) error at the plateau iterate for each optimizer and architecture. The first line reports the median error, the second line reports in smaller font and in parentheses the 2525 th and 7575 th percentiles in compact form as (a−b ×10k)(a-b\,\times 10^{k}) , and the third line reports the corresponding cumulative wall-clock time in seconds, while the fourth line reports the per-iteration cost in seconds. Each column corresponds to a different network architecture and number of quadrature points qq , as indicated in the header.
Appendix figures & tables6 assets

Supplementary material from the paper’s appendix.

Appendix

ExperimentNGD (all)(L-)BFGSSSBroydenAdamNetwork Architecture
Poisson 3D1 00050 00050 000150 000[3, 64, 64, 64, 1]
Heat 3+1D1 00050 00050 000150 000[4, 64, 64, 64, 1]
Kovasznay2 00025 00025 00075 000uu : [2, 45, 45, 45, 2]
pp : [2, 45, 45, 45, 1]     
Deep-Ritz Poisson1 00025 00025 00075 000[2, 64, 64, 64, 1]
FEINNs 2D1 00025 00025 00075 000[2, 64, 64, 64, 1]
Appendix
Table 4 . Iteration budgets for the evaluated optimizers and neural network architectures for each PDE experiment. Network shapes are denoted as [input_dim, hidden_1, ..., hidden_N, output_dim] .
ExperimentInterior / Mesh configurationBoundary PointsEvaluation set
Poisson 3D10 000 random1 000 random100 000 random
Heat 3+1D10 000 random1 000 random100 000 random
Kovasznay10 000 random1 000 random100 000 random
Deep-Ritz PoissonQuad mesh, h=1/30h=1/30 , d=4d=4 (14.5k DoF)-Gauss–Jacobi, order 1616
FEINNs 2DQuad mesh, h=1/30h=1/30 , d=4d=4 (14.5k DoF)-Gauss–Jacobi, order 1616
Appendix
Table 5 . Domain discretization, boundary conditions, and evaluation setups for each numerical experiment.
γ\boldsymbol{\gamma}Poisson 3DHeat 3+1DDeep-Ritz PoissonFEINNs 2D
0.1\mathbf{0.1}1.87×10−61.87\times 10^{-6}   
(0.886−3.27 ×10−6)(0.886-3.27\,\times 10^{-6})    
1.004×1.004\times7.42×10−67.42\times 10^{-6}   
(6.05−8.60 ×10−6)(6.05-8.60\,\times 10^{-6})    
1.023×1.023\times7.49×10−17.49\times 10^{-1}   
(3.99−10.4 ×10−1)(3.99-10.4\,\times 10^{-1})    
Appendix
Table 6 . Sensitivity analysis of γ\gamma for NyströmNGD with Gaussian test matrices. The first line reports the median relative H1(Ω)\mathrm{H}^{1}(\Omega) error, the second line reports in smaller font and in parentheses the 2525 th and 7575 th percentiles in compact form as (a−b ×10k)(a-b\,\times 10^{k}) , and the third line reports the wall-clock time relative to the default γ=10\gamma=10 .
γ\boldsymbol{\gamma}Poisson 3DHeat 3+1DDeep-Ritz PoissonFEINNs 2D
0.1\mathbf{0.1}8.42×10−38.42\times 10^{-3}   
(7.79−8.49 ×10−3)(7.79-8.49\,\times 10^{-3})    
1.062×1.062\times7.48×10−37.48\times 10^{-3}   
(5.24−9.63 ×10−3)(5.24-9.63\,\times 10^{-3})    
1.055×1.055\times2.122.12   
(1.08−2.53)(1.08-2.53)    
Appendix
Table 7 . Sensitivity analysis of γ\gamma for RPCholNGD. The first line reports the median relative H1(Ω)\mathrm{H}^{1}(\Omega) error, the second line reports in smaller font and in parentheses the 2525 th and 7575 th percentiles in compact form as (a−b ×10k)(a-b\,\times 10^{k}) , and the third line reports the wall-clock time relative to the default γ=10\gamma=10 .
r\boldsymbol{r}Poisson 3DHeat 3+1DKovasznayDeep-Ritz PoissonFEINNs 2D
0.1\mathbf{0.1}1.68×10−61.68\times 10^{-6}    
(1.55−1.75 ×10−6)(1.55-1.75\,\times 10^{-6})     
1.00×1.00\times2.03×10−62.03\times 10^{-6}    
(2.03−2.12 ×10−6)(2.03-2.12\,\times 10^{-6})     
1.00×1.00\times6.95×10−76.95\times 10^{-7}    
(6.80−7.48 ×10−7)(6.80-7.48\,\times 10^{-7})     
Appendix
Table 8 . Sensitivity analysis of the parameter rr for NyströmNGD with Gaussian test matrices. The first line reports the median relative H1(Ω)\mathrm{H}^{1}(\Omega) error, the second line reports in smaller font and in parentheses the 2525 th and 7575 th percentiles in compact form as (a−b ×10k)(a-b\,\times 10^{k}) , and the third line reports the wall-clock time relative to the default r=10r=10 .
r\boldsymbol{r}Poisson 3DHeat 3+1DKovasznayDeep-Ritz PoissonFEINNs 2D
0.01\mathbf{0.01}1.83×10−61.83\times 10^{-6}    
(1.78−1.87 ×10−6)(1.78-1.87\,\times 10^{-6})     
1.00×1.00\times2.11×10−62.11\times 10^{-6}    
(1.81−2.53 ×10−6)(1.81-2.53\,\times 10^{-6})     
1.00×1.00\times9.79×10−79.79\times 10^{-7}    
(9.14−10.3 ×10−7)(9.14-10.3\,\times 10^{-7})     
Appendix
Table 9 . Sensitivity analysis of the RPCholesky stopping tolerance scaling factor rr ( ε=r⋅p⋅ϵmach\varepsilon=r\cdot p\cdot\epsilon_{\mathrm{mach}} ) for RPCholNGD. The first line reports the median relative H1(Ω)\mathrm{H}^{1}(\Omega) error, the second line reports in smaller font and in parentheses the 2525 th and 7575 th percentiles in compact form as (a−b ×10k)(a-b\,\times 10^{k}) , and the third line reports the wall-clock time relative to the default r=1r=1 .

Explore similar work

CardsList
Jul 24, 2026cs.LG

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

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 introduce \EMNGDfull{}, 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.

Zhangyong Liang, Huanhuan Gao

National Center for Applied Mathematics Tianjin University Tianjin, 300072, China · School of Mechanical and Aerospace Engineering Jilin University Changchun, 130025, China

Neural Partial Differential Equation SolversRiemannian OptimizationRiemannian ManifoldsGradient DescentManifolds
AskPDF
Jul 27, 2026cs.LG

Global Convergence of DGM and PINN Algorithms for Solving Nonlinear PDEs

The Deep Galerkin Method (DGM) and Physics Informed Neural Networks (PINNs) have become widely-used methods for solving partial differential equations (PDEs) in the rapidly growing field of scientific machine learning. In these methods, a neural network is trained to approximate the PDE solution by using (stochastic) gradient descent to minimize the PDE residual of the neural network. Due to the non-convexity of the PDE residual objective function, the trained neural network may, in principle, only converge to a local minimizer of the objective function (which would not be a solution of the PDE). Therefore, there is a longstanding question regarding the mathematical foundations of these algorithms, and it is highly valuable to establish that the trained neural network will converge to the PDE solution. In this paper, we consider a class of semilinear PDEs with nonlinearities in the solution and its first derivative. For this class of PDEs, we prove that neural networks trained with gradient descent to minimize the PDE residual objective function will converge to the PDE solution as the network width and training time →∞\rightarrow \infty.

Justin Sirignano, Konstantinos Spiliopoulos, Samuel Cohen

Mathematical Institute, University of Oxford · Department of Mathematics and Statistics, Boston University

Partial Differential EquationsVariational FormulationConvergenceGradient Descent
AskPDF
May 24, 2026cs.LG

Mitigating Gradient Pathology in PINNs through Aligned Constraint

While Physics-Informed Neural Networks (PINNs) are powerful for solving Partial Differential Equations (PDEs), their training is often paralyzed by gradient pathology. The gradients from the PDE residuals and boundary constraints oppose each other, trapping the model in local minima. Current solutions, such as adaptive weighting or hard constraints, either fail to fundamentally resolve this ill-conditioning or are limited to simple geometries. In this study, we systematically analyze the possible causes of this gradient pathology from the perspectives of loss landscapes and optimization dynamics. Based on the obtained conclusion, we propose Constraint-Aligned loss with Manifold Lifting (CAML). By reformulating all zeroth-order terms into aligned constraints, our method effectively mitigates gradient conflicts. In addition, we introduce a delay factor to help the optimizer skip the high-curvature area. Experiments demonstrate that our CAML significantly enhances numerical stability and efficiency in highly complex PINN problems. Our code is open-sourced on https://github.com/YichenLuo-0/CAML.

Yichen Luo, Peiyu Zhu, Dongxiao Hu +5

Department of Information Science and Engineering, KTH Royal Institute of Technology, Stockholm, Sweden · School of Advanced Manufacturing and Robotics, Peking University, Beijing, China · School of Advanced Technology, Xi’an Jiaotong-Liverpool University, Suzhou, China +2

Parametric Physics-Informed Neural NetworkPartial Differential EquationsStructural ConstraintsManifolds
AskPDF