Neural Krylov Iteration for Accelerating Linear System Solving
Authors: Jian Luo
Jie Wang
Hong Wang
huanshuo dong
Zijie Geng
Hanzhu Chen
Yufei Kuang
Venue: NeurIPS 2024 spotlight (2024)
Keywords: AI4PDE
Neural Operator
Scientific Computing
Krylov Subspace
Primary area: machine_learning_for_other_sciences_and_fields

Abstract: Solving large-scale sparse linear systems is essential in fields like mathematics, science, and engineering. Traditional numerical solvers, mainly based on the Krylov subspace iteration algorithm, suffer from the low-efficiency problem, which primarily arises from the less-than-ideal iteration. To tackle this problem, we propose a novel method, namely **Neur**al **K**rylov **It**era**t**ion (**NeurKItt**), for accelerating linear system solving.
Specifically, NeurKItt employs a neural operator to predict the invariant subspace of the linear system and then leverages the predicted subspace to accelerate linear system solving. To enhance the subspace prediction accuracy, we utilize QR decomposition for the neural operator outputs and introduce a novel projection loss function for training. NeurKItt benefits the solving by using the predicted subspace to guide the iteration process, significantly reducing the number of iterations.
We provide extensive experiments and comprehensive theoretical analyses to demonstrate the feasibility and efficiency of NeurKItt. In our main experiments, NeurKItt accelerates the solving of linear systems across various settings and datasets, achieving up to a 5.5× speedup in computation time and a 16.1× speedup in the number of iterations.

Neural Krylov Iteration for Accelerating Linear System Solving
Jian Luo! Jie Wang!" Hong Wang'
Huanshuo Dong'
Zijie Geng' Hanzhu Chen' Yufei Kuang'
MoB Key Laboratory of Brain-inspired Intelligent Perception and Cognition, University of Science and Technology of China (jianluo, (onghonng700)@meee .usc.. edu. cn jiewangx@usto. edu.
Abstract
Solving large-scale sparse linear systems essential in fields like mathematics, sci- ence, and engineering. Traditional numerical solvers, mainly based on the Krylov subspace iteration algorithm, suffer from the low-efficiency problem, which primar- ily arises from the fess-than-ideal iteration. To tackle this problem, we propose a novel method, namely Neural Krylov Iteration (NeurKItt), for accelerating linear system solving. Specifically, NeurKItt employs neural operator to predict the invariant subspace of the linear system and then leverages the predicted subspace to accelerate linear system solving. To enhance the subspace prediction accuracy, we utilize QR decomposition for the neural operator outputs and introduce novel projection loss function for training. NeurKItt benefits the solving by using the predicted subspace to guide the iteration process, significantly reducing the number of iterations. We provide extensive experiments and comprehensive theoretical analyses to demonstrate the feasibility and efficiency of NeurKItt. our main experiments, NeurKItt accelerates the solving of linear systems across various settings and datasets, achieving up to 5.5x speedup in computation time and a 16. 1x speedup in the number of iterations.
Introduction
Solving linear systems the cornerstone of scientific computing, with applications in various fields including mathematics, science, and engineering(28) Traditional iteration algorithm tackle large-scale sparse linear systems[59) starts with a random initial vector, which corresponds one-dimensional linear subspace, and progressively expands this subspace to approximate the solution by iteratively minimizing the residual error.
Tradtioon solvers rely on the Krylov subspace
However, the convergence speed and stability of the Krylov subspace iteration algorithms are significantly influenced by the dimensions and characteristics of the employed subspaces in the iterations 3]. Computational inefficiency and instability often arise in scenarios involving high- dimensional problems large matrices, especially when the matrices exhibit poor conditioning 9 This inefficiency primarily arises from the less-than-ideal iteration, which leads to a higher number of iterations and consequently longer solving time. To address these challenges, our key insight that prior knowledge about the linear system's 's invariant subspace benefits the Krylov subspace iteration by guiding the iteration process, which reduces the required number of iterations [3]
Motivated by this insight, we introduce Neural Krylov Iteration (NeurKltt), a novel method that leverages neural networks to accelerate linear system solving. NeurKItt comprises modules:
"Correspondngg author. Email: jiewangx @ustc.edu.cn
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

BOOK
Thlerance namm
10
Tolerannce 12 (oom)
OIAAL CIPAS C10P5 C1MB55
CHELL BAOK
10-1
10-11
5236
10030
20020
33600
35.00
Time
Average Iterations
Figure I: The variation tolerance for NeurKItt compared to GMRES, where each line represents an experiment under solving method and specific preconditioning. Notably, the NeurKItt substantially enhances the efficiency of solving the linear systems, with reduction in the number of iterations by up to factor of 16 and achieving speed-up of up to 5.5 times.
the subspace prediction module and the acceleration module. The subspace prediction module uses a neural network to predict the invariant subspace of linear systems. Inspired by the concept that mapping linear system to its invariant subspace can be viewed an operator between two Hilbert spaces, we adopt the neural operator for subspace prediction. We also integrate thin QR decomposition and a projection loss function optimize training, improving subspace prediction performance. For the acceleration module, it leverages the property that Krylov subspace iteration approximates invariant subspace of the linear system. When partial information about this subspace is provided, our acceleration module utilizes it to reduce the iterations needed for precise solutions, effectively accelerating the process and addressing the challenges.
We provide comprehensive analyses to demonstrate NeurKItt's efficiency. Furthermore, extensive experiments conducted across various solver settings and different PDE problems validate the effectiveness of our approach. The results show that NeurKItt significantly accelerates linear systems solving, achieving up to 5.5x speedup. Both our theoretical analyses and experimental results collectively demonstrate the efficiency of KItt
We summarize our contributions as follows:
To the best of our knowledge, our work is the first to apply the data-driven approach to optimize Krylov iteration algorithms for solving generic non-symmetric linear systems.
We introduce a novel strategy that predicts the invariant subspace the linear system to accelerate Krylov iteration. To facilitate the subspace prediction, we design projection loss for efficient training, in conjunction with QR decomposition for stable outputs.
Extensive experiments and theoretical analysis demonstrate that NeurKItt reduces the computational cost and the number of iterations for solving linear systems involving non- symmetric matrices.
2 Related Work
2.1 Traditional Numerical Algorithms
In the field of computational mathematics, various algorithms have been devised to address the challenge solving systems of linear equations. Among these, methods based on the Krylov subspace have garnered attention for their efficacy in handling large matrices 33 43). These methods notably mitigate computational complexity by seeking approximate solutions within constrained, smaller subspace. Within this context, the Generalized Minimal Residual (GMRES) method [4] plays a pivotal role, particularly in addressing sonmyymmet matrices.
To enhance the efficiency of these Krylov subspace methods, range of preconditioning methods are employed[].. These methods aim to improve matrix conditioning, reduce iteration counts, and boost stability. Preconditioning varicties include Jacobi(46), Additive Schwarz Method (ASM)(64), and

---

Successive Over-Relaxation (2) They simplify the matrix structure for faster resolution, significant advantage large-scale problems.
Our methodology, NturKItt, utilizes invariant subspaces to hasten convergence in Krylov subspaces, thereby expediting linear system resolutions. It is also engineered to integrate with existing precondi- tioning techniques, boosting its efficiency and applicability in complex computational scenarios.
2.2 Learning-based Acceleration Algorithms
Although learning-based linear system solvers require training, the training time is negligible com- pared to the time saved over millions of calls. Thus, there has been a surge learning-based acceleration efforts in recent years. Research in using neural networks for accelerating linear system solving [30 23, 60) use neural networks to optimize the Conjugate Gradient algorithm, thereby accelerating the solution of symmetric positive definite lincar systems. [37] accelerated the solution of the Poisson equation. [20] focused on accelerating algorithm iterations by learning better algorithm parameters. Recent studies have utilized neural networks for matrix preconditioning. [17 38 53 have focused on improving algebraic multigrid preconditioning algorithms. [1. [5 have applied CNNs and machine learning, respectively, to optimize block Jacobi and ILU precondition. Furthermore, specialized preconditioning research is being conducted in areas like physics 77 engineering 47
and climate science [1], demonstrating the breadth of these applications.
However, these studies often face limitations due to specific matrix properties or mainly focus on reducing precisiciinn solving costs, with fewer advancements in accelerating high-precision solving. In contrast, our approach, NeurKItt, addresses all these limitations. It employs subspace prediction to refine iterative methods universally across Krylov subspace algorithms. This approach not only reduces the initial iterations but also significantly improves the speed and stability of subsequent iterations.
2.3 Neural Operators
Numerical methods for solving PDEs commonly involves solving systems of linear equations. which require extensive computational resources. Recently, Neural operators (NOs), such as the Fourier Neural Operator (FNO) 31 55 12 and Deep Operator Network (DeepONet) 35 36), have shown effectiveness in solving PDEs. Despite their effectiveness, NOs often grapple with challenges like diminished accuracy and stringent boundary condition prerequisites, which limit their applicability as standalone replacements for conventional algorithms in various scientific computing contexts [63] 19). Given that mapping a matrix to its subspace fundamentally an operator mapping task, NOs hold significant potential as tools for predicting matrix subspaces.
3 Preliminaries
3.1 Krylov Subspace Iteration
In the realm of large-scale sparse linear systems, Krylov subspace methods are frequently employed standard solution strategy [46] 16. At the heart of this methodology is the principle of utilizing the matrix involved in the linear system to iteratively construct subspace. This subspace is crafted to approximate the actual solution space, with the iterative process continuing until the generated Krylov subspace sufficiently encompasses the real solution. Such an approach enables the iterative solutions within this confined subspace to progressively approximate the accurate solution vector æ, thereby efficiently converging towards it.
Suppose we now solve the system
Ax = b, where A C"xn, The m-th Krylov subspace associated with the matrix A and the starting vector r e C" as follows:
(I)
Kmm (A,r) span Ar, A2r, A A r-).
Generally, r Am - represents the initial residual, where ap is typically generated randomly or selected as initial guess for the solution. The iterative process of the Krylov subspace leads to the Arnoldi relation:
AVm
V VmHm + Um h e
(3)
3

---

where V. (VI...., vm) E Caxm comprises unit vectors ei C for l...., m. mutually orthogonal, and V = (V vm+1) extends this set. H m E C(m++))xm is an upper Hessenberg matrix. H E Cm is the first m rows of Hm distinct due to an extra element A 1.m at the m + 1-th row, m-th column. C. is the m-th unit column vector, with its m-th element as 1.
Krylov algorithms reduce the computational effort needed for large linear systems by creating Krylov subspace. However, starting this process with a random vector can be time-consuming It's suggested that understanding the matrix's invariant subspace before beginning the iterative process, a warm-starting approach, could notably decrease the required iterations and improve convergence speed and stability.
3.2 Neural Operators
NeurKItt utilizes the framework of neural operator architectures, notably the FNO 26 31 27
We
begin by defining a spatial dimension E N and corresponding domain D C Fd Our focus is on approximating operators G A(D; Rd, - L ((D, 1 Here, a E A( and E U (D; (") denote functions that map the domain D to TR and TRd", respectively, where da du E N. Both spaces A(D; ) and (D; Rd=) are identified Banach spaces.
In line with the definition provided by [32], neural operator W A(D,RA.) 7 R(); is conceptualized as a mapping expressed by
N(a) = Q ... ol R(a),
for a designated depth L E N. In this formulation, R A(D; Rd) + U(D;Rd"), with de 2 du. functions as lifting operator, while Q L(D, R - U(D Rdv) serves as a local projection operator.
4 Method
Ne [utilizes neural operator to predict the invariant subspace K of the linear system, and then uses it to accelerate the Krylov subspace iteration. Generally, NeurKItt broadly divided into two components.
Subspace Prediction Module:
Inspired by operator learning, employ the neural operator capable of predicting subspaces at a low cost, with the help of QR dccomposition. The aim is to precisely predict the invariant subspace Æ of the linear system.
Originating from the mathematical theories of Krylov iteration, we have developed a Krylov algorithm tailored to our problem. Utilizing the invariant subspace, this method expedites Krylov Subspace convergence and enhances stability.
Acceleration Module:
4.1 Subspace Prediction Module
The mapping from the linear system to its corresponding invariant subspace can be considered an operator, i.e., mapping between two Hilbert spaces. Solving such problems involves finding the corresponding solution operator. In our subspace prediction module, we choose the Fourier Neural Operator (F1 for subspace prediction. We give brief introduction to FNO in Appendix C Generally, for a linear system Az derived from the parametric PDE problem, to predict its invariant subspace K. the input to FNO is the input function a E Rd from the given PDE, where d, N. We provide a detailed discussion in Appendix about how to build linear system problem from a PDE problem, and what is the input function.
Our task is to learn the mapping between two Hilbert spaces G R -> Cdx using FNO. For FNO, the lifting transformation R first lifts the input a to a higher dimensional representation by E Cde xe. where e is the number of channels. Then we feed the uo to Fourier layers. After T Fourier layers forward, we have Vr E Cd x from the last Fourier layer, which keeps the same shape as up. The FNO's output X (er) is the projection of or by the transformation Q Cd - ddxxn. Nt then uses QR decomposition to orthogonalize the matrix X, obtaining the predicted
FOur code is available https: h/tp://///tttb//////

---

Problem
Krylov Subspace Iteration
(a)
Heladadiz
(b)
Solution
A, b
A
K₂
K3
(c1)
(c2)
K
Solution
K
A,
A,b
K
K1
NeurKItt
subspace
Neural Operator
Krylov Subspace Iteration
Figure 2: Algorithm Flow Diagram: (a) Finding solution æ for given matrices A and b. (b) Traditional Algorithm: Krylov iterations from random initial vector. (cl) NeurKItt Subspace Prediction Module: Utilizing neural operator for estimating the invariant subspace of matrix A. (c2) NeurKItt Acceleration Module: Using the predicted invariant subspace of the matrix to guide the iteration, thereby accelerating the Krylov iteration process.
subspace I.. We provide more details about how to predict the subspace given the 2D Darcy flow problem in Appendix C
We define the one sided distance 51 from the subspace Q to the subspace C as
(Q2.) || Mell2
where II represents the projection operator for associated space. Its mathematical interpretation corresponds to the largest principal angle between the two subspaces Q and C(3), and defining Po as spectral projector onto Q. The objective the subspace prediction task is to approximate the invariant subspace of the matrix A with our predicted subspace K:
arg min (KC,
where 0 represents the parameters of NO, and o(KC S) is the distance between the two subspaces. Based on subsequent Theoretical Analysis 5.2 choose S the invariant subspace associated with the smallest n eigenvalues, assuming S span S span 82,... 8n where 8; is the eigenvector corresponding the i-th smallest eigenvalue of the given matrix A, 1,2..... ... n. To reduce computational complexity, considering the norm equivalence theorem in Banach spaces, we optimize using the following loss function:
arg min
d(C st).
This norm transforms the computation of subspace distance into the computation of the distance from a vector to a subspace. According to the properties of the Hilbert space, there exists a unique orthogonal decomposition:
= x Ex K, R.
(8)
We represent the projection operator for R by Q, with Pe QQ", leading to the relationship: a(kk, s.) ||s, QQ*QQQ*'''00
Consequently, our projection loss function I(K,, S) is defined as follows:
n
I(C, S) Σ ||st QQ* s:ll.
(10)
We train our subspace module by optimizing the above formula. Experiments have shown that as the projection loss decreases, the principal angle our predicted subspace also decreases. The principal angle is mathematical indicator derived in Theoretical Analysis 5.1 that influences the acceleration.
5

---

4.2 Acceleration Module
For Krylov algorithms,, the iteration count is critical factor influencing computational load. Drawing inspiration from the Krylov recycling algorithms 40 the predicted invariant subspace IC as the deflation space within the Krylov iteration framework. We provide the pseudocode the acceleration module in Appendix A.2 The key procedure of our implementation is outlined as follows:
13
6 58
41 24], our approach NeurKItt utilizes
Considering the linear system I we obtain a k-dimensional invariant subspace K. from the subspace prediction module. Then, NeurKItt computes matrices UR, Cje E C from ÆC and A such that AU = Cs and C C2 = Is. the specific algorithm can be found Appendix A.1 Arnoldi relation of our NeurKltt:
We can get the
(I-C )VVm V-iI==VVAAA Vm-k-1HH--
(11)
This suggests that when the k-dimensional invariant subspace is given, there is no need to start from scratch for constructing a new Krylov subspace K((Arr). r). Building upon this foundation, r((Arr) converge more rapidly the subspace where the solution lies. This significantly reduce the dimensionality of the final Krylov subspace, leading to a marked decrease in the number of iterations and resulting in accelerated performance. Compared to NeurKItt, GMRES can be intuitively conceptualized as the special case of our NeurKitt, where k is initialized at zero
Existing Krylov recycling algorithms share the same Arnoldi relation with NeurKItt. However, they can only accelerate the linear systems solving when there are strongly correlated linear systems provided, which is not common in scientific computing. For example, when solving PDEs with finite element software, matrices vary in size and element positions due to discretization methods and different grids, which makes the recycling invariant subspace from previous matrices impossible. Additionally, storing subspaces from previous solving requires significant storage when the previous linear system is large. Thus, these methods not widely used in general scenarios. In contrast, NeurKItt, which leverages neural operators to predict subspaces, is not affected by the lack of strongly correlated linear systems.
Theoretical Analysis
5.1 Convergence Analysis
Let Q be an I-dimensional invariant subspace of matrix A, and let C range(Css) be a k-dimensional space (k 2 1) selected to approximate Q. We refer to Theorem 3.1 [40 for an in-depth convergence analysis under NeurKItt.
Theorem 5.1. 140] Considering space C, define V = range dimensional Krylov subspace generated by NeurKli as in Let ro C², and rn = -I IIc)ro-
mm))) be the (m k)-
Then, for each Q such that 8(2,C) <1,
|Iro dil12 dievec
min
min (eE(I-Po))
(I Pelri del12
(12)
à Pellz ||(I 11yrril2.
where = - Pell and 5 = 8(2,C).
The left-hand side is the residual norm subsequent to m iterations of NeurKItt, employing the predictive subspace C. In contrast, the right-hand side, the initial term epitomizes the convergence of deflated problem, given that all components within the subspace Q have been eradicated 39 49]. The subsequent term the right embodies constant multiplied by the residual following m - iterations NeurKItt, when solving for Tı. Supposed that the predictive space encompass invariant subspace Q, then have 8 = 7 0 for the given Q. which ensures that the convergence rate NeurKItt matches or surpasses that of the deflated problem. In most cases, |Polla is numerically stable and not large, therefore reduced value of directly correlates with faster convergence in NttrK
Now compare two approaches to accelerating the linear systems solving 33

---

1. Providing an initial prediction for the solution, allows Krylov algorithms to start from this prediction. The solution involved might come from solutions to similar systems or predictions such as via neural networks 62].
2. Predicting the matrix invariant subspace, using this approximate subspace to speed up Krylov iterations.
For the first approach, particularly with large matrices derived from PDEs, the norm || tends to be large. Our theoretical analysis in Appendix (eindiixissuggsttt that FNO's 'ss direct solution prediction will not significantly accelerate the linear system solving, possibly only reducing a few Krylov subspace iterations. NeurKItt belongs to the second approach. This erucially speeds up the convergence of Krylov subspace iterations, markedly reducing total iterations and improving stability.
5.2 Subspace Property Analysis
The crucial questions about predicted subspace arise when accelerating Krylov subspace iterations: (1) What kind of subspace should selected for acceleration? (2) Given the computational cost associated with employing neural networks, is it necessary to expend substantial computational resources to predict a highly accurate subspace?
Regarding the first question: As indicated by our convergence analysis 5.1, effectively accelerating the solution of linear systems only requires predicted invariant subspace of matrix A. A matrix has many invariant subspaces, but which of these are easier to learn by the neural network(40)
To investigate this problem, simplify the specific scenario of the linear systems problem. The following definitions and assumptions for Theorem 5.2 are from the reference [24]. Specifically, We deal with Hermitian positive definite matrix A and corresponding Hermitian perturbation E, allowing A to have the eigendecompositionn 24):
[Q A1, A2, (i)
(13)
where Q = [Q Q2Q3 an orthogonal matrix, A₁ = (
A2 and Ag are defined in similar fashion. Let the eigenvalues satisfy the following ordering
x (2)
s.. S. x <...
< 1
(8) <... (3)
We consider the change the invariant subspace range(Q: (Q1) under a symmetric perturbation of A. Let 0 (-.-) denote the largest canonical angle between two spaces. We do not require that || |||| be small, but we assume that the projection E onto the subspace range( (2)) is small. We
assume that [2] EllF I and that is small relative to x x.
(1) We also assume that
|| EllE is small relative to xx Aj large. Also, let
(1) Note that we do not need to assume that > 2
= min(M (1(2) - - - n) - 2 - (1 + > 2e,
20
o= 1
A1 tc.
μ²2
Theorem 5.2. 124) Let A be Hermitian positive definite and have the eigendecomposition given in 13), and let E, T, ps, and is be defined as above. Then there exists matrix Q1 conforming to Q such that range( Q1) is simple invariant subspace of A + E, and
(aaigee(Q1)))) 0₁ (range( Q1), range( Q1)
The specific proof can be found in 24 A similar bound applies to the perturbation of the eigenvalues associated with Q1. In the context of Theorem anTand our NeurKItt, Q1 corresponds to Q, whereas Q2 and Q3 can be chosen to fit the theorem.
This theorem shows that if changes in a matrix occur in the subspace corresponding to larger eigenvalues, the subspace associated with the smallest eigenvalues is minimally affected, as long
For non- Hermitian matrices, the eigenvalues will be sorted by comparing their modulus, such that [^x| < Dal < DAx < |Aal < ... .. S [An] where X E C, C.i= 1,2, n. is the eigenvalue of given non ao--Hemmiiann matrix A.

---

as the changes are smaller than the gap between the smallest and larger eigenvalues. To facilitate learning of matrix invariant subspaces, we use the subspace composed of eigenvectors corresponding to the smallest eigenvalues as the prediction target for our neural operator.
Regarding the second question: The analysis of the first question, along with the convergence rate proved in Theorem [.I]] suggests that investing in substantial computational resources for predicting niih-prrecsiin subspace is unnecessary. As long as the predicted subspace shows distinct correlation with the matrix invariant subspace, particularly when ó is small, significant acceleration can be achieved without resorting to high precision [3]. This principle underpins our development of a low-precision cost-effective subspace prediction framework. Although the subspace prediction may not be precise, the final algorithm has achieved remarkable acceleration results.
6 Experiment
6.1 Experiment Settings
To comprehensively evaluate the performance of NeurKltt, we conducted extensive experiments. Our analysis centers on two primary performance metrics viewed through three perspectives. These tests conducted on 3 datasets. Specifically, the three Perspectives are: (i) Matrix preconditioning tech- niques, spanning standard methods. (ii) Accuracy criteria for linear system solutions, emphasizing 8 distinct tolerances. (iii) Different matrix sizes, considering 3 variations. For more details, please
refer to the Appendix(F. F.1
Baselines. NeurKItt focuses solving linear systems that involve large sparse non-symmetric matrices. The GMRES algorithm is widely used for non-symmetric large sparse linear system solving, which serves as the predominant solution. And we set it as the benchmark for our study. We use GMRES from PETSc (version 3.19.4).
Datasets. To investigate the algorithm's 's adaptability to various types of matrices, we examined three
different linear equation challenges, each rooted in PDE: 1. Helmholtz Equation 62
2. Darcy
Flow Problem [31] 27 36]: 3. Non-uniform Heat Conduction Equation [48] 25
18]. For an
in-depth exposition of the dataset and its generation, kindly refer the Appendix For the runtime
environment, refer to Appendix F2
Metrics. We adopt time speedup for "GMRES solving time / NeurKItt solving time" and iteration speedup for "GMRES iteration count NeurKItt iteration count". Speedup over denotes better NeurKItt performance. We also provide the average time spent and average iteration count in Appendix H We use principal angle in Equation So show how close the predicted subspace and the target subspace are.
Table 1: Comparison of our NeurKItt and GMRES computation time and iterations across datasets, preconditioning, and tolerances. The first column lists datasets with matrix size, and the next details tolerances. Results are displayed as "time speedup iteration speedup".
Dataset
Tol
None
Jacobi
BJacobi
SOR
ASM
ICC
ILU
le-2
5,55
16.13
3.66/ 11.00
1.87 / 6.82
2,20 7.78
1.91 6.64
2.64/ 6.35
2,60/ 6.35
le4
4.66 12.35
7.66
1.97 5
6.24/
2.05 2.05/5.55 /5.58
2.40/
2.36/ 5.47
Darey
1e-7
39.99
2,38 5.88
1,75 /4.46
1.98 /497
1.91 /4.60
/4.6/
2.06 24.62
160000
le-10
3.67 / 8.98
2,16/
1.65 /4.08
1.81 /4.31
1883/.113
1.86/4.16
1.87 /4.16
1e-12
3.61 8881
2.22/ 5.45
1.67 3.9.7
1.72/4.05
1.74/ 3.94
1.79/3996
3.99
le-2
3.15/9.48
2.24/6.20
1.52/5.33
5.54/5.96
1
1.90/4833
1.97 /4.83
le-4
3.23/8.55
2.03 5.27
1.77 /.79
/
1.49/4.39
1.84/.3
1.89 4,30
le-7
27.63
1.55 /4.35
1.51
1.51 /4.28
1.45/3.81
3.73
1,70
90000
le-10
2.93/7.34
/.50/4.27
4.50/
L.2/3.90
1.38
3.54
lc-12
2
1.54 /4.35
1.55 3.96
1.41 3.76
1.22/3.51
1.61 / 3,46
/.46
le-2 2.39
77.68
1.33/5.49
5.56/
1.32 4.88
1.54/4.31
1.51
lc44
2.69
2.72/7.99
1.36/4299
1.45
L.0/4.000
/
1.59/
Helmholtz
le:7
2.53/7.29
2.65/7.29
1.52/4.08
1.60
1.43
1.68 / 3.76
1.68 3.76
62500
le-10
2.52/6.74
6.77
1.56
1.77 /4.59
1.56/
1,70 3.75
1.69
2.53/6.56
2.48 6.556
1.72 4..5
1.64/3.97
1.68.

---

6.2 Main Experiment
We compare NeurKItt with GMRES under various preconditioners and tolerance, and Table presents
main experimental results. More detailed experimental results about time and number of iterations can be found in Appendix [
The results show that across all tolerances and preconditioning techniques, NeurKIt has been effective in accelerating the linear system solving and reducing the number of iterations required by the Krylov subspace method, achieving up to a 5.5x speedup in computation time and a 16.1 x speedup in the number of iterations. These results demonstrate that the invariant subspace predicted by NeurKIt greatly enhances the convergence speed of the Krylov algorithms, thereby reducing the iteration count and accelerating the solving.
With the solution accuracy inereasing, the acceleration does not significantly decrease. Though the acceleration effect of NeurKItt varies under different preconditioning, the minimum time speedup is no less than 1.67 in the Darcy flow problem, which suggests NeurKItt can be effectively combined with various preconditioning methods.
6.3 Generalization of Matrix Size
None,tal
None.tol 1e-4
None.tal 1e-8
None.tol
Jacabi,tol - 1o-4
Iteration Speedup
Jacobi.tol 1e-4
Time Speedup
Jacobi.tol
Jacobi.tol
4.5
Matrix Size
Size
Figure 3: Experiments on the Darcy Flow problem with varying matrix sizes. The results indicate that as the matrix size increases, both time speedup and iteration speedup increase.
We conducted experiments the Darcy Flow problem with varying matrix sizes under different preconditioning methods and tolerances. We provide the results in Figure The experimental results show that method consistently improves its acceleration effect as the matrix size increases, across different preconditioning methods and tolerances, which indicates that NeurKItt especially benefits solving large sparse linear systems.
6.4 Ablation Study
We conducted ablation experiments to further validate the performance of FNO and demon- strate the effectiveness of the designed projec- tion loss. The results of the ablation experiments are shown in Table 2 The results demonstrate that FNO significantly outperforms MLP when using the same projection loss. This indicates that FNO indeed capable of learning better mapping from the parameter matrix to the sub- space, resulting in predicted subspaces that are closer to the target subspaces.
Table 2: Comparison of different settings on prin- cipal Angle (rad) for the Darcy flow problem, with a matrix size of 32,400. We replace the FNO with MLP for "w/o FNO" and replace the projection loss with MSE loss for "w/o Projection Loss".
Model
Principal Anglel
NeurKItt
0.07
w/o FNO
0.54
In addition, we also attempted to utilize the Mean Squared Error (MSE) Loss for subspace prediction. The experimental results show that applying the MSE Loss on subspace learning merely ineffective. Since we report the principal angle in radians, it implies that the predicted subspace IC merely orthogonal to the target subspace S. Such subspace I cannot be effectively utilized to accelerate linear system solving. Besides, can be inferred that QR decomposition we have
w/o QR decomposition
1.57
w/o Projection Loss
1.57

---

employed plays an important role in subspace prediction. These results verifies that the projection loss and QR decomposition enable the model to better predict subspace.
7 Limitation and Conclusions
Limitation Our method only considers solving non-symmetric matrices, but for matrices with specific structures, we need to tail our method to achieve faster solving. Furthermore, we have not considered the potential impact of different preconditioning techniques on the subspace and have not optimized for different preconditioning techniques, though the influence caused by preconditioning has little impact based on the observation of our experiments.
Conclusions We propose NeurKItt, a novel method for accelerating linear systems solving through subspace prediction. To the best of our knowledge, this is the first attempt to apply a data-driven approach to optimize the Krylov subspace algorithm for non-symmetric linear systems. Specifically, we employ FNO to predict the invariant subspace and design projection Ioss function for this task. By utilizing the predicted invariant subspace, we achieve accelerated Krylov subspace itera- tions. Theoretical analysis and extensive experiments demonstrate that NeurKItt effectively reduces computational costs and iteration counts.
Acknowledgements
We would like to thank all the anonymous reviewers for their insightful comments. This work was supported in part by National Key R&D Program of China under contract 20222D0119801, National Nature Science Foundations of China grants U23A20388 and 62021001.
References
[1] Jan Ackmann, Peter D Düben, Tim N Palmer, and Piotr K Smolarkiewicz. Machine-learned preconditioners for linear solvers in geophysical fluid flows. arXiv preprini arXiv:2010.02866, 2020.
[2] Yael Azulay and Eran Treister. Multigrid-augmmnttér deep learning preconditioners for the helmholtz equation. SIAM Journal on Scientific Computing, (0):S127-S151, 2022.
[3] Christopher Beattie, Mark Embree, and John Rossi. Convergence of restarted krylov subspaces to invariant subspaces. SIAM Journal on Matrix Analysis and Applications, 1574)1007-1
2004.
[4] James V Beck, Ben Blackwell, and Charles R St Clair. Inverse heat conduction: Ill-posed problems. James Beck, 1985.
[5] Michele Benzi, Carl D Meyer, and Miroslav Tma. A sparse approximate inverse preconditioner for the conjugate gradient method. SIAM Journal on Scientific Computing. (155:1135-1149) 1996.
[6] Philipp Birken, Jurjen Duintjer Tebbens, Andreas Meister, and Miroslav Tma. Preconditioner updates applied to cfd model problems. Applied Numerical Mathematics, 18(11)1162--161, 2008.
[7] Salvatore Cali, Daniel C Hackett, Yin Lin, Phiala E Shanahan, and Brian Xiao. Neural-network preconditioners for solving the dirac equation in lattice gauge theory. Physical Review D, 107(3):335508, 2023.
[8] Luiz Mariano Carvalho, Serge Gratton, Rafael Lago, and Xavier Vasseur. A flexible generalized conjugate residual method with inner orthogonalization and deflated restarting. SIAM Journal on Matrix Analysis and Applications, 3214): 2011.
[9] Ke Chen. Matrix preconditioning techniques and applications, volume 19, Cambridge Univer- sity Press, 2005.
[10] Lu Cheng and Kuan Xu. Solving time-dependent pdes with the ultraspherical spectral method, 2023.
10

---

[11] Eric De Sturler. Truncation strategies for optimal krylov subspace methods. SIAM Journal on Numerical Analysis, 36(3):864 889, 1999.
[12] Wenhan Gao, Ruichen Xu, Hong Wang, and Yi Liu. Coordinate transform fourier neural operators for symmetries in physical modelings. Transactions Machine Learning Research, 2024.
[13] André Gaul. Recycling krylov subspace methods for sequences of linear systems. 2014.
[14] Gene H Golub and Charles F Van Loan. Matrix computations. JHU press, 2013.
[15] Markus Götz and Hartwig Anzt. Machine learning-aided numerical linear algebra: Convo- lutional neural networks for the efficient preconditioner generation. In 2018 IEEE/ACM 9th Workshop on Latest Advances in Scalable Algorithms for Large-Scale Systems (scalA), pages 49-56. IEEE, 2018.
[16] Anne Greenbaum. Iterative methods for solving linear systems. SIAM, 1997.
[17] Daniel Greenfeld, Meirav Galun, Ronen Basri, Irad Yavneh, and Ron Kimmel. Learning to optimize multigrid pde solvers. In International Conference on Machine Learning, pages 2413-2423. PMLR, 2019.
[18] Stefan Güttel and Françoise Tisseur. The nonlinear eigenvalue problem. Acta Numerica, 26:1-94, 2017.
[19] Zhongkai Hao, Songming Liu, Yichi Zhang, Chengyang Ying, Yao Feng, Hang Su, and Jun Zhu. thysics-informeet machine learning: A survey on problems, methods and applications. arXiv preprint arXiv:22/1.08066 2022.
[20] Jun-Ting Hsieh, Shengjia Zhao, Stephan Eismann, Lucia Mirabella, and Stefano Ermon. Learn- ing neural pde solvers with convergence guarantees. arXiv preprint arXiv:19066 01200, 2019.
[21] Thomas JR Hughes. The finite element method: linear static and dynamic finite element analysis. Courier Corporation, 2012.
[22] Claes Johnson. Numericul solution of partial differential equations by the finite element method. Courier Corporation, 2012.
[23] Ayano Kaneda, Osman Akar, Jingyu Chen, Victoria Alicia Trevino Kala, David Hyde, and Joseph Teran. A deep conjugate direction method for iteratively solving linear systems. In International Conference on Machine Learning, pages 15720-15736. PMLR, 2023.
[24] Misha E Kilmer and Eric De Sturler. Recycling subspace information for diffuse optical tomography. SIAM Journal on Scientific Computing, 276)):140-2216 2006.
[25] Seid Korie and Diab W Abueidda. Data-driven and physics-informed deep learning operators for solution heat conduction equation with parametric heat source. International Journal of Heat and Mass Transfer, 203: 123809, 2023.
(26) Nikola Kovachki, Samuel Lanthaler, and Siddhartha Mishra. On universal approximation and error bounds for fourier neural operators. The Journal of Machine Learning Research, 12(1311323-13331. 2021.
[27] Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Learning maps between function spaces. arXiv preprint XXi::2108.08488 2021.
[28] Steven Leon, Lisette G De Pillis, and Lisette G De Pillis. Linear algebra with applications. Pearson Prentice Hall Upper Saddle River, NJ, 2006.
[29] Randall J Vequee Finite volume methods for hyperbolic problems, volume 31. Cambridge university press, 2002.
[30] Yichen Li, Peter Yichen Chen, Wojciech Matusik, et al. Learning preconditioner for conjugate gradient pde solvers. arXiv preprint 20Xv:2305.166432, 2023.

---

[31] Zongyi LI, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differen- tial equations. arXiv preprini arXiv:2010.08895, 2020.
[32] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar, Neural operator: Graph kernel network for partial differential equations. arXiv preprint (rxxi:2003.03445 2020.
[33] Jörg Liesen and Zdenek Strakos. Krylov subspace methods: principles and analysis. Numerical Mathematics and Scie, 2013.
[34] Chih-Jen Lin and Jorge Moré. Incomplete cholesky factorizations with limited memory. SIAM Journal on Scientific Computing, 21(1):24-45 1999.
[35] Lu Lu, Pengzhan Jin, and George Em Karniadakis. Deeponet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators. arXiv preprini
1910,03/99 2019.
[36] Lu Lu, Xuhui Meng, Shengze Cai, Zhiping Mao, Somdatta Goswami, Zhongqiang Zhang, and George Em Karniadakis. A comprehensive and fair comparison of two neural operators (with practical extensions) based on fair data. Computer Methods in Applied Mechanics and Engineering, 393: 14778, 2022.
[37] Kevin Luna, Katherine Klymko, and Johannes Blaschke. Accelerating gmres with deep learning in real-time. arXiv preprint arXiv:2103. 2021.
[38] Hla Luz, Meirav Galun, Haggai Maron, Ronen Basri, and Irad Yavneh. Learning algebraic multigrid using graph neural networks. International Conference on Machine Learning. pages 6489-6499. PMLR, 2020.
[39] Ronald B Morgan. Gmres with deflated restarting. SIAM Journal on Scientific Computing, 24(1):20-37, 2002.
[40] Michael L Parks, Eric De Sturler, Greg Mackey, Duane D Johnson, and Spandan Maiti. Recy- cling krylov subspaces for sequences of linear systems. SIAM Journal on Scientific Computing, 28(5)1165111674 2006.
[41] Michael Lawrence Parks. The iterative solution of a sequence of linear systems arising from nonlinear finite element analysis, 2005.
[42] William H Press. Numerical recipes 3rd edition: The art of scientific computing. Cambridge university press, 2007.
[43] Ouyuan Qin and Kuan Xu. Solving nonlinear odes with the ultraspherical spectral method. arXiv preprini arXiv:2306. 17688, 2023.
[44] Md Ashiqur Rahman, Zachary E Ross, and Kamyar Azizzadenesheli. U-no: U-shaped neural operators. arXiv preprint 220i::2004111227, 2022.
[45] Youcef Saad and Martin H Schultz. Gmres: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM Journal on scientific and statistical computing. 733):856-86, 1986.
[46] Yousef Saad. Iterative methods for sparse linear systems. SIAM, 2003.
[47] Johannes Sappl, Laurent Seiler, Matthias Harders, and Wolfgang Rauch. Deep learning of preconditioners for conjugate gradient solvers in urban water related problems. arXiv preprint 1XXi:19966009255. 2019.
[48] Rishi Sharma, Amir Barati Farimani, Joe Gomes, Peter Eastman, and Vijay Pande. Weakly- supervised deep learning of heat transport via physics informed loss. arXiv preprint 100. 2018
[49] Valeria Simoneini and Daniel Szyld. On the occurrence of superlinear convergence of exact and inexact krylov subspace methods. SIAM review, 47(2):247-77,, 2005.
12

---

[50] Rita Stanaityte. ILU and Machine Learning Based Preconditioning for the Discretized Incom- pressible Savier-Stokes Equations. PhD thesis, University of Houston, 2020.
[51] G. W. Stewart and Ji guang Sun. Matrix Perturbation Theory. Academic Press, 1990.
[52] John Strikwerda. Finite difference schemes and partial differential equations. SIAM, 2004.
[53] Ali Taghibakhshi, Scott MacLachlan Luke Olson, and Matthew West. Optimization-ban based algebraic multigrid coarsening using reinforcement learning. Advances in neural information processing systems, 12111222212140 2021.
[54] Andrea Toselli and Olof Widlund. Domain decomposition methods--lgorrthms and theory, volume 34. Springer Science & Business Media, 2004.
[55] Alasdair Tran, Alexander Mathews, Lexing Xie, and Cheng Soon Ong. Factorized fourier neural operators. In The Eleventh International Conference Learning Representations.
[56] Lloyd N Trefethen. Spectral methods in MATLAB. SIAM, 2000.
[57] Lloyd N Trefethen. Spectra and pseudospectra: the behavior of nonnormal matrices and operators. 2020.
[58] Hong Wang. Zhongkai Hao, Jie Wang, Zijie Geng, Zhen Wang, Bin Li, and Feng Wu. Ac- celerating data generation for neural operators via krylov subspace recycling, arXiv preprint arXivv2401.09516. 2024.
[59] David S Watkins. The matrix eigenvalue problem: GR and Krylov subspace methods. SIAM, 2007.
[60] Cheng Yang. Xubo Yang. and Xiangyun Xiao. Data-driven projection method fluid simulation. Computer Animation and Virtual Worlds, 27(3-)):415 424, 2016.
[61] David Young. Iterative methods for solving partial difference equations of elliptic type. Trans- actions of the American Mathematical Society, 76(1):92-111 1954.
[62] Enrui Zhang, Adar Kahana, Eli Turkel, Rishikesh Ranade, Jay Pathak, and George Em Kami- adakis. hybrid iterative numerical transferable solver (hints) for pdes based on deep operator network and relaxation methods. arXiv preprint arXiv:2208. 01222088.33733 13273, 2022.
[63] Xuan Zhang, Limei Wang, Jacob Helwig. Youzhi Luo, Cong Fu, Yaochen Xie. Meng Liu, Yuchao Lin, Zhao Xu, Keqiang Yan, Keir Adams, Maurice Weiler, Xiner Li, Tianfan Fu, Yucheng Wang. Haiyang Yu, YuQing Xie, Xiang Fu, Alex Strasser, Shenglong Xu, Yi Liu, Yuanqi Du, Alexandra Saxton, Hongyi Ling, Hannah Lawrence, Hannes Stärk, Shurui Gui, Carl Edwards. Nicholas Gao. Adriana Ladera, Tailin Wu, Elyssa F. Hofgard, Aria Mansouri Tehrani, Rui Wang. Ameya Daigavane, Montgomery Bohde, Jerry Kurtin, Qian Huang, Tuong Phung, Minkai Xu, Chaitanya K. Joshi, Simon V. Mathis, Kamyar Azizzadenesheli, Ada Fang. Alán Aspuru-Guzik, Erik Bekkers, Michael Bronstein, Marinka Zitnik, Anima Anandkumar, Stefano Ermon, Pietro Liò, Rose Yu, Stephan Günnemann, Jure Leskovee, Heng Ji, Jimeng Sun, Regina Barzilay, Tommi Jaakkola, Connor W. Coley, Xiaoning Qian, Xiaofeng Qian, Tess Smidt, and Shuiwang Ji. Artificial intelligence for science in quantum, atomistic, and continuum systems. arXiv preprini arXiv:2307.008423, 2023.
[64] Xuejun Zhang. Multilevel schwarz methods. Numerische Mathematik, 63(1):521-539) 1992.
13

---

A Algorithmic Details
A.1 computes matrices Ug and C₁
NeurKItt can be modified solve I by carrying over subspace i span .2..... : with respect to A as follows I 40]:
[Q: R thin QR decomposition of AYk.
C Q;
U = Y.
where the matrix Y2 (y y2.
yı. The matrices Us and Cs to satisfy
AU CRi
C C2 = Is-
A.2 NeurKltt Pseudocode
Algorithm I: NeurKltt Krylov subspace iteration
Select m, the maximum size of the Krylov subspace. Define as the dimension of the predicted
invariant subspace É, where i = span y/,..... y/h Given matrix use subspace prediction module to predict the matrix Y3 [y) U2.... yk and
set tol as the convergence tolerance. Choose an initial guess En. Calculate = b - A.o. and initialize
1.
Let QQ RR] be the reduced QR-factorization of AYk-
C* = Q
5 U* Y
21 = Do + +U.Ck H
ro
Pı ro C ro
8 while |r.l22 >
tol do
i=i+1
10
Perform m Amoldi steps with the linear operator (I C2 "
A, )A letting
11-1/111;-1 and generating Vm-k+1, H Hm-k and Bm-k- 81
Let D; be diagonal scaling matrix such that UR DR where the columns of Uy have unit norm.
Vn
1
Vm
13
W = |C&
Vm-k+1]
Gm Dx
Bm-k
14
Hm-k
15
Solve min [ m+11:-11 Gm for y.
x; -:-1 Vny
17
r; = = Wm-IG
Compute the eigenvectors žn of G G, Gf W m+1 Vm, associated with smallest
magnitude eigenvalues 0, and store in PR.
Y₁ V Pe
19
20
Let Q, R] the reduced QR-factorization of G,,, Pk-
C₁ W Qm+1Q
22
U = Y
23
2; as the numerical solution, r; as the residual, as the iteration count.
14

---

B From Partial Differential Equation to Linear System
We model a steady-state PDE problem as a corresponding system of linear equations using discrete numerical methods such as FDM, FEM, and FVM [22] 21 22 29, 10 Numerical methods discretize
a partial differential equation problem by mapping it from an cehnite-dimensionnn function space to fiiie-dimennional space, resulting in system of linear equations.
We take the process of discretizing the wwo-dimensional inhomogeneous Helmholtz equation using the Finite Difference Method (FDM) as an example to illustrate how to transform a PDE into a system of linear equations:
(r, v) + (x, ()) y) = [(r,y)-
(14)
This equation could be approximated as:
u(x + At, y) + ufx - St, y) + u(x, y + At) + u(I, y- At) 4u(x, v) + 1 (x, (),()) = f(x, 9);
A/2
(15)
We use a 2 x 2 internal grid (i.c., N₂ = N2 = 2 and A = Ay), the unknowns Vi,j can be arranged in row-major order as follows: U1,1, 11,, U2,1 u2.2- For central differencing on 2 grid, the vector b will contain the values of Ji, f Mj) and the linear equation system Aæ b can be expressed as:
4 + kf,1
1
11,1
f
1
-4 +
0
11,2
fi,2
1
-4+
t/21
f2,1
0
1
4 +
22,2_
[2222]
L/2.2
In this example, we call discretization function ka(x,, y) the input function, corresponding to 2 x 2 grid. The function inputs vary from PDE problems.
C A Brief Introduction to Fourier Neural Operator
The Fourier Neural Operator (FNO) aims to map continuous input function A to its corresponding continuous output function 24 within the Fourier domain. For end-to-end training on FNO, the function pairs (A, are discretized into instance pairs (a, u). The purpose is to learn a mapping N between (a, u). which can be expressed as:
u == (aa)
(16)
The mapping A consists of several sequential steps: first, the input channel is lifted using R; next, the mapping is performed through Fourier layers (L1,62,. )r.); finally, the output is projected back to the original channel using Q.
((a) = Q L1 o..... L1 R(a),
(17)
Q and RR are pixel-wise transformations that be implemented using models like MLP.
typical Fourier layer consists of a pixel-wise linear transformation, characterized by weight W and bias 6, along with integral kernel operator K::
[(x) oW + b K(x)),
(18)
where a is the nonlinear activation function. The integral kernel operator I. composes three steps: Fast Fourier Transformation (FFT), spectral linear transformation, and inverse FFT.
In the subspace prediction task, we use projection layer Q mapping the last fno layer output to the target space. Each channel in the FNO output represents one of the basic vector of the predicted invariant subspace. We then apply QR decomposition to the output. Let the input function be Rdd xd. from the Darcy flow problem, where d, E N the resolution of the input function, which yields linear system Ax b for numerical solving, where matrix A E &dd We aim to predict subspace K, as matrix X Edaxn for matrix A in the linear system. First, the FNO's
15

---

lifting layer R transforms a to a higher- lesber-dimensional representation E Cda xe, while the 3rd dimension is channel dimension. After T Fourier layers forward, we have the ey C xd, x Because Fourier layers keep the shape unchanged, we apply transformation Q to map the vr to the desired space, K Q(eT), with Q Cd.xdaxx - Cdaxn.
In practice, transformation Q is stack of transformation layers. first flattens the first and second dimension of Ur, obtaining go x Then a fully-connected neural network (FNN) applies the mapping C + x to qn- obtaining 41 E C xn. And another FNN applies the mapping C xn Cdax= X to 91, obtaining the output X E CdAXn, Finally, we apply QR decomposition to X for orthogonalizing, obtaining K = (X).
D Analyzing the Efficiency of Using FNO Predictions as Initial Solutions in Linear System Solvers
D.1 Theoretical Analysis
The central premise of this discussion hinges on the potential application of Neural Operators (NOs), such as the Fourier Neural Operator (FNO) [31], in addressing Partial Differential Equation (PDE) problems. Specifically, explores the feasibility of utilizing these neural operators to predict solutions to PDE problems. Such predictions could serve preliminary solution vectors for Krylov subspace methods in the resolution of corresponding linear systems. This approach holds promise for expediting the solution process.
However, the practicality of this concept, particularly in contexts where FNO is commonly discussed, warrants careful consideration. This is especially true in scenarios involving large matrix linear systems that emerge from Finite Difference Methods (FDM) [52], spectral methods [56), and similar approaches. A significant obstacle to the application of this method is the disparity between the type relative error typically discussed in the context of neural operators and the relative error conventionally employed in computational mathematics. This divergence is further complicated by the impact of the norm of the PDE discretization matrix, which exhibits a tendency to increase with the refinement of the grid.
In the context of linear systems derived from PDEs, we consider the following standard formulation:
Ax b.
(19)
Within the framework of NOs, the prevalent metric for error assessment is the Mean Squared Error (MSE) of the solution function solution vector æ. designated as rlol
1 x:
vo
(20)
2
where No represents the solution as predicted by the NO, while x denotes the precise solution.
In the domain of computational mathematics, the relative residual rtol is commonly utilized as the error metric for resolving linear systems:
A& - bll [
rtol=
(21)
where è represents the solution as predicted by the Krylov subspace method.
Utilizing orthogonal polynomial grid theory, specifically Theorem 30.1 in Trefethen's work [57]. page 290, offers crucial insights for this analysis.
Theorem D.1.
[5 290] For any N, Dx > N2/3.
In this context, N denotes the count of discretization points in a one-dimensional Ordinary Differ- ential Equation (ODE), and D represents the discretization matrix constructed using Chebyshev orthogonal polynomial zeros as grid points, with the matrix dimension being N. This assertion holds similarly for other orthogonal polynomial grids, and uniformly divided grid can be analogized to grid defined by the zeros of Fourier polynomials.
16

---

The theorem suggests that with an increase in the number of grid points, resulting in a larger matrix size, the norm of the matrix || associated with the discretized PDE also escalates. Then We consider the following equation:
roolno
| x
1iNo - x: ||($
(22)
rtol
A b1l
x>
Consequently, when |||l2 is substantially large, there emerges notable discrepancy between rtol and rtol. Generally, the predictive accuracy rlolxo of the FNO lies in the range 10 -0-22pp to 10 However, if the norm |All2 approaches the magnitude of 10 the relative tolerance rtol as per computational mathematics algorithms may only range between 10 and 10 such instances, leveraging the FNO's s predicted solution as the initial vector for solving the linear system may not substantially expedite the solution process.
D.2 Experimental Analysis
To empirically corroborate the theoretical framework presented, we delve into the Heat dataset F.4.3 utilized in this study, characterized by a matrix size of 52900. The experimental results underscore several noteworthy observations.
Employing the original FNO for training, we observed rtol 6 10-3, which is line with the benchmarks reported in the foundational FNO paper. However, because || 1.2 10 it results in rtol ~ 60. This empirical finding consistent with our theoretical predictions. In such context, utilizing the FNO-derived solution as the initial vector for the linear system does not markedly expedite the resolution process.
This analysis implicitly suggests that utilizzig predicted solutions as initial vectors in iterative algorithms mandates a high level of precision the neural network model. Nonetheless, it is observed that in standard application scenarios, the impact on accelerating the process is marginal.
In contrast this method, algorithm NeurKItt, by predicting the invariant subspace, not only reduces the initial number of iterations but also, more importantly, substantially improves the convergence rate in subsequent iterations, As demonstrated in the convergence analysis [5.1] This marks a pivotal deviation from former techniques that focus on expediting the process by merely predicting initial solutions. Our theoretical analysis 5.2 shows that the predicted invariant subspace does not require high precision to accelerate the solution of linear systems and significantly reduce the number of iterations.
E Thin QR Decomposition
In our research, we utilize the Thin QR decomposition algorithm, prevalent technique in numerical linear algebra [4]. The essence of Thin QR decomposition lics in its capacity to factorize a matrix into distinct components: orthogonal matrix Q and an upper triangular matrix R.
This factorization is instrumental in diminishing both the computational complexity and the storage demands.
For a matrix A belonging to Caxm with the condition n 2 m, the Thin QR decomposition is mathematically articulated as:
A = A
(23)
where Q is an n x m orthogonal matrix that complies with Q'Q = and an m x m upper triangular matrix. In study, the integration of QR decomposition within the neural network architecture strategically employed to effectively resolve the orthogonal basis of a subspace.
Details of Experimental Data and Setup
F.1 Specific parameters of the main experiment
Model Training and Predicting: We employ the 5 FNO layers with modes from (20, 32, 40) and the width from (32, 50, 64). The learning rate is fix at 10-3 while the batch size selected from 116, 32). The number of subspace dimensions is fixed at 10 in all experiments. For darcy flow and
17

---

heat datasets, we generate 8000 samples for training and 1600 samples for testing. For Helmholtz, we generate 1000 samples for training and 200 samples for testing.
Baseline: To achieve maximum efficiency, we initiated our process by creating linear equation systems of PDEs using either Python. These subsequently processed within the C programming environment. For the GMRES solver, we employed the most recent version of PETSc, specifically 3.19.4.
Three Perspectives:
1. Precondition:
The act of preconditioning is vital when tackling large matrix equations, as it significantly speeds up the resolution process and augments algorithmic stability. Given the variety of scenarios that may arise, we explored over 10 prevalent preconditioning methods, adapting them to suit different contexts.
2. Tolerance: The precision threshold of a solution directly influences both the number of itera- tions required and the overall computational time. Different algorithms manifest varying rates of convergence, and specific NOs hold unique tolerance standards. Our comprehensive evaluation encompassed a range of tolerances, with focus on pinpointing 5-8 levels of optimal error precision for our analyses.
3. Matrix Dimensionality: Furthermore, different NOs often require matrices of various sizes. Therefore, our investigation included the study of 5 distinct matrix sizes.
The efficiency of algorithm can vary with changes in matrix size.
Two Performance Metrics: 1. Computational Time: This metric provides a clear and direct
measure of an algorithm's performance. 2. Iteration Count: The count of iterations required reflects the algorithm's numerical sensitivity and stability, offering valuable insights into its performance characteristics.
F.2 Environment
To ensure consistency in our evaluations, all comparative experiments were conducted under uniform computing environments. Specifically, the environments used are detailed follows: CPU: Intel Xeon Gold 6246R CPU 3.40GHz
GPU: NVIDIA GeForce RTX 3090
F.3
Precondition
In this study, we implemented series of preconditioning methods for our experiments, which are outlined below.
1. None: In scenarios where 'none' is selected for preconditioning. the linear system is addressed directly in its original form. This means that no preliminary transformations or adjustments are applied to the system before the iterative solving process.
Diagonal Preconditioning (Jacobi) [46]: This technique focuses solely on the diagonal components of the coefficient matrix. The preconditioner in this approach is essentially the inverse of these diagonal elements, offering a straightforward preconditioning solution.
3. Block Jacobi (BJacobi) [5 A more advanced form of Jacobi preconditioning. Block Jacobi divides the coefficient matrix into smaller, manageable blocks, each representing a specific subdomain or a separate problem. These blocks are then preconditioned individually using their diagonal elements.
Successive Over-relaxation (SOR) 611 SOR, modification of the Gauss-Seidel method, incorpo- rates a relaxation factor to enhance the process of convergence. This method reformulates the original problem by applying weight, often improving the convergence rate for certain iterative algorithms.
Additive Schwarz Method (ASM) 54 ASM operates on the principle of domain decomposition, partitioning the main problem into several smaller subdomains. Each subdomain's 's problem is solved independently, and the local solutions are then aggregated to form the overall solution.
6. Incomplete Cholesky (ICC) 34 ICC is a preconditioning method based on the Cholesky decom- position, but drops certain off-diagonal elements during the decomposition, making "incomplete". It's utilized for symmetric positive definite problems.
18

---

7. Incomplete LU (ILU) [6] ILU is based on LU decomposition, but like ICC, drops certain
off-diagonal elements during the decomposition. ILU can be applied to nonsymmetric problems.
F,4
Datasets
F.4.1
Helmholtz Equation
We consider a two-dimensional Helmholtz equation, which can be described by the following equation [62]:
Au(a, y) + (x, g) y) [(x v),
(24)
where uxx, y) the wave function, and k(x;, y) denotes the spatially varying wavenumber, with representing the source term, corresponding to origins of electromagnetic or acoustic waves. This equation is pivotal in several physical realms, notably Acousties, Electromagnetism, and Quantum Mechanics. Within these fields, the Helmholtz equation models phenomena such as wave propagation and vibration patterns.
In the context of the Helmholtz equation, the wavenumber k(z, y) is intimately connected to the frequency of the wave and the characteristics of the medium through which the wave propagates. In our experimental setup, the value of k(x,, y) is determined using the GRF method [35 while the fæ, y) the constant.
F.4.2 Darcy Flow problem
We consider two-dimensional Darcy flows, which can be described by the following equation [3 44
27 36]:
(K (K Vh(r,))) [xx, (x,y),
(25)
where K is the permeability field, h is the pressure, and is a source term which be either constant or a space-dependent function.
In our experiment, K K(r,,) y) is derived using the GRF method 35 The term f(z, v) is consistently set as a constant f.
F.4.3 Non-uniform Heat Conduction Equation
We consider a two-dimensional heat conduction equation in non-uniform medium, which is formulated as follows 18
8TT 0x2
8T
0k 9T
0k OT Dy dy + qT = f.
(26)
0y²
Ox dr
where T(x, y) is the temperature, k(r,, y) represents the spatially varying thermal conductivity, and q(x, (x) y) models internal external heat sources/sinks. The term (x, y) acts as an external heat influence.
This equation is crucial for studying heat distribution in materials with varying properties, such as in the thermal management of electronic devices. Solving this equation helps predict temperature variations, essential for designing effective cooling or heating systems.
In our experiment, k(x,y) y) is derived using the GRF 35 And q(z, y) is samples from UO, 1). The term (x, y) is consistently set constant
G Analysis of Hyperparameters
G.1 Hyperparameters of Fourier Neural Operator
Fourier Neural Operator has three key hyperparameters: model layers, mode and width for fourier layer. We conduct experiments to investigate the impacts of these hyperparameters. Results in Table suggest that proper combination of these hyper- sypar-parameters will improve the performance.
19

---

Table 8: Performance of NeurKltt with changes in layers, modes, and widths. We use Darcy Flow problem with matrix size of 32400 in these experiments, where the tolerance is fixed le-5 and the preconditioning is None. Subspace dimension is fixed to 20.
Layer
5
Mode
16
32
64
Principal Angle
0.18
0.12
0.11
0.13
Principal Angle
0.47
0.22
0.13
0.10
Time Speedup
1.47
1.68
1.69
1.54
Time Speedup
1.19
1.46
1.69
1.71
Iteration Speedup
8.15
9.86
9.45
8.04
Iteration Speedup
5.91
7.96
9.86
9.66
Width
8
16
32
64
Principal Angle
0.25
0.14
0.12
0.09
Time Speedup
0.98
1.64
1.69
1.83
Iteration Speedup
4.93
9.23
9.86
10.29
2.75
3x
235
poeticel 16.
Naneot 20-12
Romclol
Time Speedup
Jacsbiool
Iteration Speedup
25
Subspace Dimension
Subspace Dimension
Figure 4: Experiments on the Darcy Flow problem with varying subspace dimensions. The matrix size is fixed to 32400.
G.2 Impact of Invariant Subspace Dimension
We conduct extensive experiments investigate the acceleration impact of subspace dimension. Results in Figure shows that with the subspace dimension increasing, iteration speedup increase as well.
However, the time speedup does consistently increase with the subspace dimension. As shown in the left figure Figure|4) 4 the time speedup reaches its maximum when the subspace dimension 8. When the subspace dimension increase from 8, the time speedup decreases. This may be influenced by the subspace size since the subspaces are dense matrices. Recalling the Arnoldi relation in Section 4.2 larger subspaces mean larger dense matrix multiplications. When the subspace dimension is sufficiently large, the additional time cost of dense matrix operations will offset the time savings from reduced iterations, thereby reducing the overall time speedup. This indicates that larger subspace dimension is not necessarily better.
20

---

H
Detailed experimental results
H.1 Time Efficiency Analysis for Subspace Prediction Module
Table 9: The inference time spent of subspace prediction module
Helmholtz 62500
Heat 90000
Darcy 160000
Dataset
Inference Time (ms)
6.16
7.04
8.81
Table shows the time cost of the subspace prediction module for a single pass. The first row lists datasets with matrix side lengths. It's obvious that across various datasets, the computational cost of the subspace prediction module is extremely low, typically taking only a few milliseconds.
For comparison, in Darcy Flow problem, solving a linear system using GMRES without precondi- tioning at tolerance of 1e-12 requires 34.8 seconds. While the neural network's inference cost is much more lower, which is negligible.
H.2 Training Time Analysis
Table 10: The 120 epochs training cost for subspace prediction module.
Helmholtz 62500
Heat 90000
Darcy
Dataset
160000
Training Time (h)
0.51
6.48
10.47
53503
10
Derry 100000
(Tmmmmii)
Figure Training time costs for main experiments.
We train each dataset for 120 epochs, and list the training costs in Table 10 The total training time for the Helmholtz dataset is shorter compared to the other datasets. This is because the problem simple, allowing us to train on a dataset with fewer samples, and the matrix size involved in the problem is smaller than Darcy Flow and Heat problem.
As shown the Figure 5 the training typically converges within 100 min. For the Helmholtz problem, it converges in 10 min; for the Darcy Flow dataset and Heat dataset. it converges less than 100 min. However, in scientific computing software, where solving linear systems is frequently required and often time-consuming. the time saved by our approach will add-up to pay off the training costs.
21

---

H.3 Darey Flow problem
Table 11: Comparison of time and iteration counts for two algorithms, with varying preconditioners and tolerances
Propondtionngg
Mathodk
1-F-2
1.F-06
1.F-1
1.F-12
timetg)
GMRES
5.58
=0.45
8,31
+0.71
20.94
11.11
=1.42
16.09
+1.65
19.41
22.45
s1.91
27.93
52337
#5.14
34,85
NewrKhr
1.01
10:25
40.34
60.47
10.71
10.79
7.61
11.17
GMRES
5291.44
6451.50
5068 17
4636.85
#922.37
10645.95
41293.26
16400 62
41628.79
19181.24 14
21863.57
21957 72
22338.25
27531.73
43126.39
34681.57
NeurKIt
184.59
328.15
#121.09
573.21
878.1 11
#172.60
1576.28
2250.07
1278,21
1932.16
2301.87
#315.52
2063
a377.01
#477.31
3938.56
3,30
15.79
time
#097
#2.97
Incohii
NturKItt
0.50
60.17
#0.21
1.36
1,90
+0.30
3.11
3.68
20.32
4.39
5.64
20.37
7.11
#0.37
Nier
GMRES
3058.70
#448.68
2923.94
#1410.40
5213.33
#1900.35
#2011.34
9503.62
#2200 99
11837.47
22555.77
#2808.0
NeurKhr
278.29
+54.4
#93.54
+125.27
1093.00
0.73
1.00
1.24
1.64
1.83
2.04
2.43
2,86
timeçu)
GMRES
#0.05
+0.12
s0.19
+0.31
+0.36
+0,40
=0.46
#1.52
NenrKhr
10.03
2003
40.04
0.91
20.05
20.05
1.19
10.05
20.05
#0.05
GMRES
637.75
470,08
685.41
882.00
2131.37
1231.33
4277.79
1396 76
1564.51
1910.45
#414.83
2376.63
Neurktt
68.39
+4.41
110.98
+5.43
#11.13
159.66
+13.71
262.87
312.86
+14.51
363.86
#15.75
#17.33
468.04
#18.32
574,13
timets)
2000
20.37
#1173
SOR
NourKItt
0,43
+0.02
0.57
+003
0.76
00.06
40.07
1.09
1.27
+0.02
1.43
+0.07
L.76
-0.07
2.13
+0.88
OMRES
653,60
233.76
462.79
2108.1
1247.77
2228.54
4257.30
1946.7
2145.13
1345.91
2524.11
#467.13
4579.27
NeurKIit
84.07
135.66
195.51
327.29
391.85
455.74
586.11
720.40
GMRES
0.72
20.04
20.06
1.00
20.09
1.24
20.18
1.65
20.23
1.87
2.06
40.28
2,40
a0.37
2.76
#0.47
timels)
10.02
4003
20.04
20.04
10.04
10.05
20.05
40.04
GMRES
385.30
$19.85
573.78
748.66
196404
1208 00
4150.75
2184.25
1346.39
1607.4
4318.06
1873.52
Neurktt
58.03
+3.57
+501
93.66
134.21
08.53
=11.07
220.25
+11.53
262.54
304.28
+12.72
388.87
#13.53
475,30
#13.
GMRES
8.56
27.53
timele)
40.08
#7.14
ICC
NourKItt
3.24
+0.19
5.01
+0.26
7.04
+0.45
+0.56
11.14
13.16
+0.63
+0.65
15.22
40.72
19.65
+0.76
23.81
GMRES
354.05
216.03
423. 46
245.45
2101.0
1146.1 12
4133.70
1166.3 37
1278.83
#230.08
4296.41
1776.70
NeurKKit
55.75
+4.67
$9.39
128.23
269
248.35
287 80.88
=12.49
367.34
213.44
448.51
GMRES
8:51
12.42
16.35
23:44
26.66
29
35.60
41.29
Limets)
+0.38
20.66
21.05
#2.37
43.11
23,87
21.37
25.94
NeurKltt
=0.6
20.52
#1.73
GMRES
354.15
$16.03
532.55
+23.46
701 68
145.45
1101.00
1896.67
8133.70
1146.12
6166
1278,82
#230.11
1527.27
$208.46
1776.33
NeurKlit
128
209.0
287
367.32
448.45
+3.22
+4.67
+8.27
+10.40
+11.86
+12.46
+13.32
+14.
22

---

H.4 Non-uniform Heat Conduction Equation
Table 12: Comparison of time and iteration counts for two algorithms, with varying preconditioners and tolerances
Preconditioning
timets)
GMRES
40.15
1.83
40.22
2.71
40.29
3.62
30.44
5.43
40.1
6.31
50.58
7.21
40.73
9,10
$0.90
11.01
NeucKltt
0.58
0.83
1.12
1.79
2:40
2.74
3.10
3.80
+0.14
10.14
10.17
+0.23
60.28
60.31
+0.33
+0.39
2687.92
4406.
9925.35
1048541
12005.11
1536642
18445.01
=355.17
#476.97
+720.00
#963.15
#1510.96
NeueKTi
470,76
304.70
883.34
401.21
8104.10
697.08
1149.72
1150 39
4160.00
1374.07
a178.52
1594.95
#217.49
2051.06
#263.78
2539,87
GMRES
1.20
1.02
2,00
2.57
2.87
3.22
4.05
5,34
time(s)
10.13
40.27
10.44
20.78
40.93
20.96
40.95
11.05
Iscabii
#0.12
20.14
#1.14
GMRES
1797.71
2494.54
4423.53
3080.30 30
4100.70
4140060
4580.22
41553.22
5196.70
6579.74
#1727.0
8193.87
NeueKTt
285.977
429.03
596.15
897.17
105392
1215.19
1540.5
1882.29
$42.54
151.99
161.37
274.56
277.06
280.02
=75.00
182.20
1.12
tinets)
10.02
40.03
10.03
20.05
40.06
80.07
10.10
10.12
Blacobi
Neu
0.21
40.02
204.42
452.82
608 86
905.63
104487
1180.64
1494.14
1679.72
Rier
GMRES
#17.00
+28.30
#37.22
+58.9
#73.22
288.20
=134.40
=162.88
127.04
243,03
351,18
424,54
25.27
47.47
210.84
213.44
414.38
210.12
#10.99
#16.41
0.39
0.56
1.01
1,29
time(s)
20.11
60.16
49.22
SURR
NeueKTl
10.04
0.25
10.04
0.31
10.04
0.39
80.04
0.57
40.05
0.67
60.04
0.72
10.04
0.91
10.04
1.06
Ber
GMRES
+21.77
406.26
+38.00
+61.21
789.00
+131 64
1113.74
+170.0
200.51
1405 04
-86.34
1677.24
+352.50
1954.73
NeueKTu
205.04
430.65
120.55
#20.84
GMRES
11330
0.41
0.49
0.66
0.75
0.84
0.97
1.13
sime(s)
40.02
20.04
Neu-kttt
0.21
0.29
0.33
0.45
0.1
0.58
0.70
0,80
10.02
40.02
10.02
+0.00
40.03
20.03
40.03
10.00
235.58
755.72
#14.47
#17.65
#57.18
289.71
#111.94
ierr
NeucKltt
47.14
73.28
103.84
167 23
198.55
238.53
287.15
345.50
43.90
45.16
47.59
411.42
=12.52
212.87
limefe)
GMRES
10.14
2.84
10.22
4,31
10.21
5.50
10.46
8:00
30.66
9.23
80.83
10.50
21.06
12.72
11.35
14.87
1.49
5.47
7.72
ICC
+0.15
+0.20
+033
+0.34
+0.34
GMBES
+10.71
214.30
334.07
#16.08
423.08
611.00
700666
787.07
958.47
1125.41
NeueKltt
#3.69
44.41
#4.52
69.62
#7.22
98.37
158.96
29.33
#11.13
187.98
211.06
215:77
=11.75
270.79
#12.05
325,39
time(s)
GMRES
20.14
2.85
40.22
40.21
5.58
40.46
40.06
20.82
10.35
11.04
$1.33
14,64
ILU
NeueKCu
1.45
2996
4.62
5.39
6.12
7.64
GMRES
214.30
324.07
423.08
611.08
700.66
787.07
958.47
1125.41
+10.71
+16.36 66.36
+16.88
+35.00
+90.37
#62.61
+79.77
+102.30
NeueKln
69
187.98
215.77
270.79
+3.00
+4.52
+7.23
+9.38
+11.13
=11.06
=11.75
+12.06
23

---

H.S
Helmholtz Equation
Table 13: Comparison of time and iteration counts for two algorithms, with varying preconditioners and tolerances
Proconditionnng
1.F-03
0.68
1.01
133
1.94
2.23
2.55
3.18
3.79
timels)
2001
10.01
2001
10.02
20.02
10.02
10.02
None
NourKitt
0.28
0.39
0.49
+0.03
0.57
+0.04
0.98
1.26
1,50
GMRES
1578.26
2101.40
2825
4274,92
5000.12
5725.39
7176.56
8710.98
#16.11
NourKitt
179.99
262.28
353.46
565.67
686.25
811.86
1064.24
1327.96
+12.08
+17.57
#35.47
GMRES
0.20
1.01
132
1.96
2.26
3.20
3.86
timels)
40.02
20.02
20.02
Jacobi
NourKhtt
0.30
0.10
0.19
#0.03
0,72
0.85
1.00
1.27
1.56
GMRES
1378.28
2101.40
2825.34
4274,92
5000.12
5725.39
7176.56
155.41
8711.24
NourKitt
179.99
112.08
262.28
353.45
$65.66
686.26
811.87
1061.27
=34.27
1327.96
135.54
GMRES
0.1
0.27
0.34
0.46
0.56
0.67
0,76
10.00
10.00
time(s)
Blucobi
NeurKitt
0.16
0.18
0.25
40.01
0.30
20.01
0.33
0.37
0.43
0.50
GMRES
226.80
345.11
464.77
691.76
794.27
897.95
1109
1301.74
NeurKitt
41.30
71.43
108.44
168.11
194.55
231.47
276 24
332.82
+5.96
GMRES
0.2
0.31
0.38
0.56
0.63
0.70
0.85
0.98
10.00
time(s)
SOR
NeurKitt
0.15
0.22
0.27
0.34
20.01
0.39
0.42
0.48
0.57
GMRES
278.40
426.14 14
574.14
873.33
1024.37
1173.98
1472.20
1776.41
NeurKitt
49.15
92.04
131.78
193.15
225.55
257.19
$20.99
389.45
+5.08
#2.67
+3.85
=7.00
#6.26
GMRES
0.21
0.27
0.32
0.43
0.48
0.53
0.64
0,76
+0.00
time(s)
ASM
NeurKitt
0.16
+0.00
0.18
60.00
0.23
+0.00
0.28
=0.01
0.33
0.34
+0.1
0.4
0.47
GMRES
173.48
260.62
347.40
$21.99
609.52
696.81
$75.00
1056.19
+3.73
NeurKitt
35.89
55.89
86.75
134.46
155.76
176.42
221.61
266.27
#1.33
#1.29
$4.50
$4.71
$4.96
#6.07
GMRES
1.33
2.00
2.66
4.03
4.71
5.40
6.79
8.18
+0.03
+0.05
+0.03
ICC
NeurKitt
0.87
10.01
1.14
20.04
1.66
40.02
2.58
20.00
2,50
20.05
3.25
#0.06
3.99
4.86
GMRRS
146.86
223.38
300.63
=2.27
455.72
+3.05
533.56
#3.21
611.70
+3.53
768.22
+4.46
925.17
NttrKitt
34.04
53.20
80.74
125.56
141.98
164.42
204.95
249.09
+0.37
10.42
#1.81
40.75
11.39
42.42
42.97
43.62
GMRES
1.29
1.95
2.61
3.97
4.65
5.32
6.68
8.01
timels)
+0.02
s0.02
+0.03
=0.08
#0.03
#0.05
ILU
NturKitt
40.01
0.96
10.01
1.13
20.05
1.64
40.02
2.54
10.03
2.78
10.05
3.22
10.06
3.95
40.07
4,82
146.86
223.38
500.6
455.72
$33.56
611.70
768.22
925.20
GMRES
#2.27
43.05
$3.21
43.53
$4.46
NeurKitt
34.04
$3.20
80.74
125. 135.56
141.98
164.42
204.95
240.08
40.37
11.81
10.75
11.39
42.42
12.97
43.61
24

---

I Do we need to learn the neural network for matrices of a fixed size?
Table 14: Experiments on using subspace prediction module trained 9604 Darcy Flow dataset to accelerate the datasets with different matrix sizes. Results indicate that NeurKItt able to accelerate linear systems with different sizes of the matrix when the parametric PDE is given.
Dataset
Tol
None
Jacobi
BJacobi
SOR
ASM
ICC
ILU
3.05
8.75/8.75
1.27/4.85
1.21
1..7
2.15/344
2.14/
le-4
8.19/8.19
1.9966116
1.42/4.47
1,43/4.45
1.19/3.51
1.29/3.09
1.28/3.09
Darcy 39204
le-7
2.65
1.72/4.91
1.45/3.92
1.88
1.27
1.27
1.26/2.87
lc-10
6.79
1.0/4,27
1.27
1.37
1.25/3.12
1.26/2880
1.24
le-12
2.36/ 6.60
1.54/4.10
/3.57
1.36/
1,27 /3.06
1.27 /2.80
1.26/2880
le-2
4.83 / 13.80
2.81/8.21
6.52/6.59
2.25
1.8266.16
2.24
2
le-4
41.2/11.25
2,42/6,49
2.09 5.06
2.14/5.01 14,9
/5.01
Darcy 158404
le-7
3,50/8,92
2.06/
1.12/4188
1.87 14.7
1.66
1.10/4117
1.78/4.1
1c-10
3.15
192/4.93
/3.89
1.69/4.11
1.61 /3.75
1..22
1.64 1,,22382
Ic-12
3.07 /7.99
2.01/5.13
1.58
1.64
1.56
/.629
1.60/ 3.69
An important question is whether model trained on a fixed-size dataset can be used to predict other matrices under the same PDE problem, which common scientific computing. this section, we attempt to address this question using Darey Flow dataset. First, we downsample the input function of the larger Darey Flow dataset generate new Darcy Flow datasets according the discretization method mentioned in B Specifically, we downsample the original dataset with a matrix size of 158,404 obtain Darcy Flow datasets with matrix sizes of 39,204 and 9,604.
Next, we use the subspace prediction module trained on the Darcy Flow dataset with a matrix size of 9,604 predict the subspaces of larger datasets. The input to the module the input function from the larger matrix dataset, which is downsampled to the same size as the input function of the 9,604 matrix dataset, and the output is the predicted subspace of the original matrix dataset size, recovered using cubic spline interpolation.
In experiments, the principal angles obtained through this method are generally small. For example, the model trained on the 9,604 matrix dataset predicted principal angle of 0.07 for the test set of the 39,204 matrix dataset and 0.08 for the test set of the 158,404 matrix dataset. The experimental results in Table 14 are based on the acceleration results using the predicted invariant subspaces.
These results indicate that our method can scale across different matrix sizes for given parametric PDE, which demonstrates that NeurKItt is practical. Additionally, using this approach for subspace prediction significantly reduces the required training time.
25

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The main claims in the abstract and introduction accurately reflect the contri- bution and scope of the paper.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification:
We refer Section for the discussion of our limitations.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g.. independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
the paper. The authors should use their best
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
26

---

We have cite the specific paper with the assumptions and proofs for the Theorems used our paper.
Justification:
Guidelines:
The answer NA means that the paper does not include theoretical results. All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
Justification: We provide the pseudocode for our algorithm in Appendix A and provide the details about experimental settings and hyperparameter of our model in Appendix P The dataset generation method is provided in Appendix
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer this question will not perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution be necessary to either make it possible for others to replicate the model with the same dataset, or provide access the model. In general. releasing code and data often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
specific model and empirical evaluation, may
While NeurIPS does not require releasing code, the conference does require submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how reproduce that algorithm.
(b) If the contribution is primarily
model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e.g., a large language model), then there should either be way access this model for reproducing the results way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reprodueibility may be tricky in some cases, in which authors are welcome to describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers have some path reproducing verifying the results.
5. Open access to data and code
27

---

Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
Justification:
We refer to Section for code repo
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips.co cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (c.g., for a new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https
//nips //ips.cc//pill wttpcccpppppppppeeeeeeeeeeeeeeee hotaa//gggidee//eeeeabiiiisssss /Cc////oddiiisssseeeeeee ssiiiiiiccyyyyyyy for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary understand the results?
Answer: [Yes]
Justification:
We provide details about our experiments and settings in Appendix
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core the paper to a level of detail that is necessary appreciate the results and make sense of them.
The full details can be provided either with the code, appendix, as supplemental material.
7. Experiment Statistical Significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer:
[Yes]
Justification: We provide the stand deviation of average time cost and average iteration count in Appendix H
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall with given experimental conditions).
28

---

The method for calculating the error bars should be explained (closed form formula, call a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma enror bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are range (e.g. negative enror rates).
If error bars reported in tables plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification:
We provide the details about experiments compute recources in Appendix F22
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics http:///wwwrgp...cppaaaaa/accaaa
Answer: [Yes]
Justification: The research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special cireumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g.. if there is special consid- eration due to laws or regulations their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: Contents in our paper focus on accelerating the linear system solving, which could have potential help the industry or other field involving linear system solving.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
29

---

Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), faireess considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: Our data and model focus linear system solving, especially PDE problems, which have no high risk for misuse.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification:
We have checked this and confirmed the paper poses no such risks.
Guidelines:
The answer NA means that the paper does use existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
30

---

If assets are released, the license, copyright information, and terms of use in the thcode. has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
package should be provided. For popular datasets, paperswi
latasets.paperssit
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information not available online, the authors are encouraged to reach out to the asset's s creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification:
We provide the instruction in supplemental material about our code and data.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moodl as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsoureing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then much detail as possible should be included in the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country institution) were obtained?
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
31

---

We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
32


Reviews:
[{"review_id": "Ay5FHTO9IO", "reviewer": "Reviewer_CSA7", "summary": "The paper develops a new method, NeurKItt, for accelerating the solution of non-symmetric linear systems. NeurKItt constructs an approximate invariant subspace of non-symmetric matrix A using the Fourier Neural Operator. This invariant subspace is then used as the initial subspace within GMRES. The paper provides theoretical support for the proposed method and empirical results showing the benefits of NeurKItt.", "questions": "See weaknesses.", "limitations": "The authors have adequately addressed the current limitations of the method.", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "1.) The paper develops an interesting new method for solving non-symmetric linear systems based on operator learning. In particular, I have not seen the Fourier Neural Operator used in this way before. I believe the idea has the potential to be quite useful and could have valuable implications for other problems.\n\n2.)  The authors provide theoretical analysis supporting the method.\n\n3.) The method performs well in practice, which is what is most important. NeurKItt seems to yield significant reductioning in the number of GMRES iterations, and also yields good wall-clock time speed-ups (which is particularly impressive as these matrices are very sparse). So I think the method has the potential to be very useful in practice. \n\nOverall, the paper proposes an exciting new idea that appears to work quite well and has the potential to be quite useful.\nMoreover, the paper addresses an important hard problem: improving the solutions of non-symmetric linear systems. This setting is quite challenging, as the non-symmetry makes things difficult. Often, methods for accelerating GMRES are highly problem-dependent, so a method that can lead to generic improvements is significant. I'm happy to recommend its acceptance to NeurIPS provided an issue given below in the weaknesses section is properly addressed.", "weaknesses": "I'd say the main weakness of the paper is its presentation. In particular, a significant issue is that the paper lacks a precise description of how FNO is implemented to learn the invariant subspace. The appendix only briefly describes FNO but does not explain how the paper applies it. The current discussion in lines 148-151 is vague and unclear.  In particular, the output of FNO is the output of an operator on a function, not a subspace. So, how are you getting the predictive subspace from the output of the FNO? I have ideas of how this is done, but rather than guessing, I would like the authors to explain the precise procedure adequately. This should then be added to the paper, ideally in the main body, but at least the supplement with clear pointers in the main paper of where to find the details. \n\nI'm willing to consider raising my score if the authors address this point well. If they don't, I will lower my score, as I cannot recommend acceptance of a paper for which a significant part of the methodology is unclear.  \n\nAside from this, the authors should do a careful spell-check, as there are many typos throughout, which is somewhat distracting.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the insightful and valuable comments. We respond to your comments as follows and sincerely hope that our rebuttal could properly address your concerns. If so, we would deeply appreciate it if you could raise your score. If not, please let us know your further concerns, and we will continue actively responding to your comments and improving our work.\n\n## Weaknesses\n\n> **How are you getting the predictive subspace from output of the FNO? I have ideas of how this is done, but rather than guessing, I would like the authors to explain the precise procedure adequately. This should then be added to the paper, ideally in the main body, but at least the supplement with clear pointers in the main paper of where to find the details.**\n\n- **Here, we use the 2D Darcy flow problem to demonstrate how the Fourier Neural Operator (FNO) predicts the subspace given the input function.**\n  - Let the input function be $a\\in\\mathbb{R}^{d_a\\times d_a}$ from the 2D Darcy flow problem, where $d_a\\in\\mathbb{N}$ is the resolution of the input function, which yields a linear system $Ax=b$ for numerical solving, where matrix $A\\in\\mathbb{R}^{d_A\\times d_A}$. We aim to predict a subspace $\\mathcal{\\hat{K}}$, as a matrix $X\\in\\mathbb{C}^{d_A\\times n}$, for matrix $A$ in the linear system. First, the FNO's lifting layer $\\mathcal{R}$ transforms $a$ to a higher-dimensional representation $v_0\\in\\mathbb{C}^{d_a\\times d_a\\times c}$, while the 3rd dimension is channel dimension. After $T$ Fourier layers forward, we have the $v_T\\in\\mathbb{C}^{d_a\\times d_a\\times c}$. Because Fourier layers keep the shape unchanged, we apply a transformation $Q$ to map the $v_T$ to the desired space, $\\mathcal{\\hat{K}}=Q(v_T)$, with $Q:\\mathbb{C}^{d_a\\times d_a\\times c}\\rightarrow\\mathbb{C}^{d_A\\times n}$.\n  - In practice, transformation $Q$ is a stack of transformation layers. It first flattens the first and second dimension of $v_T$, obtaining $q_0\\in\\mathbb{C}^{d_a^2\\times c}$. Then a fully-connected neural network (FNN) applies the mapping  $\\mathbb{C}^{d_a^2\\times c}\\rightarrow\\mathbb{C}^{d_a^2\\times n}$ to $q_0$, obtaining $q_1\\in\\mathbb{C}^{d_a^2\\times n}$. And another FNN applies the mapping $\\mathbb{C}^{d_a^2\\times n}\\rightarrow\\mathbb{C}^{d_A\\times n}$ to $q_1$, obtaining the output $X\\in\\mathbb{C}^{d_A\\times n}$. Finally, we apply QR decomposition to $X$ for orthogonalizing, obtaining $\\mathcal{\\hat{K}}=\\mathrm{span}\\{X\\}$.\n- **To address the problems identified in the reviewers' comments, we will make several changes to our paper.** First, we will replace lines 148-151 with more precise content to clarify the subspace prediction task. The revisions include the problem setup and the detailed description of how to predict subspace by FNO. Second, we will add the 2D Darcy flow example mentioned above to Appendix C to give an image of how to predict the subspace by the FNO. The sentences quoted below are the revised version of lines 148-151:\n\n  - \"Generally, for a linear system $Ax=b$ derived from the parametric PDE problem, to predict its invariant subspace $\\mathcal{\\hat{K}}$, the input to FNO is the input function $a\\in\\mathbb{R}^{d_a}$ from the given PDE, where $d_a\\in\\mathbb{N}$. We provide a detailed discussion in Appendix B about how to build a linear system problem from a PDE problem, and what is the input function. Our task is to learn the mapping between two Hilbert space $\\mathcal{G}:\\mathbb{R}^{d_a}\\rightarrow \\mathbb{C}^{d\\times n}$ using FNO. \n\n  - For FNO, the lifting transformation $\\mathcal{R}$ first lifts the input $a$ to a higher dimensional representation $v_0\\in\\mathbb{C}^{d_a\\times c}$, where $c$ is the number of channels. Then we feed the $v_0$ to Fourier layers. After $T$ Fourier layers forward, we have  $v_T\\in\\mathbb{C}^{d_a\\times n}$ from the last Fourier layer, which keeps the same shape as $v_0$. The FNO's output $X=Q(v_T)$ is the projection of $v_T$ by the transformation $Q:\\mathbb{C}^{d_a\\times c}\\rightarrow\\mathbb{C}^{d_A\\times n}$. NeurKItt then uses QR decomposition to orthogonalize the matrix $X$, obtaining the predicted subspace $\\mathcal{\\hat{K}}$. We provide more details about how to predict the subspace given the 2D Darcy flow problem in Appendix C.\"\n- Due to NeurIPS rebuttal restrictions, we cannot submit a revised version during the rebuttal. We will incorporate these modifications into our paper in the future version. If these modifications do not adequately address the issues raised, please let us know your further concerns, and we will continue to actively respond to your comments.\n\n> **Aside from this, the authors should do a careful spell-check, as there are many typos throughtout, which is somewhat distracting.**\n\n- Thank you for bringing this to our attention. We will conduct a thorough spell-check and correct any typos in the future version.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "LGfpnDGYE1", "reviewer": "Reviewer_XS9H", "summary": "This paper proposes a data-driven approach to accelerate solving linear systems. Linear Systems are widespread in scientific computing applications like solving partial differential equations (PDEs), nonlinear systems, etc. so this method has potential for major impact. The proposed method builds upon the idea of neural operators that learn mapping between function spaces.\n\nWhile most prior methods have attempted to accelerate solving linear systems by predicting a better initial guess, this work instead predicts the matrix invariant subspace. It uses this subspace to accelerate the Krylov Iterations. The method is validated on linear systems originating from PDEs and achieves speedups of around $5.5\\times$ in computation time.", "questions": "1. Are there any cases where GMRES (with no-preconditioners) fails to converge while NeurKItt can solve it because it is data-driven?\n2. Is it necessary to learn the NO for each system? Or is it possible to learn a single network and reuse it?\n3. (Maybe I missed this) Is it possible to include the training times for the neural networks?", "limitations": "N/A", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "1. The presented method is novel, and to the best of my knowledge there doesn't exist a method that uses a similar approach\n2. Experiments are well designed and use strong baseline methods (GMRES from PETSc). Speedup in computation time and iteration count over such a strong baseline presents a strong case for the proposed method.", "weaknesses": "1. The Neural Operators need to be trained on each individual problem. So while there is an improved convergence speed, there is an increased overall time when the training time is taken into account. However, this is a pretty common shortcoming across neural operator approaches and hence I haven't used this shortcoming as part of my overall scoring.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the insightful and valuable comments. We respond to your comments as follows and sincerely hope that our rebuttal could properly address your concerns. If so, we would deeply appreciate it if you could raise your score and your confidence. If not, please let us know your further concerns, and we will continue actively responding to your comments and improving our work.\n\n## Weaknesses\n\n> **The Neural Operators need to be trained on each individual problem. So while there is an improved convergence speed, there is an increased overall time when the training time is taken into account. However, this is a pretty common shortcoming across neural operator approaches and hence I haven't used this shortcoming as part of my overall scoring.**\n\n- We agree with the reviewer that training the neural network involves additional costs, and it is a shortcoming of any neural operator approach. But considering the benefits of acceleration over millions of linear system solving, such training costs could be negligible due to the time saved. This is one of the reasons why we claim the feasibility of NeurKItt in practice.\n\n## Questions\n\n>  **Are there any cases where GMRES (with no-preconditioners) fails to converge while NeurKItt can solve it because it is data-driven?**\n\n- There might be the case where GMRES fails but NeurKItt can solve.\n- Specifically, NeurKItt consists of the subspace prediction module and acceleration module. The subspace prediction module is data-driven, while the acceleration module uses the predicted subspace to improve the Krylov iteration convergence speed. The case mentioned is due to the acceleration module that replaces the original Arnoldi relation with $(I-C_kC_k^H)AV_{m-k}=V_{m-k+1}\\underline{H}_{m-k}$, which takes the invariant subspace into iteration. It can also be considered as a preconditioned linear system solving, which replaces the original linear system problem $Ax=b$ with $PAx=Pb$, where $P=(I-C_kC_k^H)$.  This change directly improves the convergence speed while lowering the condition number of the given matrix $A$, which leads to the case that GMRES fails to but NeurKItt can solve. The data-driven module (Subspace Prediction Module) makes the improvement possible, but it's not the direct answer.\n\n>  **Is it necessary to learn the NO for each system? Or is it possible to learn a single network and reuse it?**\n\n- We'd like to give answers given different contexts.\n- **First, given the case that systems derived from different PDEs**, we did not conduct relevant experiments. But we think it is possible to learn a single network and reuse it in the given context, while the network is pre-trained. Recent works like \\[1\\], \\[2\\], and \\[3\\], give experimental proofs that the neural operator can be a pre-train model, which can be reused for different PDEs with only fine-tuning.\n- **Second, given the case that systems derived from the same PDE but with different scales**, like different sizes of the linear systems, we have conducted an experiment to show the potential of NeurKItt in Appendix I. In particular, we train the neural operator on linear systems with a fixed matrix size but predict linear systems with larger sizes, then use the predicted subspace to accelerate the solving for those linear systems with larger sizes. The experimental results show that NeurKItt successfully accelerates the solving, which indicates the potential to learn a single network and reuse it for the linear systems derived from the same PDE but with different scales.\n\n> **(Maybe I missed this) Is it possible to include the training times for the neural networks?**\n- Yes, we have reported the training time for each dataset in Appendix H.2. In our experiments, all the training converges within 150 minutes, and keeping training has no additional performance improvement once reaching the convergence.\n\n[1] Yang, Liu, et al. \"In-context operator learning with data prompts for differential equation problems.\" *Proceedings of the National Academy of Sciences* 120.39 (2023): e2310142120.\n\n[2] Hao, Zhongkai, et al. \"DPOT: Auto-regressive denoising operator transformer for large-scale pde pre-training.\" *arXiv preprint arXiv:2403.03542* (2024).\n\n[3] Zhou, Anthony, et al. \"Strategies for Pretraining Neural Operators.\" *arXiv preprint arXiv:2406.08473* (2024).", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "U3P7wu0h06", "reviewer": "Reviewer_iSRc", "summary": "The manuscript introduces a novel method, referred to as NeurKItt, which combines neural network techniques with Krylov subspace methods to accelerate the solution of linear systems derived from partial differential equations (PDEs). The core innovation of NeurKItt lies in its capability to predict invariant subspaces associated with the matrices defining these linear systems, thereby significantly reducing the number of iterations required for convergence.", "questions": "In reviewing the theoretical analysis provided between lines 237 and 248, I observed a notable similarity in both the phrasing and the mathematical formulations with those presented in reference [23]. Could the authors clarify the extent of originality in these sections? It is crucial for academic integrity to distinguish clearly between direct quotes and original analysis.", "limitations": "The authors discussed limitations in Section 7 Limitation and Conclusions.", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "1. **Originality:**\nThe manuscript introduces a novel approach by integrating neural network techniques with Krylov subspace methods to predict invariant subspaces of linear systems, significantly enhancing the efficiency of these traditional methods.\n\n2. **Quality:**\nThe manuscript effectively targets the critical issues of computational inefficiency and instability in the context of high-dimensional and poorly conditioned linear systems. The choice of problem is highly relevant to both academic research and practical applications in scientific computing, making the work valuable to a wide audience.\n\n3. **Clarity:** \nThe manuscript is well-structured, presenting complex concepts in an understandable manner. It effectively communicates the challenges of existing methods and how NeurKItt addresses these issues, providing clear explanations and logical progression from problem statement to solution. However, improvements in typographical accuracy and some clarifications in methodological descriptions could further enhance clarity.\n\n4. **Significance:**\nThe significance of this work lies in its potential impact on the efficiency of solving large-scale linear systems, which are crucial in various scientific and engineering applications. By reducing computational time and resource consumption, NeurKItt could significantly benefit fields reliant on large-scale computations, making this contribution highly relevant to both academic research and industry applications.", "weaknesses": "1. The sentence \"Research in using neural networks for accelerating linear system solving\" on line 75 appears abruptly terminated and lacks a verb or continuation that would complete the thought.\n\n2. On line 81, the term \"precondition\" is likely a typographical error. The correct term, given the context, could be \"preconditioning\".\n\n3. There appears to be a typographical error in Equation (8), where the variable $z$ should presumably be $y$.\n\n4. Equation (13) on page 7 starting with an unnecessary 's' before the equal sign seems to be a typo.\n\n5. In the manuscript, Section 4.1 discusses the use of Fourier Neural Operators (FNO) for parametric PDE problems and introduces the subsequent subspace prediction concept. It employs a subspace $S$ defined as the span of vectors $\\{s_1, s_2, \\cdots, s_n\\}$. The manuscript mentions that $s_i$ is the eigenvector corresponding to the $i$-th smallest eigenvalue, but it does not specify whether these eigenvalues are of the matrix  or another related matrix. This omission could lead to confusion about the origin and relevance of these eigenvectors, particularly in how they are integrated into the model and influence the outcome of the subspace prediction. \n\n6. In Section 4.1, the manuscript proposes using invariant subspaces associated with the smallest eigenvalues to train the subspace prediction module. However, an oversight in the theoretical analysis (Section 5.2) arises from applying this concept to potentially non-Hermitian matrices. The section (section 5.2) is based on Hermitian positive definite matrices, whose eigenvalues are all real numbers that can be orderly compared. Non-Hermitian matrices, on the other hand, may exhibit complex eigenvalues, complicating or nullifying the notion of \"smallest\" eigenvalues due to their lack of a natural order. This discrepancy raises concerns about the applicability of the method to non-Hermitian matrices, which are prevalent in many practical applications.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the insightful and valuable comments. We respond to your comments as follows and sincerely hope that our rebuttal could properly address your concerns. If so, we would deeply appreciate it if you could raise your score and your confidence. If not, please let us know your further concerns, and we will continue actively responding to your comments and improving our work.\n\n## Weaknesses\n\n> **Typos mentioned in Weakness 1-4**\n\n- Thank you for pointing out these typos. We will fix them in the future version.\n\n\n> **The manuscript mentions that $s_i$ is the eigenvector corresponding to the $i$-th smallest eigenvalue, but it does not specify whether these eigenvalues are of the matrix or another related matrix. This omission could lead to confusion about the origin and relevance of these eigenvectors, particularly in how they are integrated into the model and influence the outcome of the subspace prediction.**\n\n- Given a linear system $Ax=b$, these eigenvalues are derived from the same matrix $A$. Thank you for pointing out this potentially misleading statement. We will fix it in the future version. In particular, we will replace the sentence in lines 164-165 \"where $s_i$ is the eigenvector corresponding to the $i$-th smallest eigenvalue, $i = 1, 2, \\dots, n$.\" with \"where $s_i$ is the eigenvector corresponding to the $i$-th smallest eigenvalue of **the given matrix $A$**, $i = 1, 2, \\ldots, n$.\"\n\n\n> **In Section 4.1, the manuscript proposes using invariant subspaces associated with the smallest eigenvalues to train the subspace prediction module. However, an oversight in the theoretical analysis (Section 5.2) arises from applying this concept to potentially non-Hermitian matrices. The section (section 5.2) is based on Hermitian positive definite matrices, whose eigenvalues are all real numbers that can be orderly compared. Non-Hermitian matrices, on the other hand, may exhibit complex eigenvalues, complicating or nullifying the notion of \"smallest\" eigenvalues due to their lack of a natural order. This discrepancy raises concerns about the applicability of the method to non-Hermitian matrices, which are prevalent in many practical applications.**\n\n- Thank you for pointing out this problem. Analyzing perturbations for non-Hermitian matrices involves pseudo-spectral methods, which makes it difficult to investigate. To simplify the analysis, Section 5.2 explores how to select an appropriate subspace in the Hermitian positive definite case. In practice, for Hermitian positive definite matrices, the smallest eigenvalues are compared following the natural order. But **for non-Hermitian matrices, the eigenvalues will be sorted by comparing their modulus**, such that: $$|\\lambda_1|\\le|\\lambda_2|\\le|\\lambda_3|\\le|\\lambda_4|\\le\\dots\\le|\\lambda_n|$$\n  where $\\lambda_i\\in\\mathbb{C}$, $i=1,2,\\dots,n$, is the eigenvalue of a given non-Hermitian matrix $A$. We will add these details to Section 5.2 to avoid the discrepancy.\n\n## Questions\n\n> **In reviewing the theoretical analysis provided between lines 237 and 248, I observed a notable similarity in both the phrasing and the mathematical formulations with those presented in reference [23]. Could the authors clarify the extent of originality in these sections? It is crucial for academic integrity to distinguish clearly between direct quotes and original analysis.**\n\n- We are sorry for the misunderstanding caused by our improper citation. We fully adhere to academic ethics, and this issue was due to a writing oversight. Sentences in lines 237-248 provide the definitions and assumptions serving as the preliminaries for Theorem 5.2, which is mentioned in Theorem 5.2. **To keep Theorem 5.2, which is cited properly, coherent with the original text, we reused the definitions and assumptions from reference [23], and this might potentially lead to disputes.**\n- Due to NeurIPS rebuttal restrictions, we are not allowed to submit a revised paper during rebuttal. **We will make the quotes clear by adding the sentence \" The following definitions and assumptions for Theorem 5.2 are from the reference [23].\" at line 238 to indicate that content in lines 237-248 comes from the reference [23].**", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "yoGXDuIQQr", "reviewer": "Reviewer_g2gN", "summary": "The authors use a neural operator approach to generate subspaces that are used for the acceleration of Krylov subspace methods for several partial differential equations setups.", "questions": "Why the smallest eigenvalues in line 163 on page 5?\nI assume that Y_k in Algorithm 1 is the predicted subspace from the neural network, but then explain the derivation used for equation (11)? How is this motivated?\nHow is the method of the authors combined with preconditioning? The details are not given but the choice of the preconditioner is crucial and how would this depend on the learned subspace the method provides?", "limitations": "The method is limited to linear problems but this is a very broad and general class. If the method was better explained and the details clear and convincing there would be no limitations for these kinds of problems.", "rating": 5, "confidence": 3, "soundness": 2, "presentation": 1, "contribution": 2, "strengths": "The authors are interested in a very relevant problem of computational science and engineering of solving linear systems of equations. The attempt of providing a convergence analysis for the proposed method.", "weaknesses": "The details of the method remain hidden to me and it seems difficult to grasp the details of the method and many details and even the major workings of the method are unclear.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the patience to read through our paper. We are pleased to re-introduce our work and respond to your comments as follows. We sincerely hope that our rebuttal could properly address your concerns. If so, we would deeply appreciate it if you could raise your score and your confidence. If not, please let us know your further concerns, and we will continue actively responding to your comments and improving our work.\n\n---\n## Major Workings of Our Method\n\n- We propose a novel method, namely NeurKItt, for accelerating linear systems solving. NeurKItt consists of two modules, one is the subspace prediction module, while the other is the acceleration module. \n- The subspace prediction module employs the neural operator for predicting the invariant subspace $\\hat{\\mathcal{K}}$ of given linear systems $Ax=b$, as a matrix $X$ that contains eigenvectors corresponding to the smallest eigenvalues from $A$.\n- The acceleration module uses the predicted subspace for accelerating the linear system solving. It is the Krylov subspace algorithm with predicted subspace involved for iteration. \n---\n## Subspace Prediction Module\n\n- We use FNO to predict the subspace of the matrix $A$ from the given linear system $Ax=b$. We give the introduction to FNO in Appendix C.\n- The reason why we select the neural operator for subspace prediction is that the mapping from the matrix $A$ of linear system $Ax=b$ to its corresponding invariant subspace $\\mathcal{\\hat{K}}$ can be considered an operator, i.e., a mapping between two Hilbert spaces, thus solving such problems involves finding the corresponding operator.\n- We employ the projection loss function for training the FNO and apply thin QR decomposition to the output of the FNO. Projection loss is designed for learning the subspace. The idea behind it is that the predicted subspace should contain the basis vectors in the target subspace. We employ the thin QR decomposition for orthogonalizing the FNO output. We provide the details about thin QR decomposition in Appendix E.\n\n> **Why the smallest eigenvalues in line 163 on page 5?**\n\n- The FNO learns the subspace $\\mathcal{S}$ that associates with the $n$ smallest eigenvalues, assuming $\\mathcal{S}=\\mathrm{span}\\{ s_1, s_2, \\ldots, s_n \\}$, where $s_i$ is the eigenvector corresponding to the $i$-th smallest eigenvalue of the given matrix $A$, $i = 1, 2, \\ldots, n$.  We think such subspace $\\mathcal{S}$ is easy to learn because Theorem 5.2 gives proof that if changes in a matrix $A$ occur in the subspace corresponding to larger eigenvalues, the subspace associated with the smallest eigenvalues is minimally affected, as long as the changes are smaller than the gap between the smallest and larger eigenvalues. Theorem 5.2 can be found in Section 5.\n---\n## Acceleration Module\n\n- The acceleration module is the Krylov iteration algorithm with predicted subspace involved, inspired by Krylov recycling algorithms. The detailed pseudo-code about the acceleration module can be found in Appendix A.2. The key of the acceleration module is the Arnoldi relation $(I−C_kC^H_k )AV_{m−k} = V_{m−k+1}\\underline{H}_{m−k}$, which takes the subspace into iteration.\n- Here we give an intuitive idea about the acceleration module. Krylov subspace iteration solves the linear system problem by iteratively approximating the invariant subspace of matrix $A$ from the linear system $Ax=b$. When the subspace information about $A$ is given, we don't need to start the approximation from scratch, thus the acceleration module can improve the convergence speed.\n\n> **I assume that Y_k in Algorithm 1 is the predicted subspace from the neural network, but then explain the derivation used for equation (11)? How is this motivated?**\n\n- $Y_k$ is the matrix of predicted subspace $\\mathcal{\\hat{K}}$ from the neural operator. In particular, we obtain $C_k$ and $U_k$ from the thin QR decomposition on $AY_k$, such that $AU_k=C_k$ and $C_k^{H}C_k=I_k$. \n- Equation $C_k^{H}C_k=I_k$ allows us to find the basis of $\\mathcal{\\hat{K}}$, which enables us to remove the subspace $\\mathcal{\\hat{K}}$ in Equation (11)  $(I−C_kC^H_k)AV_{m−k} = V_{m−k+1}\\underline{H}_{m−k}$ by $(I-C_kC_k^H)$. This approach makes the Krylov subspace iteration more efficient because we only need to search the rest of the invariant subspace for matrix $A$.\n\n- Equation $AU_k=C_k$ makes sure the matrix $\\underline{H}_{m-k}$ is an upper Hessenberg matrix, which speedups the solving of the linear system in the later iteration procedure.\n\n> **How is the method of the authors combined with preconditioning? The details are not given but the choice of the preconditioner is crucial and how would this depend on the learned subspace the method provides?**\n\n- Supposing the original Arnoldi relation $AV=VH_0$ and the preconditioner matrix $P$, the preconditioned Arnoldi relation is $PAV=VH_1$. So, for NeurKItt, after the combination, the Arnoldi relation is $P(I-C_kC_k^H)AV=VH_2$. Here, $H_i$ is the upper Hessenberg matrices.\n- In our experiments, the learned subspace for any preconditioning method is the subspace composed of eigenvectors corresponding to the smallest eigenvalues of the given matrix $A$.\n---\nWe also provide a theoretical analysis of the acceleration module in Section 5.1, which gives the convergence analysis and explains our choice to predict the subspace rather than the solution. In Section 5.2, we further justify our choice of using the subspace composed of eigenvectors corresponding to the smallest eigenvalues for prediction.\n\nIf you are confused about any of the content in our paper, like what kind of obstacle hinders the understanding of our work, please let us know your further concerns, and we will continue to actively respond to your comments.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: Keylov subspace method is a standard choice in solving linear system, it suffers from low efficiency when the iterates is less-than-ideal, making it not practical for large-scale sparse linear system. Regarding this difficulty, the submission proposes a neural operator based method, named Neural Keylov Iteration (NeurKItt), which uses a neural operator to predict the invariant subspace of the linear system and in turn to speed up solving the problem. Theoretical analysis is provided, and extensive numerical experiments demonstrate the advantages of the method. 


From the reviews, the reviewers reach a consensus that the paper contains significant contribution with two ratings as strong accept. In particular 
 - Novel idea: solving linear systems via operator learning is new to the literature, and novel contribution to the community. Theoretical guarantees are also provided. 
 - Strong performance: the method manages to significantly reduce both number of iteration and wall-clock time, implying promising applications. 

During the author-reviewer discussions, the authors managed to addressed the reviewers' concerns, and they are easy to incorporate into the final version of the paper. Regarding Reviewer g2gN, though his rating is the lowest, his review and comment are less convincing that the others, implying that he may fail to fully follow the paper. As a result, his rating is not fully taken into my evaluation (but his comments together with Reviewer iSRc's show that the authors need to work on the presentation of the paper). Given the novel idea and strong performance of the method, I recommend for a spotlight.
Decision: Accept (spotlight)