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_cc_2024_Conference
Published: 2024-05-15 16:40:24

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.

Main paper: 
 introduction: 
 Solving linear systems is the cornerstone of scientific computing, with applications in various fields including mathematics, science, and engineering [28] . Traditional solvers rely on the Krylov subspace iteration algorithm to tackle large-scale sparse linear systems [59] . It starts with a random initial vector, which corresponds to a one-dimensional linear subspace, and progressively expands this subspace to approximate the solution by iteratively minimizing the residual error.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 highdimensional problems or 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 is that prior knowledge about the linear system's invariant subspace benefits the Krylov subspace iteration by guiding the iteration process, which reduces the required number of iterations [13] .Motivated by this insight, we introduce Neural Krylov Iteration (NeurKItt), a novel method that leverages neural networks to accelerate linear system solving. NeurKItt comprises two modules: Figure 1 : The variation in tolerance for NeurKItt compared to GMRES, where each line represents an experiment under a solving method and a specific preconditioning. Notably, the NeurKItt substantially enhances the efficiency of solving the linear systems, with a reduction in the number of iterations by up to a factor of 16 and achieving a 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 a linear system to its invariant subspace can be viewed as 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 to optimize training, improving subspace prediction performance. For the acceleration module, it leverages the property that Krylov subspace iteration approximates an 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 a 5.5× speedup. Both our theoretical analyses and experimental results collectively demonstrate the efficiency of NeuralKItt.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 of the linear system to accelerate Krylov iteration. To facilitate the subspace prediction, we design a 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 nonsymmetric matrices.2 Related WorkQuestion: 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 contribution and scope of the paper. Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in 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. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper. 
traditional numerical algorithms: 
 In the field of computational mathematics, various algorithms have been devised to address the challenge of 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 a constrained, smaller subspace. Within this context, the Generalized Minimal Residual (GMRES) method [45, 14] plays a pivotal role, particularly in addressing non-symmetric matrices.To enhance the efficiency of these Krylov subspace methods, a range of preconditioning methods are employed [9] . These methods aim to improve matrix conditioning, reduce iteration counts, and boost stability. Preconditioning varieties include Jacobi [46] , Additive Schwarz Method (ASM) [64] , and 3 PreliminariesThe precision threshold of a solution directly influences both the number of iterations 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 a focus on pinpointing 5-8 levels of optimal error precision for our analyses.Question: Does the paper discuss the limitations of the work performed by the authors? Answer: [Yes] Justification: We refer to Section 7 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 in 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 on 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 a few datasets or with a 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, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a 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 be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
krylov subspace iteration: 
 In the realm of large-scale sparse linear systems, Krylov subspace methods are frequently employed as a 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 a 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 x, thereby efficiently converging towards it.Suppose we now solve the system Ax = b, (1) where A ∈ C n×n . The m-th Krylov subspace associated with the matrix A and the starting vector r ∈ C n as follows:K m (A, r) = span{r, Ar, A 2 r, . . . , A m-1 r}.(2) Generally, r = Ax 0 -b represents the initial residual, where x 0 is typically generated randomly or selected as an initial guess for the solution. The iterative process of the Krylov subspace leads to the Arnoldi relation:EQUATIONwhere V m = (v 1 , . . . , v m ) ∈ C n×m comprises unit vectors v i ∈ C n for i = 1, . . . , m, mutually orthogonal, and V m+1 = (V m , v m+1 ) extends this set. H m ∈ C (m+1)×m is an upper Hessenberg matrix. H m ∈ C m×m is the first m rows of H m , distinct due to an extra element h m+1,m at the m + 1-th row, m-th column. e m 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 a 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.NeurKItt utilizes the framework of neural operator architectures, notably the FNO [26, 31, 27] . We begin by defining a spatial dimension d ∈ N and a corresponding domain D ⊂ R d . Our focus is on approximating operators G :A(D; R da ) → U(D; R du ).Here, a ∈ A(D; R da ) and u ∈ U(D; R du ) denote functions that map the domain D to R da and R du , respectively, where d a , d u ∈ N. Both spaces A(D; R da ) and U(D; R du ) are identified as Banach spaces.In line with the definition provided by [32] , a neural operator N : A(D; R da ) → U(D; R du ) is conceptualized as a mapping expressed byEQUATIONfor a designated depth L ∈ N. In this formulation, R :A(D; R da ) → U(D; R dv ), with d v ≥ d u ,functions as a lifting operator, while Q : U(D; R dv ) → U(D; R du ) serves as a local projection operator.The efficiency of an algorithm can vary with changes in matrix size. Furthermore, different NOs often require matrices of various sizes. Therefore, our investigation included the study of 5 distinct matrix sizes.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.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof? Answer: [Yes]Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material?Answer: [Yes] Justification: We refer to Section 4 for code repo.Guidelines:• The answer NA means that paper does not include experiments requiring code. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
method: 
 NeurKItt ‡ utilizes a neural operator to predict the invariant subspace K of the linear system, and then uses it to accelerate the Krylov subspace iteration. Generally, NeurKItt is broadly divided into two components.Subspace Prediction Module: Inspired by operator learning, we employ the neural operator capable of predicting subspaces at a low cost, with the help of QR decomposition. The aim is to precisely predict the invariant subspace K of the linear system.Acceleration Module: 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.The mapping from the linear system to its corresponding invariant subspace can be considered an operator, i.e., a 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 (FNO) [31] for subspace prediction. We give a brief introduction to FNO in Appendix C. Generally, for a linear system Ax = b derived from the parametric PDE problem, to predict its invariant subspace K, the input to FNO is the input function a ∈ R da from the given PDE, where d a ∈ 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 spaces G : R da → C d×n using FNO. For FNO, the lifting transformation R first lifts the input a to a higher dimensional representation v 0 ∈ C da×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 ∈ C da×c from the last Fourier layer, which keeps the same shape as v 0 . The FNO's outputX = Q(v T ) is the projection of v T by the transformation Q : C da×c → C d A ×n .NeurKItt then uses QR decomposition to orthogonalize the matrix X, obtaining the predicted We define the one -sided distance [51] from the subspace Q to the subspace C asEQUATIONwhere Π represents the projection operator for the associated space. Its mathematical interpretation corresponds to the largest principal angle between the two subspaces Q and C [3] , and defining P Q as the spectral projector onto Q. The objective of the subspace prediction task is to approximate the invariant subspace of the matrix A with our predicted subspace K :EQUATIONwhere θ represents the parameters of NO, and δ( K, S) is the distance between the two subspaces.Based on subsequent Theoretical Analysis 5.2, we choose S as the invariant subspace associated with the smallest n eigenvalues, assuming S = span S = span{s 1 , s 2 , . . . , s n }, where s i is the eigenvector corresponding to the i-th smallest eigenvalue of the given matrix A, i = 1, 2, . . . , n. To reduce computational complexity, considering the norm equivalence theorem in Banach spaces, we optimize using the following loss function:EQUATIONThis 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:EQUATIONWe represent the projection operator for K by Q, with P Q = QQ * , leading to the relationship:EQUATIONConsequently, our projection loss function l( K, S) is defined as follows:EQUATIONWe train our subspace module by optimizing the above formula. Experiments have shown that as the projection loss decreases, the principal angle of our predicted subspace also decreases. The principal angle is a mathematical indicator derived in Theoretical Analysis 5.1 that influences the acceleration.For Krylov algorithms, the iteration count is a critical factor influencing computational load. Drawing inspiration from the Krylov recycling algorithms [40, 13, 6, 58, 41, 24] , our approach NeurKItt utilizes the predicted invariant subspace K as the deflation space within the Krylov iteration framework. We provide the pseudocode of the acceleration module in Appendix A.2. The key procedure of our implementation is outlined as follows:Considering the linear system (1), we obtain a k-dimensional invariant subspace K from the subspace prediction module. Then, NeurKItt computes matricesU k , C k ∈ C n×k from K and A such that AU k = C k and C H k C k = I k, the specific algorithm can be found in Appendix A.1. We can get the Arnoldi relation of our NeurKItt:EQUATION)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(A, r). Building upon this foundation, K(A, r) can converge more rapidly to the subspace where the solution x lies. This can 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 [8, 39, 40] .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 are 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.5 Theoretical Analysis 
convergence analysis: 
 Let Q be an l-dimensional invariant subspace of matrix A, and let C = range(C k ) be a k-dimensional space (k ≥ l) selected to approximate Q. We refer to Theorem 3.1 in [40] for an in-depth convergence analysis under NeurKItt. Theorem 5.1. [40] Considering a space C, define V = range(V m-k+1 H m-k ) be the (m -k)dimensional Krylov subspace generated by NeurKItt as in (11) . Let r 0 ∈ C n , andEQUATION)where γ = ∥(I -Π C )P Q ∥ 2 and δ = δ(Q, C).The left-hand side is the residual norm subsequent to m -k iterations of NeurKItt, employing the predictive subspace C. In contrast, on the right-hand side, the initial term epitomizes the convergence of a deflated problem, given that all components within the subspace Q have been eradicated [39, 49] . The subsequent term on the right embodies a constant multiplied by the residual following m -k iterations of NeurKItt, when solving for r 1 . Supposed that the predictive space C encompass an invariant subspace Q, then we have δ = γ = 0 for the given Q, which ensures that the convergence rate of NeurKItt matches or surpasses that of the deflated problem. In most cases, ∥P Q ∥ 2 is numerically stable and not large, therefore a reduced value of δ directly correlates with faster convergence in NeurKItt.Now we 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 ∥A∥ 2 tends to be large. Our theoretical analysis in Appendix D suggests that FNO's 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 crucially speeds up the convergence of Krylov subspace iterations, markedly reducing total iterations and improving stability.The crucial questions about predicted subspace arise when accelerating Krylov subspace iterations:(1) What kind of subspace should be 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 a 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, we 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 a Hermitian positive definite matrix A and a corresponding Hermitian perturbation E, allowing A to have the eigendecomposition [24] :EQUATIONwhereEQUATIONj1 ), Λ 2 and Λ 3 are defined in a similar fashion. Let the eigenvalues satisfy the following ordering * :λ (1) 1 ≤ . . . ≤ λ (1) j1 < λ (2) 1 ≤ . . . ≤ λ (2) j2 < λ (3) 1 ≤ . . . ≤ λ (3) j3 .We consider the change in the invariant subspace range(Q 1 ) under a symmetric perturbation E of A. Let θ 1 (•, •) denote the largest canonical angle between two spaces. We do not require that ||E|| F be small, but we assume that the projection of E onto the subspace range([Q 1 Q 2 ]) is small. We assume that ||[Q 1 Q 2 ] H E|| F ≤ ϵ and that ϵ is small relative to λ (2) 1 -λ (1)j1 . We also assume that η = ||Q H 3 E|| F is small relative to λj1 . Note that we do not need to assume that λj2 -λ(1) j1 is large. Also, let µ ≡ min(λEQUATIONEQUATIONEQUATIONj1 + ϵ.Theorem 5.2. [24] Let A be Hermitian positive definite and have the eigendecomposition given in (13) , and let E, ϵ, η, µ, and μ be defined as above. Then there exists a matrix Q1 conforming to Q 1 such that range( Q1 ) is a simple invariant subspace of A + E, andtan θ 1 (range(Q 1 ), range( Q1 )) ≤ ϵ μ .The specific proof can be found in [24] . A similar bound applies to the perturbation of the eigenvalues associated with Q 1 . In the context of Theorem 5.1 and our NeurKItt, Q 1 corresponds to Q, whereas Q 2 and Q 3 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 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 error convergence rate proved in Theorem 5.1, suggests that investing in substantial computational resources for predicting a high-precision subspace is unnecessary. As long as the predicted subspace shows a distinct correlation with the matrix invariant subspace, particularly when δ is small, significant acceleration can be achieved without resorting to high precision [13] . 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 
experiment settings: 
 To comprehensively evaluate the performance of NeurKItt, we conducted extensive experiments. Our analysis centers on two primary performance metrics viewed through three perspectives. These tests are conducted on 3 datasets. Specifically, the three Perspectives are: (i) Matrix preconditioning techniques, spanning 7 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.1.Baselines. NeurKItt focuses on 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 adaptability to various types of matrices, we examined three different linear equation challenges, each rooted in a PDE: 1. Helmholtz Equation [62] . 2. Darcy Flow Problem [31, 44, 27, 36] ; 3. Non-uniform Heat Conduction Equation [48, 25, 4, 18] . For an in-depth exposition of the dataset and its generation, kindly refer to the Appendix F.4. For the runtime environment, refer to Appendix F.2.Metrics. We adopt time speedup for "GMRES solving time / NeurKItt solving time" and iteration speedup for "GMRES iteration count / NeurKItt iteration count". Speedup over 1 denotes better NeurKItt performance. We also provide the average time spent and average iteration count in Appendix H. We use principal angle in Equation 5to show how close the predicted subspace and the target subspace are. 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.5× speedup in computation time and a 16.1× 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 increasing, 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. We conducted experiments on the Darcy Flow problem with varying matrix sizes under different preconditioning methods and tolerances. We provide the results in Figure 3 . The experimental results show that our 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. We conducted ablation experiments to further validate the performance of FNO and demonstrate the effectiveness of the designed projection 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 is indeed capable of learning a better mapping from the parameter matrix to the subspace, resulting in predicted subspaces that are closer to the target subspaces.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 is merely ineffective. Since we report the principal angle in radians, it implies that the predicted subspace K is merely orthogonal to the target subspace S. Such subspace K cannot be effectively utilized to accelerate linear system solving. Besides, it can be inferred that QR decomposition we have 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.Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes]Justification: We provide details about our experiments and settings in Appendix F Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
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 a projection loss function for this task. By utilizing the predicted invariant subspace, we achieve accelerated Krylov subspace iterations. Theoretical analysis and extensive experiments demonstrate that NeurKItt effectively reduces computational costs and iteration counts.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, confidence intervals, or statistical significance tests, at 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 run with given experimental conditions).• The method for calculating the error bars should be explained (closed form formula, call to 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. 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 circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
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 to the industry or other field involving linear system solving.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.• Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness 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 a 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 or 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 a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: Our data and model focus on linear system solving, especially PDE problems, which have no high risk for misuse.Guidelines:• The answer NA means that the paper poses no 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 is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: 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?Answer: [Yes]Justification: We have checked this and confirmed the paper poses no such risks.Guidelines:• The answer NA means that the paper does not 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 a 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.• If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • 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 is not available online, the authors are encouraged to reach out to the asset's creators. 
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/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. • At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file. 
crowdsourcing 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 as 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 in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as 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. 
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 an equivalent approval/review based on the requirements of your country or institution) were obtained?Answer: [NA] Justification:Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Depending on 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.• We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to 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. 
 
 Appendix: A.1 computes matrices U k and C kNeurKItt can be modified to solve (1) by carrying over subspace K = span{y 1 , y 2 , . . . , y k } with respect to A as follows [11, 40] :[Q, R] = thin QR decomposition of AY k , C k = Q; U k = Y k R -1 ,where the matrix Y k = [y 1 , y 2 , . . . , y k ]. The matrices U k and C k to satisfyAU k = C k ; C H k C k = I k . A.2 NeurKItt Pseudocode Algorithm 1: NeurKItt Krylov subspace iterationSelect m, the maximum size of the Krylov subspace. Define k as the dimension of the predicted invariant subspace K, where K = span{y 1 , y 2 , . . . , y k }. Given matrix A, use subspace prediction module to predict the matrix Y k = [y 1 , y 2 , . . . , y k ] and set tol as the convergence tolerance. Choose an initial guess x 0 . Calculate r 0 = b -Ax 0 , andinitialize i = 1. Let [Q, R] be the reduced QR-factorization of AY k . C k = Q U k = Y k R -1 x 1 = x 0 + U k C k H r 0 r 1 = r 0 -C k C k H r 0 while ∥r i ∥ 2 > tol do i = i + 1Perform m -k Arnoldi steps with the linear operator(I -C k C H k )A, letting v 1 = r i-1 ∥r i-1 ∥ 2 and generating V m-k+1 , H m-k and B m-k .Let D k be a diagonal scaling matrix such that U k = U k D k , where the columns of U k have unit norm.V m = [ U k V m-k ] W m+1 = [C k V m-k+1 ] G m = D k B m-k 0 H m-k Solve min ∥ W H m+1 r i-1 -G m y∥ 2 for y. x i = x i-1 + V m y r i = r i-1 -W m+1 G m y Compute the k eigenvectors z i of G H m G m z i = θ i G H m W H m+1 V m z i associated with smallest magnitude eigenvalues θ i and store in P k . Y k = V m P k Let [Q, R] be the reduced QR-factorization of G m P k . C k = W m+1 Q U k = Y k R -1x i as the numerical solution, r i as the residual, i as the iteration count.We model a steady-state PDE problem as a corresponding system of linear equations using discrete numerical methods such as FDM, FEM, and FVM [52, 21, 22, 29, 10] . Numerical methods discretize a partial differential equation problem by mapping it from an infinite-dimensional function space to a finite-dimensional space, resulting in a system of linear equations.We take the process of discretizing the two-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:EQUATIONThis equation could be approximated as:EQUATIONWe use a 2 × 2 internal grid (i.e., N x = N y = 2 and ∆x = ∆y), the unknowns u i,j can be arranged in row-major order as follows:u 1,1 , u 1,2 , u 2,1 , u 2,2 .For central differencing on a 2 × 2 grid, the vector b will contain the values of f i,j = f (x i , y j ) and the linear equation system Ax = b can be expressed as:    -4 + k 2 1,1 1 1 0 1 -4 + k 2 1,2 0 1 1 0 -4 + k 2 2,1 1 0 1 1 -4 + k 2 2,2        u 1,1 u 1,2 u 2,1 u 2,2    =    f 1,1 f 1,2 f 2,1 f 2,2    .In this example, we call the discretization of function k 2 (x, y) the input function, corresponding to a 2 × 2 grid. The function inputs vary from PDE problems.The Fourier Neural Operator (FNO) aims to map a continuous input function A to its corresponding continuous output function U within the Fourier domain. For end-to-end training on FNO, the function pairs (A, U) are discretized into instance pairs (a, u). The purpose is to learn a mapping N between (a, u), which can be expressed as:EQUATIONThe mapping N consists of several sequential steps: first, the input channel is lifted using R; next, the mapping is performed through L Fourier layers {L 1 , L 2 , . . . , L L }; finally, the output is projected back to the original channel using Q.EQUATIONQ and R are pixel-wise transformations that can be implemented using models like MLP.A typical Fourier layer consists of a pixel-wise linear transformation, characterized by weight W and bias b, along with an integral kernel operator K:EQUATION)where σ is the nonlinear activation function. The integral kernel operator K 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 a ∈ R da×da from the Darcy flow problem, where d a ∈ N is the resolution of the input function, which yields a linear system Ax = b for numerical solving, where matrix A ∈ R d A ×d A . We aim to predict a subspace K, as a matrix X ∈ C d A ×n , for matrix A in the linear system. First, the FNO's lifting layer R transforms a to a higher-dimensional representation v 0 ∈ C da×da×c , while the 3rd dimension is channel dimension. After T Fourier layers forward, we have the v T ∈ C da×da×c . Because Fourier layers keep the shape unchanged, we apply a transformation Q to map the v T to the desired space,K = Q(v T ), with Q : C da×da×c → C d A ×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 ∈ C d 2 a ×c . Then a fully-connected neural network (FNN) applies the mapping C d 2 a ×c → C d 2 a ×n to q 0 , obtaining q 1 ∈ C d 2 a ×n . And another FNN applies the mapping C d 2 a ×n → C d A ×n to q 1 , obtaining the output X ∈ C d A ×n . Finally, we apply QR decomposition to X for orthogonalizing, obtaining K = span{X}.Linear System SolversThe 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, it explores the feasibility of utilizing these neural operators to predict solutions to PDE problems. Such predictions could serve as 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 of 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:EQUATION)Within the framework of NOs, the prevalent metric for error assessment is the Mean Squared Error (MSE) of the solution function or solution vector x, designated as rtol N O :EQUATIONwhere xNO 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:EQUATIONwhere x 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. In this context, N denotes the count of discretization points in a one-dimensional Ordinary Differential Equation (ODE), and D N 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 a uniformly divided grid can be analogized to a grid defined by the zeros of Fourier polynomials.The theorem suggests that with an increase in the number of grid points, resulting in a larger matrix size, the norm of the matrix ∥A∥ 2 associated with the discretized PDE also escalates. Then We consider the following equation:EQUATIONConsequently, when ∥A∥ 2 is substantially large, there emerges a notable discrepancy between rtol N O and rtol. Generally, the predictive accuracy rtol N O of the FNO lies in the range of 10 -2 to 10 -4 . However, if the norm ∥A∥ 2 approaches the magnitude of 10 5 , the relative tolerance rtol as per computational mathematics algorithms may only range between 10 0 and 10 2 . In such instances, leveraging the FNO's predicted solution as the initial vector for solving the linear system may not substantially expedite the solution process.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 N O ∼ 6 * 10 -3 , which is in line with the benchmarks reported in the foundational FNO paper. However, because ∥A∥ 2 ∼ 1.2 * 10 6 , it results in rtol ∼ 60. This empirical finding is consistent with our theoretical predictions. In such a 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 utilizing predicted solutions as initial vectors in iterative algorithms mandates a high level of precision in the neural network model. Nonetheless, it is observed that in standard application scenarios, the impact on accelerating the process is marginal.In contrast to this method, our 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.In our research, we utilize the Thin QR decomposition algorithm, a prevalent technique in numerical linear algebra [14] . The essence of Thin QR decomposition lies in its capacity to factorize a matrix into two distinct components: an 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 C n×m with the condition n ≥ m, the Thin QR decomposition is mathematically articulated as:EQUATION)where Q is an n × m orthogonal matrix that complies with Q * Q = I, and R is an m × m upper triangular matrix. In our study, the integration of QR decomposition within the neural network architecture is strategically employed to effectively resolve the orthogonal basis of a subspace.F.1 Specific parameters of the main experimentModel 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 1 × 10 -3 while the batch size is selected from {16, 32}. The number of subspace dimensions is fixed at 10 in all experiments. For darcy flow and 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 were 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.To ensure consistency in our evaluations, all comparative experiments were conducted under uniform computing environments. Specifically, the environments used are detailed as follows: CPU: Intel ® Xeon ® Gold 6246R CPU @ 3.40GHz GPU: NVIDIA GeForce RTX 3090In this study, we implemented a 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.2. 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.4. Successive Over-relaxation (SOR) [61] : SOR, a modification of the Gauss-Seidel method, incorporates a relaxation factor to enhance the process of convergence. This method reformulates the original problem by applying a weight, often improving the convergence rate for certain iterative algorithms.5. Additive Schwarz Method (ASM) [54] : ASM operates on the principle of domain decomposition, partitioning the main problem into several smaller subdomains. Each subdomain'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 decomposition, but drops certain off-diagonal elements during the decomposition, making it "incomplete". It's utilized for symmetric positive definite problems.7. Incomplete LU (ILU) [46] : ILU is based on LU decomposition, but like ICC, drops certain off-diagonal elements during the decomposition. ILU can be applied to nonsymmetric problems.We consider a two-dimensional Helmholtz equation, which can be described by the following equation [62] :EQUATIONwhere u(x, y) is the wave function, and k(x, y) denotes the spatially varying wavenumber, with f representing the source term, corresponding to origins of electromagnetic or acoustic waves. This equation is pivotal in several physical realms, notably in Acoustics, 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(x, 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 (x, y) is the constant.We consider two-dimensional Darcy flows, which can be described by the following equation [31, 44, 27, 36] :EQUATIONwhere K is the permeability field, h is the pressure, and f is a source term which can be either a constant or a space-dependent function.In our experiment, K(x, y) is derived using the GRF method [35] . The term f (x, y) is consistently set as a constant f .We consider a two-dimensional heat conduction equation in a non-uniform medium, which is formulated as follows [4, 18] :EQUATIONwhere T (x, y) is the temperature, k(x, y) represents the spatially varying thermal conductivity, and q(x, y)T (x, y) models internal or external heat sources/sinks. The term f (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) is derived using the GRF [35] . And q(x, y) is samples from U (0, 1). The term f (x, y) is consistently set as a constant f .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 8 suggest that a proper combination of these hyper-parameters will improve the performance. We conduct extensive experiments investigate the acceleration impact of subspace dimension. Results in Figure 4 shows that with the subspace dimension increasing, iteration speedup increase as well.However, the time speedup does not consistently increase with the subspace dimension. As shown in the left figure of Figure 4 , the time speedup reaches its maximum when the subspace dimension is 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 a larger subspace dimension is not necessarily better. Table 9 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 preconditioning at a tolerance of 1e-12 requires 34.8 seconds. While the neural network's inference cost is much more lower, which is negligible. 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 two datasets. This is because the problem is 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 in 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. I Do we need to learn the neural network for matrices of a fixed size? An important question is whether a model trained on a fixed-size dataset can be used to predict other matrices under the same PDE problem, which is common in scientific computing. In this section, we attempt to address this question using Darcy Flow dataset. First, we downsample the input function of the larger Darcy Flow dataset to generate new Darcy Flow datasets according to the discretization method mentioned in B. Specifically, we downsample the original dataset with a matrix size of 158,404 to 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 to predict the subspaces of larger datasets. The input to the module is 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 our experiments, the principal angles obtained through this method are generally small. For example, the model trained on the 9,604 matrix dataset predicted a 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 the given parametric PDE, which demonstrates that NeurKItt is practical. Additionally, using this approach for subspace prediction significantly reduces the required training time.‡ Our code is available at https://github.com/smart-JLuo/NeurKItt 


Decision: 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.

Reviews: review: paper_summary: paper_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.
strength_weakness: strength_weakness: 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. 

I'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.  

Aside from this, the authors should do a careful spell-check, as there are many typos throughout, which is somewhat distracting.
questions: questions: See weaknesses.
limitations: limitations: The authors have adequately addressed the current limitations of the method.
paper_summary: paper_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.

While 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.
strength_weakness: strength_weakness: 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.
questions: questions: 1. Are there any cases where GMRES (with no-preconditioners) fails to converge while NeurKItt can solve it because it is data-driven?
2. Is it necessary to learn the NO for each system? Or is it possible to learn a single network and reuse it?
3. (Maybe I missed this) Is it possible to include the training times for the neural networks?
limitations: limitations: N/A
paper_summary: paper_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.
strength_weakness: strength_weakness: 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.

2. On line 81, the term "precondition" is likely a typographical error. The correct term, given the context, could be "preconditioning".

3. There appears to be a typographical error in Equation (8), where the variable $z$ should presumably be $y$.

4. Equation (13) on page 7 starting with an unnecessary 's' before the equal sign seems to be a typo.

5. 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. 

6. 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.
questions: 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: limitations: The authors discussed limitations in Section 7 Limitation and Conclusions.
paper_summary: paper_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.
strength_weakness: strength_weakness: 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.
questions: questions: Why the smallest eigenvalues in line 163 on page 5?
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?
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?
limitations: 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.

Comments: title: Thanks for Your Positive Feedback
Thanks for Your Positive Feedback
Thanks for Your Positive Feedback
Thanks for your positive feedback
Thanks for the rebuttal
comment: We appreciate the positive feedback and the time you taken to evaluate our work. To improve clarity, we will include these discussions to help readers grasp the key aspects of our study.
We appreciate your positive feedback and the efforts you made in reviewing our rebuttal. Thanks for your constructive comments and valuable suggestions. We will incorporate them into our paper in the future version.
We sincerely appreciate your positive feedback and the time you’ve dedicated to reviewing our rebuttal. Thank you for appreciating the merits of our work. Your constructive comments are invaluable to us.
Thanks for the response of the authors.

All of my previous concerns have been addressed in the rebuttal, I think this is a good work and I would like to raise my score from 5 to 6 for this submission.
I would like to thank the authors for their detailed responses to my questions. Thank you for pointing out the relevant sections in the appendix which clearly answer my remaining questions.

Considering this and the response to Reviewer CSA7 that improves upon the explanation of the method, I have increased the score to 8.
Dear Reviewer iSRc,

We are writing as the authors of the paper titled "Neural Krylov Iteration for Accelerating Linear System Solving" (ID: 18161). Thanks again for your valuable comments and constructive suggestions, which are of great help to improve the quality of our work. We have followed your suggestions to significantly enhance the quality of our work. As the deadline for the author-reviewer discussion period is approaching (due on Aug 13), we are looking forward to your further comments and/or questions.

We sincerely hope that our rebuttal has properly addressed your concerns. If so, we would deeply appreciate it if you could raise your scores. If not, please let us know your further concerns, and we will continue actively responding to your comments and improving our work.

Best regards,

Authors
Dear Reviewer XS9H,

We are writing as the authors of the paper titled "Neural Krylov Iteration for Accelerating Linear System Solving" (ID: 18161). Thanks again for your valuable comments and constructive suggestions, which are of great help to improve the quality of our work. We have followed your suggestions to significantly enhance the quality of our work. As the deadline for the author-reviewer discussion period is approaching (due on Aug 13), we are looking forward to your further comments and/or questions.

We sincerely hope that our rebuttal has properly addressed your concerns. If so, we would deeply appreciate it if you could raise your scores. If not, please let us know your further concerns, and we will continue actively responding to your comments and improving our work.

Best regards,

Authors
Thank you for your question. We address it as follows:

- We agree with the reviewer about how preconditioners work. However, our method is based on the predicted subspace to improve the convergence speed, which avoids approximating the subspace from scratch for Krylov subspace iteration. Thus, preconditioners that disrupt the invariant subspace, not eigenvalues, will have a negative influence on our method, like ICC and ILU.

- Despite this potential influence, extensive experiments show that these disruptions do not significantly affect our improvement in practice. In particular, our method accelerates the solving across all preconditioners, especially achieving up to 2.64$\times$ time speedup under ICC and 2.6$\times$ time speedup under ILU.
- Besides, optimizations for the approximated subspace could be our future work, which takes the negative influence of different preconditioners on subspace into consideration.
Thank you for your positive feedback and for taking the time to read and respond to our rebuttal. We appreciate your constructive comments and concrete suggestions. We will incorporate them into our paper in the future version.
I appreciate the authors' detailed reply to my concerns. This adequately addresses the issues I had with the presentation. I trust the authors to include this in the final revision, as they have promised here. Given this I will raise my score from 7 to 8, and presentation from 1 to 3.
Thank you for the comments. Regarding your last point on the use of the preconditioner, the main idea of applying a preconditioner P is to cluster the eigenvalues of the preconditioned matrix PA. One typically hopes that the n previously smallest eigenvalues would now be clustered in a tighter region, ideally the perfect preconditioner guarantees only very few distinct eigenvalues What does this mean for your method and the approximated subspace?