Neural network-based solvers for partial differential equations (PDEs), such as Physics-Informed Neural Networks (PINNs), offer mesh-free approximations and hold great promise for high-dimensional problems. However, their practical use is often limited by optimization difficulties due to the ill-conditioning of the loss landscape. Recently, second-order optimizers derived from preconditioning the underlying infinite-dimensional PDE operator have shown promise, particularly Natural Gradient Descent (NGD) [3]. NGD uses the Gramian matrix as a preconditioner in gradient descent, but suffers from high computational cost due to the cubic complexity of solving associated linear systems with direct methods. 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 problem classes and propose using Randomized Nyström preconditioning [2] to accelerate convergence of the inner CG solver. The resulting algorithm leverages the rapid spectral decay of the Gramian to construct an efficient low-rank approximation, which then serves as a preconditioner for the CG iterations. This turns the strong spectral decay of the Gramian from a challenge to a strength. Our approach significantly reduces the cost of the inner solver, achieving faster convergence in both iterations and wall-clock time. We validate our method on different PDE problems discretized using PINNs, demonstrating that it consistently outperforms previous NGD variants by achieving higher accuracy at lower computational cost. References [1] Bioli, C. Marcati, and G. Sangalli. Accelerating Natural Gradient Descent for PINNs with Randomized Numerical Linear Algebra. [2] Z. Frangella, J. A. Tropp, and M. Udell. “Randomized Nystr¨ om Preconditioning”. In: SIAM Journal on Matrix Analysis and Applications 44.2 (). Publisher: Society for Industrial and Applied Mathematics, pp. 718–752. [3] J. Müller and M. Zeinhofer. Position: Optimization in SciML Should Employ the Function Space Geometry. arXiv: 2402.07318.