An overview of machine learning-enhanced iterative methods for systems of linear and nonlinear equations

Yuhuang Meng; Jing Zhao; Alexander Heinlein

math.NAOct 5, 2026

An overview of machine learning-enhanced iterative methods for systems of linear and nonlinear equations

Authors: Yuhuang Meng, Jing Zhao, Alexander Heinlein

Organizations: Delft Institute of Applied Mathematics, Delft University of Technology, Mekelweg 4, Delft, 2628 CD, the Netherlands · Department of Hydrodynamics and Forecasting, Deltares, Boussinesqweg 1, Delft, 2629 HV, the Netherlands

Abstract

Systems of equations arise in a wide range of scientific and engineering applications. The present work focuses on solvers for general systems of equations, including but not limited to those arising from partial differential equations. These systems can be broadly categorized into linear and nonlinear problems. For large linear systems, iterative solvers are generally preferred over direct methods due to the latter's superlinear growth of computational costs. Although convergence theory is well-developed under certain assumptions on the coefficient matrix, many classes of systems still pose open challenges. These difficulties become even more severe for systems of nonlinear equations, where nonlinear solvers typically rely on repeated linearization. For example, Newton's method may even converge quadratically near the solution; it can also converge slowly or diverge when the initial guess is not chosen appropriately. A wide range of solvers with diverse variants and hyperparameter settings exists, and the development of efficient and robust iterative methods remains an active area of research. Recently, machine learning (ML) techniques have been applied to enhance the efficiency of classical iterative methods while preserving their interpretability and reliability. We refer to these ML-enhanced iterative methods as hybrid iterative methods, in the sense that they combine classical iterative methods with ML. This paper provides a comprehensive overview of state-of-the-art approaches to constructing hybrid iterative methods for systems of both linear and nonlinear equations, while also discussing open challenges and outlining potential directions for future research.

Figures & tables

Figure 1 : Overview of ML-enhanced iterative methods, referred to in this paper as hybrid iterative methods (HIMs) . (a) General workflow of HIMs for solving linear and nonlinear systems of equations. (b) Classification of HIMs for linear systems into four categories, as discussed in Section 3 . The classification of HIMs for nonlinear systems, discussed in Section 4 , follows a similar structure.

Figure 2 : Schematic of a fully connected feedforward neural network (FNN) with two hidden layers. Nonlinearity is introduced through neuron activation functions, while the output layer is often linear.

Figure 3 : Illustration of a 2D convolution and max pooling. The input

X

has size

5\times 5

and is extended to

7\times 7

by zero padding, analogous to homogeneous Dirichlet boundary conditions, whereas padding with repeated boundary values resembles homogeneous Neumann boundary conditions. A convolution with a

2\times 2

kernel

K

produces a

6\times 6

feature map, which is then downsampled by taking the maximum over each

3\times 3

subblock.

Figure 4 : Illustration of the U-Net architecture, taken from PlotNeuralNet [ 42 ] . The architecture consists of an encoder path and a decoder path, where corresponding refinement levels are connected by skip connections. This structure resembles the process in a multigrid V-cycle.

Figure 5 : Applications of CNNs and GNNs to the linear system

A\mathbf{x}=\mathbf{b}

. (a) The coefficient matrix

A

is represented as an image and processed by a CNN, where blank cells denote zero entries of

A

. (b) The matrix

A

is represented as a graph, in which matrix entries define edge connections and attributes, and is processed by a GNN. (c) For a PDE discretized on a uniform Cartesian grid, the solution

\mathbf{x}

is represented as a grid-based field and can be naturally processed by CNNs. (d) For a PDE discretized on an irregular (also called unstructured) mesh, the solution is represented on a graph and can be handled more flexibly by GNNs. The red star denotes a representative degree of freedom, and the red nodes and connecting edges indicate its local neighborhood. Thus, CNNs exploit the fixed spatial arrangement of regular grids, whereas GNNs aggregate information over general connectivity patterns.

Figure 6 : Illustration of a recurrent neural network (RNN). At each step

k

, the shared network

\mathcal{N}_{\theta}

updates the hidden state

\mathbf{z}_{k}

from the previous hidden state

\mathbf{z}_{k-1}

and the current input

\mathbf{x}_{k}

, and produces the output

\mathbf{y}_{k}

. Here,

\mathcal{N}_{\theta}

represents the mappings

g_{1}

and

g_{2}

in Equation 4 .

method class	initial guess	update function	parameters
stationary methods	determines the initial error	$\mathbf{x}_{k+1}=\mathbf{x}_{k}+M^{-1}(\mathbf{b}-A\mathbf{x}_{k})$	relaxation/acceleration parameters, stopping tolerance
Krylov subspace methods	determines $\mathbf{r}_{0}$ and $\mathcal{K}_{k}(A,\mathbf{r}_{0})$	e . g ., CG for SPD systems, GMRES for non-SPD systems	restart length, truncation or recycling parameters, and stopping tolerances

Table 1 : Core components of stationary iteration and Krylov subspace methods.

Figure 7 : Overview of core components learning ( Section 3.1 ). (a) Core components of a general iterative method: the initial guess

\mathbf{x}_{0}

, the update function

\Phi_{k}

, and parameters

\gamma_{k}

. (b) Initial guess learning : ML predicts an approximate solution that serves as the initial guess. (c) Update function learning : ML replaces or augments the classical update function. (d) Parameter learning : ML learns algorithmic parameters from problem- or iteration-dependent information.

ref.	tested problem	ML	iter. method	main motivation
[ 34 ]	Poisson eq.	FNN	Jacobi	from ML: spectral bias of neural networks
[ 84 ]	Poisson eq.	CNN	multigrid/BiCGSTAB	from ML: improve CNN predictions
[ 85 ]	Poisson eq. from incompressible flow	CNN	Jacobi	from ML: generalization issue of the CNN
[ 86 ]	Poisson eq.	CNN	GMRES	from classical iter.: reduce iterations
[ 87 ]	Poisson eq. from plasma simulations	GraphSAGE	GMRES	from classical iter.: reduce iterations
[ 88 ]	Parametrized elasticity & Biot problems	CAE, FNN	POD-2G	from classical iter.: reduce iterations

Table 2 : Comparison of existing approaches for initial guess learning ( Section 3.1.2 ). The last column summarizes our interpretation of the main motivation behind each approach.

Figure 8 : Overview of preconditioner learning ( Section 3.2 ). (a) Given a coefficient matrix $A$ , preconditioner learning seeks to construct or enhance a preconditioner $M^{-1}$ that approximates $A^{-1}$ or its action. The arrows indicate the workflow of constructing or enhancing $M^{-1}$ from $A$ and subsequently incorporating it into a classical iterative method. We summarize two main approaches and one usage strategy. (b) Learning the action of the inverse ( Section 3.2.2 ), in which ML models directly approximate the mapping induced by $A^{-1}$ . (c) Enhancing classical preconditioners ( Section 3.2.3 ), in which ML models learn components of classical preconditioners. (d) Alternating hybrid iterative methods ( Section 3.2.4 ), which interleave classical iterative methods with ML-preconditioned steps. In (d) , a representative instance is shown in which a classical stationary iteration is alternated with an ML-preconditioned Richardson iteration. More generally, both schemes in (d) can be replaced by preconditioned Krylov methods.

ref.	tested problem	$\mathcal{N}$	ML	loss function	iter. method
[ 137 ]	shallow water equations (time-dependent)	$\mathbf{d}=\mathcal{N}(\mathbf{r})$	FNN	$\left\lVert\mathcal{N}(\mathbf{r})-\mathbf{d}\right\rVert_{2}$	GCR
[ 138 ]	elliptic PDEs (different resolutions)	$\mathbf{d}=\mathcal{N}(\mathbf{r})$	SNO	$\left\lVert\mathcal{N}(\mathbf{r})-\mathbf{d}\right\rVert_{A}/\left\lVert\mathbf{d}\right\rVert_{A}$	FCG
[ 139 , 140 ]	parametric Helmholtz equations [ 139 ] , heterogeneous fluid PDEs [ 140 ]	$\mathbf{x}=\mathcal{N}(\mathbf{b},\mu)$	U-Net	$\left\lVert A\mathcal{N}(\mathbf{b},\mu)-\mathbf{b}\right\rVert_{2}^{2}/\left\lVert\mathbf{b}\right\rVert_{2}^{2}$	FGMRES
[ 134 ]	Poisson equations from incompressible flow	$\mathbf{d}=\mathcal{N}(\mathbf{r})$	CNN	$\left\lVert\alpha A\mathcal{N}(\mathbf{r})-\mathbf{r}\right\rVert_{2}$	PSDO
[ 136 ]	Poisson equations with mixed boundary conditions	$\mathbf{d}=\mathcal{N}(\mathbf{r},\mathcal{I})$	CNN	$\left\lVert A\mathcal{N}(\mathbf{r},\mathcal{I})-\mathbf{r}\right\rVert_{2}$	PSDO
[ 135 ]	Poisson equations from incompressible flow	$\mathbf{d}=\mathcal{N}(\mathbf{r})$	DeepONet	$\left\lVert\alpha A\mathcal{N}(\mathbf{r})-\mathbf{r}\right\rVert_{2}$	Richardson

Table 3 : Comparison of several existing approaches for learning the action of the inverse ( Section 3.2.2 ). Here,

\mu

denotes the problem parameters, and

\mathcal{I}

denotes the geometric information. For [ 134 , 135 ] ,

\alpha

denotes the iteration step size. For [ 134 , 136 ] , PSDO stands for “preconditioned steepest descent with orthogonalization.” For [ 133 ] ,

\sin\left(\angle(\mathbf{r},\mathcal{W})\right)

denotes the sine of the principal angle between the residual

\mathbf{r}

and

\mathcal{W}

, the image of the Krylov subspace.

	ref.	ML model	final prec.	loss function	iter. method
IC	[ 147 ]	$L=\mathrm{GNN}(A,b)$	$M=LL^{T}$	$\\|LL^{T}\mathbf{x}-\mathbf{b}\\|_{2}^{2}$	PCG
	[ 52 ]	$L=\mathrm{GNN}(A)$	$M=LL^{T}$	$\\|LL^{T}\mathbf{w}-A\mathbf{w}\\|_{2}^{2}$	PCG
	[ 148 ]	$L_{\mathrm{IC}}=\mathrm{IC}(A)$ , $L_{\text{corr}}=\mathrm{GNN}(A)$	$M=LL^{T}$ , $L=L_{\mathrm{IC}}+L_{\text{corr}}$	$\\|LL^{T}\mathbf{w}-A\mathbf{w}\\|_{2}^{2}$	PCG
	[ 131 ]	$L_{\mathrm{IC}}=\mathrm{IC}(A)$ , $L_{\text{corr}}=\mathrm{GNN}(L_{\mathrm{IC}})$	$M=LL^{T}$ , $L=L_{\mathrm{IC}}+\alpha L_{\text{corr}}$	$\\|LL^{T}\mathbf{x}-\mathbf{b}\\|_{2}^{2}$	PCG
ILU	[ 149 ]	$\operatorname{S}_{L,U}$ = $\mathrm{U}\text{-}\mathrm{Net}(A)$	$M=LU$ , $(L,U)=\mathrm{ILU}(A;\operatorname{S}_{L,U})$	binary cross-entropy	prec. GMRES
	[ 150 ]	$(L,U)=\mathrm{GNN}(A)$	$M=LU$	$\\|LU\mathbf{x}-\mathbf{b}\\|_{2}$	prec. LGMRES

Table 4 : Comparison of existing approaches for enhancing classical preconditioners ( Section 3.2.3 ). For IC-based preconditioners,

L_{\mathrm{IC}}

denotes the factor obtained from classical IC, and

L_{\text{corr}}

denotes a learned correction term. For ILU-based preconditioners,

\operatorname{S}_{L,U}

denotes the sparsity patterns of

L

and

U

. For SPAI-based preconditioners,

T

and

D

denote a strictly lower triangular matrix and a diagonal matrix, respectively, and

\kappa(\cdot)

denotes the condition number. The vector

\mathbf{w}

denotes a random vector used in loss functions.

ref.	tested problem	classical iter.	ML-preconditioned iter.	learned prec. type
[ 37 ]	Poisson and Helmholtz eq.	Jacobi, Gauss-Seidel	DeepONet-prec. Richardson	learning the action of $A^{-1}$
[ 178 ]	Darcy flow and linear elasticity	Gauss-Seidel	DeepONet-prec. Richardson	learning the action of $A^{-1}$
[ 179 ]	Helmholtz eq.	Gauss-Seidel, GMRES	DeepONet-prec. Richardson	learning the action of $A^{-1}$
[ 135 ]	Poisson eq. from incompressible flow	CG or PCG	DeepONet-prec. Richardson (DLSM)	learning the action of $A^{-1}$
[ 180 ]	diffusion and Helmholtz eq.	Jacobi	DeepONet-prec. Richardson/Krylov	enhancing class. (subspace corr.)
[ 181 ]	Helmholtz scattering	Stationary methods	DeepONet-prec. Richardson	enhancing class. (subspace corr.)

Table 5 : Comparison of existing approaches for alternating hybrid iterative methods ( Section 3.2.4 ). The last column shows the type of learned preconditioner: learning the action of

A^{-1}

( Section 3.2.2 ) or enhancing classical preconditioners ( Section 3.2.3 ). For the latter, the underlying classical preconditioner is given in parentheses. In [ 146 ] , in addition to Jacobi and ML-preconditioned iterations, a physics-aware Anderson acceleration (PA-AA) strategy is introduced.

	multigrid	domain decomposition	multiscale	deflation/augmentation
Sec. 3.3.2	smoothers	subdomain solvers	–	–
Sec. 3.3.2	coarse-grid solvers	coarse-grid solvers	–	–
Sec. 3.3.3	coarsening in AMG	–	–	–
Sec. 3.3.4	transfer operators / coarse spaces / deflation (or augmentation) subspaces
Sec. 3.3.5	parameters	parameters	✗	✗

Table 6 : Outline of topics in Section 3.3 on enhancements of subspace correction methods . The symbol “–” indicates that the component is not applicable to the corresponding method, and “✗” indicates that we found no related work on ML enhancement in the reviewed literature.

Panel of Figure 9 — Figure 9 : Subspace correction structures in (a) MG and (b) DD methods. MG methods combine fine-grid smoothing with coarse-grid correction, whereas DD methods combine local subdomain corrections with a coarse space correction. Both approaches accelerate convergence through corrections on suitably constructed subspaces.

Figure 10 : Schematic of automated solver selection based on numerical and structural features. The general workflow consists of two stages: feature extraction and classification. The upper and lower paths can be used independently or in combination. In the upper path, selected numerical features, such as the rank and the number of nonzero entries, are fed into a classifier, typically implemented using classical ML methods such as ADT and SVM. In the lower path, CNNs or GNNs are used to extract structural features from image- or graph-based representations of the input matrix, and an FNN performs the classification. When the two paths are combined, the extracted numerical and structural features are fused into a unified representation and fed into a classifier for solver selection.

ref.	features		classifier		solver search space (# of candidates used for different problems)
ref.	numerical (#)	structural	numerical	structural
[ 290 ]	not specified	✗	ADT, SVM	✗	Krylov methods + preconditioners (48, 242)
[ 294 ]	$57\xrightarrow{\text{reduced}}27$	✗	ADT, SVM	✗	Krylov methods + preconditioners (80)
[ 295 ]	$66\xrightarrow{\text{reduced}}8$	✗	RF, XGBoost, GBDT	✗	Krylov solver + restart parameter (6)
[ 175 ]	35	✗	FNNs	✗	ILUTP_Mem parameters (72)
[ 177 ]	32	✗	reinforcement learning	✗	preconditioned solvers (576)

Table 7 : Comparison of existing approaches for automated solver selection based on numerical and structural features ( Section 3.4.2 ).

Appendix figures & tables6 assets

Supplementary material from the paper’s appendix.

Appendix

package name	description	download link
PETSc	the portable, extensible toolkit for scientific computation	https://petsc.org/release/
Trilinos	collection of open-source packages for scientific application development	https://trilinos.github.io/
Hypre	high-performance preconditioners	https://www.llnl.gov/casc/hypre/
Ginkgo	high-performance numerical linear algebra library	https://gko-project.org/
SuiteSparse Matrix Collection	formerly known as the University of Florida Sparse Matrix Collection	https://sparse.tamu.edu/

Appendix

Table 8 : Packages on classical sparse iterative solvers and sparse matrix datasets.

ref.	description	download link
[ 84 ]	Poisson CNN, learning initial guesses for multigrid & BiCGSTAB	https://github.com/aligirayhanozbay/poisson_CNN
[ 86 ]	learning initial guesses for GMRES	https://github.com/ML4FnP/GMRES-Learning
[ 89 ]	NOWS, neural operator warm starts	https://github.com/eshaghi-ms/NOWS
[ 90 ]	a pretraining-finetuning computational framework for material homogenization	https://github.com/yizheng-wang/HomoGenius
[ 91 ]	a pretraining and warm-start framework for PDEs	https://github.com/yizheng-wang/PFEM
[ 97 ]	enhancing classical iterative methods for PDEs	https://github.com/ermongroup/Neural-PDE-Solver

Appendix

Table 9 : Open-source code for core components learning ( Section 3.1 ), including initial guess learning , update function learning , and parameter learning .

ref.	description	download link
[ 137 ]	machine-learned preconditioners for GCR	https://github.com/JanAckmann/MLPrecon
[ 138 ]	FCG-NO, neural operator-based preconditioner for FCG	https://github.com/arudikov/FCG-NO
[ 130 ]	GNP, Graph neural preconditioners for FGMRES	https://github.com/jiechenjiechen/GNP
[ 134 ]	DCDM, deep conjugate direction method	https://github.com/ayano721/2023_DCDM
[ 136 ]	neural-preconditioned steepest descent with orthogonalization	https://github.com/kai-lan/MLPCG
[ 135 ]	DLSM and HyDEA, alternating hybrid iterative method	https://github.com/HMB9666/HyDEA

Appendix

Table 10 : Open-source code for preconditioner learning ( Section 3.2 ).

ref.	description	download link
[ 208 ]	learning multigrid smoothers	https://github.com/jerryhuangru/Learning-optimal-multigrid-smoothers-via-neural-networks
[ 211 ]	multigrid-augmented deep learning preconditioners	https://github.com/BGUCompSci/HelmholtzAccCNN.jl
[ 214 ]	multigrid-augmented deep learning preconditioners using compact implicit layers	https://github.com/BGUCompSci/CompactImplicitHelmholtz
[ 219 ]	algebraic multigrid coarsening using reinforcement learning	https://github.com/compdyn/rl_grid_coarsen
[ 225 ]	AMG-ANN, learning the strong threshold parameter using CNNs	https://github.com/MatteoCaldana/AMG-ANN
[ 227 ]	AutoAMG( $\theta$ ), learning the strong threshold parameter using GNNs	https://github.com/zhf-0/autoamg

Appendix

Table 11 : Open-source code for enhancements of subspace correction methods ( Section 3.3 ).

ref.	description	download link
[ 102 ]	PCGBandit, learning preconditioners for PCG using online learning	https://github.com/mkhodak/PCGBandit
[ 310 ]	automated solver selection for simulation of multiphysics processes in porous media	https://github.com/Yuriyzabegaev/solver_selector
[ 424 ]	embedding-based methods for linear solver performance prediction	https://gitlab.com/h.liu/embedding-based-predictor

Appendix

Table 12 : Open-source code for automated solver selection ( Section 3.4 ).

ref.	description	download link
[ 334 , 341 ]	hybrid Newton’s method	https://github.com/PINN-Well-opening-and-closing-events/Yads
[ 346 ]	DeepPhysics: a physics aware deep learning framework for real-time simulation	https://github.com/AlbanOdot/DeepPhysics-article
[ 347 ]	TherINO: Thermodynamically-informed iterative neural operators	https://github.com/conlain-k/therino
[ 351 ]	Newton informed neural operator for computing multiple solutions of nonlinear PDEs	https://github.com/xlliu2017/Newton-Informed-Neural-Operator
[ 355 ]	adaptive learning acceleration for nonlinear PDE solvers	https://github.com/vlssanonymous/mlsolveracc
[ 111 ]	improving pseudo-time stepping convergence for CFD simulations with neural networks	https://github.com/searhein/nn-pseudo-time-stepping-data

Appendix

Table 13 : Open-source code for ML-enhanced nonlinear iterative methods ( Section 4 ).

Explore similar work

Learning universal approximations for partial differential equations with Physics-Informed Broad Learning System
Jun 18, 2026Zhiwen Yu, Derong Yang, Liujian Zhang +5Partial Differential Equations Physics-Informed Learning
School of Computer Science and Engineering, South China University of Technology, Guangzhou 510006, China · Peng Cheng Laboratory, Shenzhen 518000, China · School of Future Technology, South China University of Technology, Guangzhou 510006, China +2
Towards Universal Convergence of Backward Error in Linear System Solvers
Apr 17, 2026Michał Dereziński, Yuji Nakatsukasa, Elizaveta RebrovaSolver Iterations Convergence
University of Michigan, Ann Arbor, MI · University of Oxford, Oxford, UK · Princeton University, Princeton, NJ
Graph Neural Multilevel Preconditioners for Iterative Solvers
Jul 30, 2026Zechen Zhang, Rui Peng Li, Yousef SaadAlgebraic Multigrid Spectral Preconditioning
University of Minnesota Minneapolis, MN, USA · CASA, Lawrence Livermore National Laboratory Livermore, CA, USA

Jun 18, 2026, cs.LGJ/K move · Enter open · S save

Learning universal approximations for partial differential equations with Physics-Informed Broad Learning System

Zhiwen Yu, Derong Yang, Liujian Zhang, Kaixiang Yang+4

School of Computer Science and Engineering, South China University of Technology, Guangzhou 510006, China · Peng Cheng Laboratory, Shenzhen 518000, China · School of Future Technology, South China University of Technology, Guangzhou 510006, China+2

Partial Differential Equations Physics-Informed Learning Least Squares Approximation

Open paper AskPDF arXiv

Partial differential equations (PDEs) play a central role in modeling complex physical, biological, and engineering systems. While traditional numerical solvers are robust, they often incur prohibitive computational costs due to mesh dependencies, whereas recent Physics-Informed Neural Networks (PINNs) offer a mesh-free alternative but frequently suffer from slow convergence and optimization instability. To bridge this gap, this article proposes the Physics-Informed Broad Learning System (PIBLS), a novel backpropagation-free framework that reformulates PDE solving as a direct least-squares optimization. We improved an algorithm within this framework to handle nonlinear PDEs efficiently and provide a rigorous mathematical proof establishing the universal approximation property of PIBLS for these equations. Experiments on linear and nonlinear PDEs demonstrate that PIBLS is one to three orders of magnitude faster than conventional PINNs while achieving significantly higher solution accuracy. This framework provides a computationally efficient paradigm for scientific machine learning, offering a practical, high-speed alternative for real-time simulation and design optimization tasks.

Abstract

Explore similar work

Learning universal approximations for partial differential equations with Physics-Informed Broad Learning System

Towards Universal Convergence of Backward Error in Linear System Solvers

Graph Neural Multilevel Preconditioners for Iterative Solvers