Accelerating Data Generation for Neural Operators via Krylov Subspace Recycling
Authors: Hong Wang
Zhongkai Hao
Jie Wang
Zijie Geng
Zhen Wang
Bin Li
Feng Wu
Venue: ICLR 2024 spotlight (2024)
Keywords: AI4PDE; Neural Operator; Data Generation; Krylov Subspace
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)
TL;DR: This paper proposes a novel method, namely SKR, to boost the data generation efficiency for neural operator training based on Krylov subspace recycling.

Abstract: Learning neural operators for solving partial differential equations (PDEs) has attracted great attention due to its high inference efficiency.
However, training such operators requires generating a substantial amount of labeled data, i.e., PDE problems together with their solutions.
The data generation process is exceptionally time-consuming, as it involves solving numerous systems of linear equations to obtain numerical solutions to the PDEs.
Many existing methods solve these systems independently without considering their inherent similarities, resulting in extremely redundant computations.
To tackle this problem, we propose a novel method, namely **S**orting **K**rylov **R**ecycling (**SKR**), to boost the efficiency of solving these systems, thus significantly accelerating data generation for neural operators training.
To the best of our knowledge, SKR is the first attempt to address the time-consuming nature of data generation for learning neural operators.
The working horse of SKR is Krylov subspace recycling, a powerful technique for solving a series of interrelated systems by leveraging their inherent similarities.
Specifically, SKR employs a sorting algorithm to arrange these systems in a sequence, where adjacent systems exhibit high similarities.
Then it equips a solver with Krylov subspace recycling to solve the systems sequentially instead of independently, thus effectively enhancing the solving efficiency.
Both theoretical analysis and extensive experiments demonstrate that SKR can significantly accelerate neural operator data generation, achieving a remarkable speedup of up to 13.9 times.

Published as a conference paper at ICLR 2024
ACCELERATING DATA GENERATION FOR NEURAL OPERATORS VIA KRYLOV SUBSPACE RECYCLING
Hong Wang! Zhongkai Hao*r, Jie Wang!2t, Zijie Geng Zhen Wang! Bin Li Feng Wu
CAS Key Laboratory of Technology in GIPAS, University of Science and Technology of China MoE Key Laboratory of Brain-inspired Intelligent Perception and Cognition, University of Sci- ence and Technology of China
7000 ustoggs, wangzhen051 8 @mmail wsstc.edu.... [jiewangx, binli, fengwu Bust edu.cn
3 Tsinghua University
jzj21 imails. tsinghua..... edu.cn
ABSTRACT
Learning neural operators for solving partial differential equations (PDEs) has attracted great attention due to its high inference efficiency. However, training such operators requires generating a substantial amount of labeled data, i.c. PDE problems together with their solutions. The data generation process is exception- ally tone-conuumnng,
involves solving numerous systems of linear equations to obtain numerical solutions to the PDEs. Many existing methods solve these systems independently without considering their inherent similarities, resulting extremely redundant computations. To tackle this problem, we propose a novel method, namely Sorting Krylov Recycling (SKR), to boost the efficiency of solv- ing these systems, thus significantly accelerating data generation for neural op- erators training. To the best of our knowledge, SKR is the first attempt to ad- dress the time-consuming nature of data generation for learning neural operators. The working horse of SKR is Krylov subspace recycling, a powerful technique for solving a series of interrelated systems by leveraging their inherent similari- ties. Specifically, SKR employs a sorting algorithm to arrange these systems in a sequence, where adjacent systems exhibit high similarities. Then it equips a solver with Krylov subspace recycling to solve the systems sequentially instead of independently, thus effectively enhancing the solving efficiency. Both theoret- ical analysis and extensive experiments demonstrate that SKR can significantly accelerate neural operator data generation, achieving a remarkable speedup of up to 13.9 times.
INTRODUCTION
Solving Partial Differential Equations (PDEs) plays fundamental and crucial role in various scien- tific domains, including physics, chemistry, and biology (Zachmanoglou & Thoe 1986). However, traditional PDE solvers like the Finite Element Method (FEM) (Thomas) (0013) often involve sub- stantial computational costs. Recently, data-driven approaches like neural operators (NOs) (Lu al, 2019) have emerged as promising alternatives for rapidly solving PDEs (Zhang et al.) 2023). NOs can be trained on pre-generated datasets as surrogate models. During practical applications, they only require a straightforward forward pass to predict solutions of PDEs which takes only several milliseconds. Such efficient methodologies for solving PDEs hold great potential for real-time pre- dictions and addressing both forward and inverse problems in climate (Pathak et al. 2022). fluid dynamics (Wen al. 2022), and electromagnetism (Augenstein et al. 2023).
Despite their high inference efficiency, a primary limitation of NOs is the need for a significant volume of labeled data during training. For instance, when training Fourier Neural Operator
"Equal contribution.
1 Corresponning author.

---

Published as a conference paper at ICLR 2024
le-1
1e-2
tolerance (2 norm)
9 de
1e.3
SKR[purs] None
solver
SKRtours] jacobi
1e.4
SKR(purs) Blacobi
PDE generator
PD
solver
DATA SET
1e.5
SkRtours] ASM
SKRlours] SOR
4
GMRES
GMRES Jacabi
(o
solver
1e.6
GMRES Bjacobi
EMRES SOR
1e-7
GMRES ASM
3 time
Left. Generation process of the NO dataset. 1. Generate a set of random parameters from NO. 2. Export the corresponding PDE based on these parameters. 3. Transform the PDE into a system of linear equations using discretization methods. 4. Invoke linear equation solvers to solve these systems independently. Obtain solutions to the linear systems and convert them into PDE solutions. 6. Assemble into dataset. Right. The accuracy change over time for our SKR algorithm and the basclinn GMRES algorithm. As seen, the SKR algorithm can significantly accelcrate the solution of the system of linear equations, achieving speed-up of up to 13.9 times.
Figure 1:
(FNO) (Li et al. 2020) for Darcy flow problems, thousands of PDEs and their corresponding solutions under various initial conditions are often necessary. Acquiring this data typically en- tails running numerous traditional simulations independently, leading to elevated computational ex- penses. For solvers such as Finite Element Method (FEM) Veque 2002), accuracy often in- creases quadratically or even cubically with the mesh grid count. Depending on the desired accu- racy, dataset generation can range from several hours to thousands of hours. Furthermore, the field of PDEs is unique compared to other domains, given the significant disparities between equations and the general lack of data generalization. Specifically, for each type of PDE, dedicated dataset needs recreation when training NOs. The computational cost for generating data for large-scale problems has become bottleneck hindering the practical usage of NOs (Zhang al. 2023). While certain methodologies propose integrating physical priors or incorporating physics-informed informed loss to potentially increase data efficiency, these approaches are still in their primary stages and are not satisfactory for practical problems (Hao et al. 2022). Consequently, designing efficient strategies to generate the data for NOs remains a foundational and critical challenge (Geng et al., 2023).
The creation of the NO dataset is depicted in Figure [] Notably, the computation in the 4-th step Potentially constitutes approximately 95% of the entire process (Hughes) (012), underscoring as prime candidate for acceleration. While conventional approaches tackle each linear system independently during dataset formulation, we posit that there is an inherent interconnectivity among them. Specifically, shown in Figure and systems originating from similar category of PDEs often exhibit comparable matrix structures, eigenvector subspaces, and solution vectors (Parks et al. 2006). Addressing them in isolation leads considerable computational redundancy.
To address the prevalent issue of computational redundancy in conventional linear system solutions, introduce a pioneering and efficient algorithm, termed Sorting Krylov Recycling (SKR). SKR is ingeniously designed to optimize the process from the ground up. Its initial phase involves a Sorting algorithm, which serializes the linear systems. This serialization strategy carefully crafted to augment the correlations between successive systems, setting the stage for enhanced computa- tional synergy. In its subsequent phase, SKR delves into the Krylov subspace established during the linear system resolutions. Through the application of a 'Recycling" technique, it capitalizes on eigenvectors and invariant subspaces identified from antecedent solutions, thereby accelerating the convergence rate and substantially reducing the number of iterations and associated computation time. Central to its design, SKR discerns and leverages the inherent interrelations within these lin- systems. Rather than approaching each system as a discrete entity, SKR adeptly orchestrates their sequential resolution, eliminating considerable redundant calculations that arise from similar structures. This refined methodology not only alleviates the computational demands of linear sys- tem solutions but also markedly hastens the creation of training data for NOs. Codes are available at https://giihub..o///aggg/ge http:///..o..ooooooooooo/oooaaao 770/No-aaaaaainnnSnSS.
2

---

Published as a conference paper at ICLR 2024
2
RELATED WORK
2.1 DATA-EFFICIEN NEURAL OPERATORS AND LEARNED PDE SOLVERS
Neural operators, such as the Fourier Neural Operator (FNO) Li et al. 2020) and Deep Operator Network (DeepONet) (Lu et al. 2019), effective models for solving PDEs. However, their train- ing requires large offline paired parametrized PDE datasets. To improve data efficiency. research has integrated physics-informed loss mechanisms similar to Physics Informed Neural Networks (PINNs) (Raissi et al. 2017). This loss function guides neural operators to align with PDEs, cutting down data needs. Moreover, specific architectures have been developed to maintain symmetries and conservation laws (Brandstetter et al. 20022] Liu et al. 2023), enhancing both generalization and data efficiency. Yet, these improvements largely focus on neural operators without fundamentally revising data generation.
Concurrently, there's push to create data-driven PDE solvers. For example, [Hsieh et al. (2019) suggests data-optimized parameterized iterative scheme. Additionally, combining neural operators with traditional numerical solvers is an emerging trend. An example is the hybrid iterative numerical transferable solver, merging numerical methods with neural operators for greater accuracy (Zhang) itt al. 2022). This model incorporates Temporal Stencil Modeling. Notably, numerical iterations can also be used as network layers, aiding in solving PDE-constrained optimization and other challenges.
2.2 KRYLOV SUBSPACE RECYCLING
The linear systems produced by the process of generating NOs training data typically exhibit prop- erties of sparsity and large-scale nature (Hao et al.) 2022 Zhang et al. 2023). In general, for non- symmetric matrices, Generalized minimal residual method (GMRES) (Saad & Schultz, 1986 Qin &
Xu 2023 of the Krylov algorithm is often employed for generating training data pertinent to NOs.
The concept of recycling in Krylov algorithms has found broad applications across multiple disci- plines. For instance. the technigue has been leveraged to refine various matrix algorithms (Wang al. 2007 Mehrmann Schröder, 2011). Furthermore, it has been deployed for diverse matrix problems (Parks et al. [006]] Gaul 2014) Soodhalter et al. (020)). Moreover, the recycling method- ology has also been adapted to hasten iterative solutions for nonlinear equations (Kelley) [2033] Deuflhard 2005) Gaul & Schlömer 2012 Gaul & Schlömer). Additionally, it has been instrumental in accelerating the resolution of time-dependent PDE equations (Meurant 2001; Bertaccini et al.) 2004 [Birken et al.] 2008). We have designed a sorting algorithm tailored for NO, building upon the original recycling algorithm, making more suitable for the generation of NO datasets.
3 PRELIMINARIES
3.1 DISCRETIZATION OF PDES IN NEURAL OPERATOR TRAINING
We primarily focus on PDE NOs for which the time overhead of generating training data is sub- stantial. As illustrated in Figure [] the generation of this training data necessitates solving the corresponding PDEs. Due to the inherent complexity of these PDEs and their intricate boundary conditions, discretized numerical algorithms like FDM, FEM, and FVM are commonly employed
for their solution (Strikwerda) 2004; (ughes 2012 Johnson 2012
Le Veque
2002 Cheng & Xu
2023
These discretized numerical algorithms embed the PDE problem from an infmite-dimennsionll Hilbert function space into suitable finite-dimensional space, thereby transforming the PDE is- sue into a system of linear equations. We provide a straightforward example to elucidate the process under discussion. The detailed procedure for generating the linear equation system available the Appendix A Specifically, we discuss solving two-dimensional Poisson equation using the FDM to transform it into a system of linear equation:
uu( v) = f(x,
using a 2 2 internal grid (i.e., N₂ = N2, = 2 and Ar = Ay), the unknowns ut,, can be arranged in row-major order follows: 41,1, 21,2; 12,1, 22 For central differencing on 2 x grid, the vector
3

---

Published as a conference paper at ICLR 2024
b will contain the values of f2 y/)
and the linear equation system Ax b can be expressed as:
-4
1
0
[1,1
f((11, y1)
4
0
U1,2
f(r1 y2)
1
0
4
1
12,1
f(r2, y1)
1
[2,,2]
[f y2)
By employing various parameters to generate f. such utiliziig Gaussian random fields (GRF) or truncated Chebyshev polynomials, we can derive Poisson equations characterized by distinct parameters.
Typically, training a NO requires generating between 10 and 10 numerical solutions for PDEs Lu al., 2019). Such multitude of linear systems, derived from the same distribution of PDEs, nat- urally exhibit a high degree of similarity. It is precisely this similarity that is key to the effective acceleration of our SKR algorithm's recycle module. (Soodhalter et al. 2020). We conceptu- alize this as the task of solving sequential series of linear equations:
A (1) g)
(())
i= 1,,,,
(I)
where the matrix A ( E C and the vector (+) E C vary based on different PDEs.
3.2 KRYLOV SUBSPACE METHOD
For solving large-scale sparse linear systems, Krylov subspace methods are typically used (Saad, 2003 Greenbaum 1997). The core principle behind this technique is leveraging the matrix from the linear equation system to produce lower-dimensional subspace, which aptly approximates the true space. This approach enables the iterative solutions within this subspace to gravitate towards the actual solution x.
Suppose we now solve the i-th system and remove the superscript of A The m-th Krylov sub- space associated with the matrix A and the starting vector r as follows:
Kmn ( (A, r) span(r, Ar, A²r, rr.
Typically, r chosen as the initial residual in linear system problem. The Krylov subspace iteration process leads the Arnoldi relation:
AV.. H
where Vm E C and H E C(n+T)xm is upper Hessenberg. Let Hm. E Cmxm denote the first m rows of Hm By solving the linear equation system related to Hm. we obtain the solution within the Km (A, r), thereby approximating the true solution. If the accuracy threshold is not met, we'll expand the dimension of the Krylov subspace. The process will iteratively transitions from Km,(A, (A,r) r) to K (A,r), continuing until the desired accuracy is achieved (Arnoldi, 1951 Saad 2003]. For the final Krylov subspace, we omit the subscript m and denote it simply as K(A,r)
For any subspace S C C", E S is Ritz vector of A with Ritz value 0, if
Ay Oy 1 w
Vi E
For iterative methods within the Krylov subspace, choose S Km (A,( r) and the eigenvalues of Hm are the Ritz values of A.
Various Krylov algorithms exist for difffrrnnt types of matrices. The linear systems generated by the NOs we discuss are typically non-symmetric. The most widely applied and successful algorithm for solving large-scale sparse non-symmetric systems of linear equations is GMRES (Saad & Schultz, 1986 Morgan 2002), which serves as the baseline for this paper.
METHOD
As shown in the Figure 2 we aim to accelerate the solution by harnessing the intrinsic correla- tion among these linear systems. To elaborate, when consecutive linear systems exhibit significant correlation, perhaps due to minor perturbations, it's plausible that certain information from the prior solution, such Ritz values and Ritz vectors, might not be fully utilized. This untapped information
4

---

Published as a conference paper at ICLR 2024
Data set
PDE generator
Solution
Solution2
Solution3
Linear systems
2
x
Krylov algorithm
Krylov algorithm
Krylov algorithm
Sort
SKR
(1)
Wz
Recycle
Recycle
d1
Krylov algorithm
d2
Krylov algorithm
d5
Krylov algorithm
Solution1
Solution2
Solution3
Data set
Figure 2: Algorithm Flow Diagram: a. Derive the PDEs to be solved from the NO. b. Transform these PDEs into system of linear equations. c. Apply the SKR algorithm to sort the linear systems, obtaining sequence with strong correlations. d. Traditional algorithms independently solve the linear systems from step b using Krylov methods. d1, d2, d3. The SKR algorithm utilizes 'recycling' sequentially solve the linear systems, reducing the dimension of the Krylov subspace. Obtain the solutions and assemble them into dataset.
could potentially expedite the resolution of subsequent linear systems. This concept, often termed 'recycling', has been well-established in linear system algorithms.
While the notion of utilizing the recycling concept is appealing. the PDEs and linear systems de- rived from NO parameters are intricate and convoluted. This complexity doesn't 'tt guarantee strong correlation between sequential linear equations, therehy limiting the effectiveness of this approach. To address this issue, we'v 'e formulated sorting algorithm tailored for the PDE solver. Our algo- rithm SKR, engineered for efficiency, ensures sequence of linear equations where preceding and subsequent systems manifest pertinent correlations. This organized sequence enables the fruitful application of the recycling concept to expedite the solution process.
4.1 THE SORTING ALGORRIHM
For linear systems derived from PDE solvers, the generation of these equations is influenced by specific parameters. These parameters and their associated linear equations exhibit continuous vari- ations (Lu et al. 2022, Li et al. 2020). Typically, these parameters are represented as parameter matrices. Our objective is to utilize these parameters assess the resemblance among various linear systems, and then sequence them to strengthen the coherence between successive systems.
Subsequent to in-depth Theoretical Analysis 5.2 arrive at the conclusion that the chosen recycling component largely resilient to matrix perturbations. It means there no necessity for excessive computational resources for sorting. As result, incorporated Sorting Algorithm based on the greedy algorithm, serving as the step e in Figure
In typical training scenarios for NOs, the number of data points, represented as num, be quite vast, frequently falling between 10 and 10a. To manage this efficiently, one could adopt a cost- effective sorting strategy. Initially, divide the data points into smaller groups, each containing either
5

---

Published as a conference paper at ICLR 2024
Algorithm 1: The Sorting Algorithm
Input: Sequence of linear systems to be solved A E C b(6) E C", corresponding
parameter matrix p()) E CPxP and i 1, N
Output:
Sequence for solving systems of linear equations selmat
Initialize the list with sequence seqo 2,, , seqmat is an empty list;
2 Set io = 1 as the starting point. And remove 1 from sego and append I to segmati do Refresh dis and set it to a large number, e.g., 1000; for each j in sego do dis, = the Frobenius norm of the difference between p(to and p();
for i 1,.......
if diss then dis diss jmin = j:
Remove Jmin from sequ and append jmin seq and set ig Imini
Get the sequence for solving linear systems seqmat-
10 to 10 data points, based their coordinates. Then, use the greedy algorithm to sort within these groups Once sorted, these smaller groups can be concatenated.
4.2 KRYLOV SUBSPACE RECYCLING
For sequence of linear systems, there exists inherent correlation between successive systems. Consequently, we posit that by leveraging information from the solution process of prior systems, can expedite the iterative convergence of subsequent system, thereby achieving notable accel- eration in computational performance.
For various types PDE problems, the linear systems generated exhibit matrices with distinct struc- tural characteristics. These unique matrix structures align with specific Krylov recycling algorithms. To validate our algorithm's efficacy, we focus primarily on the most prevalent and general scenario within the domain of NOs, where matrix A is nonsymmetric. Given this context, Generalized Con- jugate Residual with Deflated Restarting (GCRO-DR stands out as the optimal Krylov recycling algorithm. It also forms the core algorithm for steps (d1, d2, d3) in Figure 2
GCRO-DR employs deflated restarting, building upon the foundational structure GCRO as illus- trated in (de Sturler 1996). Essentially, GCRO-DR integrates the approximate eigenvector derived from the solution of a preceding linear system into subsequent linear systems. This is achieved through the incorporation of a deflation space preconditioning group solution, leading to a pro- nounced acceleration effect (Soodhalter et al. 2020; Nicolaides
1987 Erlangga & Nabben 2008
The detailed procedure is delineated as follows:
Suppose that we have solved the system of I with GCRO-DR, and we retain approximate eigenvectors, Yx [ .... ÿt}. Then, GCRO-DR computes matrices Uv. C₂ E C from Y₁ and A (i+1) such that A 1 U2 C₁ and C C I.. The specific algorithm can be found in Appendix B.1 By removing the superscript of A (i+1) we get the Amoldi relation of GCRO-DR:
(I ).C.CH)AAA
Vn H
Due the correlation between the two consecutive linear systems, Cj encompasses informa- tion from the previous system A (+), This means that when constructing new Krylov subspace ((A (i+1), ((+1)),, there is no need start from scratch; the subspace span Yy3 composed of
some eigenvectors from (+) already exists. Building upon this foundation, KC K(A (++1), y (f+1)) can converge more rapidly to the subspace where the solution lics. This can significantly reduce the dimensionality of the final Krylov subspace, leading to marked decrease in the number of iterations and resulting in accelerated performance.
Compared to the Arnoldi iteration presented in 3 used by GMRES, the way GCRO-DR leverages solution information from previous systems accelerate the convergence of linear sequence be- comes clear. GMRES can be intuitively conceptualized as the special case of GCRO-DR where is initialized at zero (Carvalho al. 2011 Morgan 2002 Parks et al. 2006).

---

Published as a conference paper at ICLR 2024
The effectiveness of acceleration within the Krylov subspace recycling is highly dependent on lecting Yú in manner that boosts convergence for the following iterations of the linear systems. Proper sequencing can amplify the impact of Yx, thereby reducing the number of iterations. This underlines the importance of sorting. For detailed understanding of GCRO-DR, kindly consult the pseudocode provided in the Appendix B.2
THEORETICAL ANALYSIS
5.1
CONVERGENCE ANALYSIS
When solving a single linear system, GCRO-DR and GMRES-DD DR are algebraically equivalent. The primary advantage of GCRO-DR is its capability for solving sequences of linear systems (Carvalho et al., 2011: Morgan 2002 Parks et al. 2006). Stewart (2001): Embree (2022) provide a good framework analyze the convergence of this type of Krylov algorithm. Here we quote Theorem 3.1 in Parks al. (2006) to make detailed analysis of GCRO-DR. We define the one sided distance from the subspace Q to the subspace as
8(2.C) = |I(I - Tee Tell,,
(5)
where II represents the projection operator for the associated space. Its mathematical interpretation corresponds to the largest principal angle between the two subspaces Q and (Beattie al. 2004).
and defineing Pe as the spectral projector onto Q.
Theorem 1. Given space C range(Ca,) , let V range m Hm be the (m k) dimen- sional Krylov subspace generated by GCRO-DR as in 44) Let To E C", and let Tı = r.
Then, for each Q such that 6(2,C) 1,
min dievec
||ro
2 1
min dsE(1-Po)))
I - Po)ri - d2l12
(6)
1 Pollz ||(I Ily)ril12
where 7 = - Pell
In the aforementioned bounds, the left-hand side signifies the residual norm subsequent to m- k iter- ations of GCRO-DR, employing the recycled subspace C. Contrastingly, 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 (Morgan 2002, Simoncini Szyld 2005 Van der Vorst Vuik) 1993). The subsequent term on the right embodies a constant multiplied by the residual follow- ing m 1 iterations of GCRO-DR, when solving for T1- Should the recycling space encompass an invariant subspace Q, then = y 0 for the given Q, ensuring that the convergence rate of GCRO-DR matches or surpasses that the deflated problem. In most cases, |PPll is numerically stable and not large, therefore reduced value of directly correlates with faster convergence in GCRO-DR.
5.2 THE RATIONALITY OF SORTING ALGORITHMS
critical inquiry arises: How much computational effort should be allocated sorting algorithm? Given the plethora of sorting algorithms available, is it necessary to expend significant computational resources identify the optimal sorting outcome?
"When the magnitude of the change smaller than the gap between smallest and large eigenvalues, then the invariant sub- space associated with the smallest cigenvalues is not significantly altered." We derive the following insights: The perturbation of invariant subspaces associated with the smallest eigenvalues when the change in the matrix is concentrated in an invariant subspace corresponding to large eigenvalues (Si- fuentes et al., 2013 Parks et al., (006). This effectively implies that the pursuit of optimal sorting isn't imperative. Because we typically recycle invariant subspaces formed by smaller eigenvec- tors. As long as discernible correlations persist between consecutive linear equations, appreciable acceleration effects attainable (Gaul 2014). This rationale underpins decision to develop a cost-efficient sorting algorithm. Despite its suboptimal nature, the resultant algorithm exhibits commendable performance.
Drawing from Theorem 3.2 as presented in (Kilmer & De Sturler 2006):
7

---

Published as a conference paper at ICLR 2024
EXPERIMENT
6.1
SET UP
To comprehensively evaluate the performance of SKR in comparison to another algorithm, we con- ducted nearly 3,000 experiments. The detailed data is available in the Appendix D.6 Each experi- ment utilized a data set crafted to emulate an authentic NO training data set. Our analysis centered on two primary performance metrics viewed through three perspectives. These tests spanned four different datasets, with SKR consistently delivering commendable results. Specifically, the three Perspectives are: Matrix preconditioning techniques, spanning 7 to 10 standard methods. 2. Ac- curacy criteria for linear system solutions, emphasizing 5 to 8 distinct tolerances. 3. Different matrix sizes, considering computational time overhead. 2. Mean iteration count. For a deeper dive into the specifics, please consult the Appendix D.I
variations. Our primary performance Metrics encompassed: I. Average
Baselines. As previously alluded to, our focus revolves around system of linear equations, which consists of large sparse non- symmmetttii matrices. The GMRES algorithm serves as the predominant solution and sets the benchmark for our study. We utilized the latest version from PETSc 3.19.4 for GMRES.
Datasets. To probe the algorithm's adaptability across matrix types, we delved into four distinct linear equation challenges, each rooted in PDE: Darcy Flow Problem (Li al.) [2020] Rahman
al. 2022 Kovachki al. 202 Lu et al. 2022); 2. Thermal Problem (Sharma al. 2018 Koric
Abueidda 2023); 3. Poisson Equation (Hsieh et al.) 2019 Zhang et al., y022); Helmholtz
Equation (Zhang al. 2022). For an in-depth exposition of the dataset and its generation, kindly
refer to the Appendix D.2, For the runtime environment, refer to Appendix D.4] D.4]
6.2 QUANTITATIVE RESULTS
Table 1: Comparison of our SKR and GMRES computation time and iterations across datasets, preconditioning, and tolerances. The first column lists datasets with matrix side lengths, the next details tolerances. The data is displayed "computation time speed-up ratio/iteration count speed- up ratio". GMRES/SKR ratio over denotes better SKR performance.
Dataset
Time/her
None
Jacobi
BJacobi
SOR
ASM
ICC
ILU
le-2
2.62/19.2
2.88/22.6
3.21/23.6
2.69/23.3
2.23/13.4
1.97/9.55
1.93/9.44
Darcy
le-5
2.92/21.1
3.42/24.5
4.07/28.6
3.45/27.9
3.66/22.2
3.18/14.9
3.08/14.4
6400
le-8
2.70/19.1
3.09/22.9
4.00/27.5
3.54/27.5
4.53/25.8
4.11/18.9
3.70/17.2
le-5
4.53/20.8
3.32/15.0
2,38/10.3
1.96/8.76
2.46/10.3
2.40/10.3
2.35/10.3
Thermal
le-8
5,35/23.6
3.06/13.5
2.73/11.1
2.34/9.30
2.83/11.1
2.77/11.1
2.69/11.1
11063
lc-11
5.47/24.9
2.93/12.8
3.05/11.7
2.62/10.2
3.14/11.7
3.08/11.7
3.01/11.7
le-5
1.27/4.69
1,28/4.68
1.13/3887
1.19/4.05
1.46/3.93
0.99/3.92
0.98/3.92
Poisson 71313
lc-8
1.74/6.29
1.75/6.30
1.19/3.97
1.35/4.45
1.94/4.83
1.32/4,83
1.30/4.87
lc-11
1.90/6.83
1.91/6.82
1.19/3.95
1,33/4.35
2.12/5.11
1.42/5.13
1.38/5.13
le-2
7.74/17.3
8.44/20.3
8.83/21.5
8.37/22.3
10.6/25.4
4.77/21.4
3,88/17.6
Helmholtz
le-5
7.61/16.5
8.62/20.0
11.4/26.5
10.5/26.2
13.1/29.3
6.38/28.3
6.13/27.2
10000
1o-7
6.47/13.68
7.96/18.1
11.3/26.2
10.6/25.9
13.9/30.0
6.72/29.3
6.34/28.0
Table [holccssssassee selected experimental data. From this table, we can infer several conclusions:
Firstly, across almost all accuracy levels and preconditioning techniques, our SKR method consis- tently delivers impressive acceleration. It is most noticeable in the Helmholtz dataset. Depending on the convergence accuracy and preconditioning technique, our method reduces wall clock time by factors of 2 to 14 and requires up to 30 times fewer iterations. This suggests that linear sys- tems derived from similar PDEs have inherent redundancies in their resolution processes. Our SKR algorithm effectively minimizes such redundant computations, slashing iteration counts and signifi- cantly speeding up the solution process, ultimately accelerating the generation of datasets for neural operators.
Secondly, in most cases, as the demanded solution accuracy intensifies, the acceleration capability of algorithm becomes pronounced, For instance, in the Darey flow problem, the acceleration
8

---

Published as a conference paper at ICLR 2024
for high accuracy level (le-8) is 50% 100% more than for lower accuracy (1e-2). On one hand, high- segh-preciiio data demands larger Krylov subspace and more iterations. SKR's s recycling tech- nology leverages information from previously similar linear systems, enabling faster convergence in the Krylov subspace. Thus, its performance is especially superior at higher precisions. Detailed subsequent analyses further corroborate this observation, as seen in Appendix [D.5.1 D.5.2 On the other hand, Theoretical Analysis 5.2 indicates that the information recycled in the SKR algorithm possesses strong anti-perturbation capabilities, ensuring maintained precision. Its acceleration effect becomes even more evident at higher accuracy requirements.
Thirdly, we observe that our SKR algorithm can significantly reduce the number of iterations, with reductions up to 30 times, indirectly highlighting the strong algorithmic stability of the SKR method. By analyzing the count of experiments reaching the maximum number of iterations, we can conclude that the stability of the SKR algorithm is superior to that of GMRES. Detailed experimental analyses can be found in the Appendix DD5.3
Lastly, it is worth noting that outcomes vary based on the preconditioning methods applied. Our algorithm demonstrates remarkable performance when no preconditioner, Jacobian, or Successive Overrelaxation Preconditioner (SOR) methods are used. However, its efficiency becomes less prominent for low-accuracy Poisson equations that employ incomplete Cholesky factorization (ICC) incomplete lower-upper factorization (ILU) preconditioners. The underlying reason is that al- gorithms like ICC and ILU naturally operate by omitting smaller matrix elements or those with negligible impact on the preconditioning decomposition, and then proceed with Cholesky or LU de- composition. Such an approach somewhat clashes with the recycling method. While SKR leverages information from prior solutions of similar linear equation systems, ICC and ILU can disrupt this similarity. As result, when employing ICC or ILU preconditioning, the improvement in SKR's speed less noticeable lower accuracy levels. However, at higher precisions, our SKR algorithm still exhibits sigmificant acceleration benefits with ICC and ILU preconditioning
6.3 ABLATION STUDY
To assess the impact of 'sort' on SKR algorithm, we approached the analysis from both theoretical and exper- imental perspectives:
Table 2: Comparison of algorithm performance with and without sort in Darcy flow problem, using SOR precon- ditioning, matrix size 10 and compu- tational tolerance le
Theoretical Perspective: Based on previous Theoretical insights 5.1 the pivotal metric to gauge the convergence speed of the Krylov recycle algorithm is 8. Hence, we ex- amined the variations in the ő metric before and after sort- ing. Experimental Perspective: We tested the SKR algo- rithm both with and without the 'sort' feature, examining the disparities in computation time and iteration count.
Time(s)
Iter
SKR(sort) SKR(nosort)
0.101
183.9
0.90
0.114
202.5
0.95
As illustrated in the Table 2 The metric declined by 5% after employing the 'sort' algorithm. As illustrated in Figure 8 this effectively demonstrates that the 'sort' algorithm can enhance the correlation between consecutive linear equation sets, achieving the initial design goal of 'sort.. 2. Using 'sort' enhances the SKR algorithm's 's computational speed by roughly 13% and decreases its iterations by 9.2%. This implies that by increasing the coherence among sequential linear equation sets, the number of iterations needed is minimized, thus hastening the system's solution process.
7 LIMITATION AND CONCLUSIONS
Our paper presents the SKR algorithm as a significant advancement in generating neural operator datasets, especially for linear PDEs, yet recognizes areas for improvement. Future work could explore adapting the SKR for non-linear PDEs, refining the sorting algorithm's 's distance metrics for better system correlation, and extending the recycling concept to further data generation realms. This pioncering effort not only enhances computational efficiency but also makes neural operators more accessible, indicating a major leap in their development and application.

---

Published as a conference paper at ICLR 2024
REPRODUCIBILITY STATEMENT
For the sake of reproducibility, we have included our codes within the supplementary materials. However, it's worth noting that the current code version lacks structured organization. Should this paper be accepted, we commit to reorganizing the codes for improved clarity. Additionally, in the Appendix D we provide in-depth description of our experimental setups and detailed results.
ACKNOWLEDGEMENTS
The authors would like 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 0022ZD0119801 National Nature Science Foundations of China grants U23A20388, U19B2026, U19B2044, and 62021001.
REFERENCES
Walter Edwin Arnoldi. The principle of minimized iterations in the solution of the matrix eigenvalue problem. Quarterly of applied mathematics, 91):177-9, 1951.
Yannick Augenstein, Taavi Repan, and Carsten Rockstuhl. Neural operator-based surrogate solver for free-form electromagnetic inverse design. ACS Photonics, 2023.
Christopher Beattie, Mark Embree, and John Rossi. Convergence of restarted krylov subspaces to invariant subspaces. SIAM Journal Malrix Analysis and Applications, 25(4):1074-1109 2004.
Michele Benzi, Carl D Meyer, and Miroslav Tma. A sparse approximate inverse preconditioner for the conjugate gradient method. SIAM Journal on Scientific Computing, 17(5):1 1996.
Daniele Bertaccini et al. Efficient preconditioning for sequences of parametric complex symmetric linear systems. Electronic Transactions on Numerical Analysis, 18:49-64, 2004.
Philipp Birken, Jurjen Duintjer Tebbens, Andreas Meister, and Miroslav Tma. Preconditioner up- dates applied to cfd model problems. Applied Numerical Mathematics, 58(1 (8(11))1628-1641 2008
Johannes Brandstetter, Rianne den Berg, Max Welling, and Jayesh K Gupta. Clifford neural layers for pde modeling. arXiv preprint wwXi:2220.044933 2022.
Luiz Mariano Carvalho, Serge Gratton, Rafael Lago, and Xavier Vasseur. A flexible generalized conjugate residual method with inner orthogonaliizaton and deflated restarting. SIAM Journal on Matrix Analysis and Applications, 32(4):1212-1235, 2011.
Lu Cheng and Kuan Xu. Solving time-dependent pdes with the ultraspherical spectral method, 2023.
Philippe G Ciarlet. The finite element method for elliptic problems. SIAM, 2002.
Eric de Sturler. Nested krylov methods based ger. Journal of Computational and Applied Math-
ematics, 67(1):15-41, 1996.
Eric De Sturler. Truncation strategies for optimal krylov subspace methods. SIAM Journal Numerical Analysis, 36(3):864-889, 1999.
Peter Deuflhard. Newton methods for nonlinear problems: affine invariance and adaptive algo- rithms, volume 35. Springer Science & Business Media, 2005.
Tobin A Driscoll, Nicholas Hale, and Lloyd N Trefethen. Chebfun guide, 2014.
Mark Embree. How descriptive are gmres convergence bounds? arXiv preprint arXiv:2209.01231, 2022.
Yogi A Erlangga and Reinhard Nabben. Deflation and balancing preconditioners for krylov subspace methods applied to nonsymmetric matrices. SIAM Journal on Matrix Analysis and Applications, 30(2):684-699 2008.
10

---

Published as a conference paper at ICLR 2024
A Gaul and N Schlömer, pynosh: Python framework for nonlinear schrödinger equations. july 2013. URI: https:/bbituckkt. org/nschloe/pynosh.
André Gaul. Recycling krylov subspace methods for sequences of linear systems. 2014.
André Gaul and Nico Schlömer. Preconditioned recycling krylov subspace methods for self-adjoint problems. arXiv preprint arXiv: /208.. 02644, 2012.
André Gaul, Martin H Gutknecht, Jorg Liesen, and Reinhard Nabben. A framework for deflated and augmented krylov subspace methods. SIAM Journal Matrix Analysis and Applications, 34(2): 495-518, 2013.
Zijie Geng, Xijun Li, Jie Wang, Xiao Li, Yongdong Zhang, and Feng Wu. A deep in- stance generative framework for milp solvers under limited data availability. In A. Oh, T. Neumann, A. Globerson, K. Saenko. M. Hardt, and S. Levine (eds.). Advances in Neu- ral Information Processing Systems, volume 36, pp. 26025 00025-26047, Curran Associates, Inc., 2023. URL pttps:://proceed.ngg. #LL:://9/eSSacbb2222ftteeteeeee-
http://oocceddd..............oee neurips cc cc/ pwwer_appprrpeeer httpc/paarr__feeppaapppeeeeee htt.fffll/e/ppp//22/33333
Paper Conference.
Anne Greenbaum. Iterative methods for solving linear systems. SIAM, 1997.
Daniel Greenfeld, Meiray Galun, Ronen Basri, Irad Yavneh, and Ron Kimmel. Learning to optimize multigrid pde solvers. In International Conference on Machine Learning, pp. 2415-2423. PMLR, 2019.
Arne S Gullerud and Robert H Dodds Jr. Mpi-based implementation of a peg solver using ebe architecture and preconditioner implicit, 3-d finite element analysis. Computers Structures, 79(5):553-573 2001.
Zhongkai Hao, Songming Liu, Yichi Zhang, Chengyang Ying, Yao Feng, Hang Su, and Jun Zhu. Physics-informed machine learning: A survey on problems, methods and applications. arXiv preprint arXiv:2211.0800 10.0644 2022.
Jun-Ting Hsieh, Shengjia Zhao, Stephan Eismann, Lucia Mirabella, and Stefano Ermon, Learning neural pde solvers with convergence guarantees. arXiv preprint 190..0000000000, 2019.
Thomas JR Hughes. The finite element method: linear static and dynamic finite element analysis. Courier Corporation, 2012.
Claes Johnson. Numerical solution of partial differential equations by the finite element method. Courier Corporation, 2012.
Carl T Kelley. Solving nonlinear equations with Newton's method. SIAM, 2003.
Misha E Kilmer and Eric De Sturler. Recycling subspace information for diffuse optical tomography.
SIAM Journal on Scientific Computing. 27(6):2140-211 22662216.6220 2006.
Dana A Knoll and David E Keyes. Jacobian-free newton-kryloo methods: survey of approaches and applications. Journal of Computational Physics, 193(2):357-397, 2004.
Seid Koric and Diab W Abueidda. Data-driven and physics-informed deep learning operators for solution of heat conduction equation with parametric heat source. International Journal of Heat and Mass Transfer, 203:123809, 2023.
Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Avizzadenesheli, Kaushik Bhattacharya, An- drew Stuart, and Anima Anandkumar. Neural operator: Learning maps between function spaces. arXiv preprint arXiv:2108.0844 2021.
Randall J Vequee Finite volume methods for hyperbolic problems, volume 31. Cambridge univer- sity press, 2002.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, An- drew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. arXiv preprint 100.0220 2020.

---

Published as a conference paper at ICLR 2024
Chih-Jen Lin and Jorge J Moré. Incomplete cholesky factorizations with limited memory. SIAM Journal on Scientific Computing, 21(1):24-45 45, 1999.
Ning Liu, Yue Yu, Huaiqian You, and Neeraj Tatikola. Ino: Invariant neural operators for learning complex physical systems with momentum conservation. In International Conference on Artifi- cial Intelligence and Statistics, pp- 6822-6838. PMLR, 2023.
Lu Lu, Pengzhan Jin, and George Em Karniadakis. Deeponet: Learning nonlinear operators for iden- tifying differential equations based the universal approximation theorem of operators. arXiv preprint arXiv:1910.0311 2019.
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 Applied Mechanics and Engineer- ing. 393: 14778, 2022.
Volker Mehrmann and Christian Schröder. Nonlinear eigenvalue and frequency response problems in industrial practice. Journal of Mathematics in Industry, 1:1-18, 2011.
Gérard Meurant. On the incomplete cholesky decomposition of class of perturbed matrices. SIAM Journal on Scientific Computing, 23(2):419-429, 2001.
Igor Moret. note on the superlinear convergence of gmres. SIAM journal on numerical analysis, 34(2):513-516, 1997.
Ronald B Morgan. Gmres with deflated restarting. SIAM Journal on Scientific Computing, 24(1): 20-37, 2002.
Shin Muroya, Atsushi Nakamura, Chiho Nonaka, and Tetsuya Takaishi. Lattice qed at finite density: introductory review. Progress of theoretical physics, 110(4):615-66) 668, 2003.
Roy A Nicolaides. Deflation of conjugate gradients with applications to boundary value problems. SIAM Journal on Numerical Analysis, 24(2):355-365 1987.
Michael L. Parks, Eric De Sturler, Greg Mackey, Duane Johnson, and Spandan Maiti. Recycling krylov subspaces for sequences of linear systems. SIAM Journal Scientific Computing, 28(5): 1651-1674, 2006.
Jaideep Pathak, Shashank Subramanian, Peter Harrington, Sanjeev Raja, Ashesh Chattopadhyay. Morteza Mardani, Thorsten Kurth, David Hall, Zongyi Li, Kamyar Azizzadenesheli, et al. Four- castnet: A global data-driven high-resolution weather model using adaptive fourier neural opera- tors. arXiv preprint arXiv: 11214, 2022.
Ouyuan Qin and Kuan Xu. Solving nonlinear odes with the ultraspherical spectral method. arXiv preprint 10X 17688, 2023.
Md Ashiqur Rahman, Zachary Ross, and Kamyar Azizzadeneshelii U-no: U-shaped neural oper- ators. arXiv preprint arXiv:2204.11122 2022.
Maziar Raissi, Paris Perdikaris, and George Em Karniadakis. Physics informed deep learn- ing (part i): Data-driven solutions of nonlinear partial differential equations. arXiv preprint arXiv:1711. 10561, 2017.
Youcef Saad and Martin H Schultz. Gmres: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM Journal on scientific and statistical computing. 7(3):856- 869. 1986.
Yousef Saad. Iterative methods for sparse linear systems. SIAM, 2003.
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 1807. 13374, 2018.
12

---

Published as a conference paper at ICLR 2024
Josef A Sifuentes, Mark Embree, and Ronald B Morgan. Gmres convergence for perturbed co- efficient matrices, with application to approximate deflation preconditioning. SIAM Journal on Matrix Analysis and Applications, 34(3):1066-1000 2013.
Valeria Simoncini and Daniel В Szyld. On the occurrence of superlinear convergence of exact and inexact krylov subspace methods. SIAM review, 47(2):247-272, 2005.
Gordon D Smith. Numerical solution of partial differential equations: finite difference methods.
Oxford university press, 1985.
Kirk M Soodhalter, Eric Sturler, and Misha E Kilmer. survey of subspace recycling iterative methods. GAMM-Mitteilungen, 43(4)::202000016, 2020.
Gilbert W Stewart. Matrix Algorithms: Volume II: Eigensystems. SIAM, 2001.
John Strikwerda. Finite difference schemes and partial differential equations. SIAM, 2004.
James William Thomas. Numerical partial differential equations: finite difference methods, vol- ume 22, Springer Science & Business Media, 2013.
Andrea Toselli and Olof Widlund. Domain decomposition methodss-algorihhms and theory, vol- ume 34. Springer Science & Business Media, 2004.
Henk A Van der Vorst and Vuik. The superlinear convergence behaviour of gmres. Journal of computational and applied mathematics, 48(3):327-341, 1993.
Shun Wang. Eric de Sturler, and Glaucio Paulino. Large-scale topology optimization using pre- conditioned krylov subspace methods with recycling. International journal for numerical meth- ods in engineering, (1112)2:44112468 2007.
Gege Wen, Zongyi Li, Kamyar Azizzadenesheli, Anima Anandkumar, and Sally M Benson. U- fno enhanced fourier neural operator-based deep-learning model for multiphase flow. Ad- vances in Water Resources, 163: 1000 2022.
David Young. Iterative methods for solving partial difference equations of elliptic type. Transactions of the American Mathematical Sociery, 76(1):92-111, 1954.
Eleftherios Zachmanoglou and Dale W Thoe. Introduction to partial differential equations with applications. Courier Corporation, 1986.
Enrui Zhang, Adar Kahana, Eli Turkel, Rishikesh Ranade, Jay Pathak, and George Em Karniadakis. A hybrid iterative numerical transferable solver (hints) for pdes based on deep operator network and relaxation methods. arXiv preprint arXiv:2208. 2tt.:22.8..327... 2022.
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, Liu, Yuanqi Du, Alexan- dra Saxton, Hongyi Ling, Hannah Lawrence, Hannes Stärk, Shurui Gui, Carl Edwards, Nicholas Gao, Adriama Ladera, Tailin Wu, Elyssa F. Hofgard, Aria Mansouri Tehrani, Rui Wang, Ameya Daigavane, Montgomery Bohde, Jerry Kurtin, Qian Huang, Tuong Phung, Minkai Xu, Chai- tanya K. Joshi, Simon Mathis, Kamyar Azizzadeneshelii Ada Fang, Alán Aspuru-Guzik, Erik Bekkers, Michael Bronstein, Marinka Zitnik, Anima Anandkumar, Stefano Ermon, Pietro Liø, Rose Yu, Stephan Günnemann, Jure Leskovec, Heng Ji, Jimeng Sun, Regina Barzilay, Tommi Jaakkola, Connor W. Coley, Xiaoning Qian, Xiaofeng Qian, Tess Smidt, and Shuiwang Ji. Ar- tificial intelligence for science in quantum, atomistic, and continuum systems. arXiv preprint arXiv:2307.08423. 2023.
13

---

Published as a conference paper at ICLR 2024
A CONVERTING PDES TO LINEAR SYSTEMS: AN EXAMPLE
A.1 OVERVIEW
The general approach for solving Partial Differential Equations (PDEs) using techniques like Finite Difference Method (FDM), Finite Element Method (FEM), and Finite Volume Method (FVM) can be broken down into the following key steps (Strikwerda 2004 Hughes
2012.
Johnson 2012
Veque 2002):
1. Mesh Generation: Partition the domain over which the PDE is defined into grid of specific shapes, such as squares or triangles. As illustrated, the Figure shows the FEM finite element mesh for the Poisson equation with square boundary.
Block With Finite Element Mesh Displayed
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0.2
0.4
0.6
0.8
Figure 3: FEM mesh for the Poisson equation
2. Equation Discretization: Utilize the differential form of operators or select basis functions com- patible with the grid to transform the PDE into a discrete problem. Essentially, this maps the PDE from an infinte-dimmeniionll Hilbert space to a finite-dimensional representation.
Matrix Assembly: If dealing with linear PDEs, the discretized PDE and boundary conditions can be converted into system of linear equations in the fiite-dimeenioonl space. The system assembled based on the chosen grid type. For nonlinear PDEs, methods similar to Newton's iteration typically employed, transforming the problem into series of linear systems to be solved iteratively (Ciarlet 2002 Knoll & Keyes 2004).
4. Applying Boundary Conditions: Discretize the boundary conditions the grid and incorporate them into the linear systems.
5. Solving the System of Linear Equations: This step is generally the most time-consuming part of the entire algorithm.
6. Obtaining the Numerical Solution: Map the solution of the system of linear equations back onto the domain of the original PDE using the specific grid and corresponding basis functions, thereby obtaining the numerical solution.
A.2 SPECIFIC EXAMPLES
To illustrate how the FDM can transform the Poisson equation into a system of linear equations, let's consider concrete and straightforward example (Smith) 1985). Assume aim solve the Poisson equation in a two-dimensional space:
Vuu(,, v) fxx, v)-
14

---

Published as a conference paper at ICLR 2024
This equation is defined square domain (a, b] {c, d] with given boundary conditions. The process can be broken down into the following steps:
1. Mesh Generation: In the x-direction, select points 2; y-direction, select points yj = c jAy, where j = ,, Ny.
iAz, where = 0, 1, ................................ No. In the
2. Equation Discretization: Discretize the Poisson equation using central difference scheme. For a grid point (zi, Vj), the discretized equation becomes:
1,j. 2 1i + Wi.+1 2ui.j
Vi fū.j
Ar
Ay²
Here, May is approximation of uus;, Wb), and Jig [3 Uj)-
This equation be rewritten
1
1 Ay² fa.j
-2u4
11 A. + u4 Ay² 2a
Ax
Ay²
A
Matrix Assembly: Each equation corresponding to an internal point (i.j) contributes one row to the matrix A and the vector b. The corresponding row in A contains all the coefficients for Wijj while the element is si..
4. Applying Boundary Conditions: Implement the given boundary conditions, which may be of Dirichlet type (specifying u values) or Neumann type (specifying the derivative of u). Adjust the corresponding elements in the vector b by adding or subtracting terms related to these boundary conditions.
5. Solving the System of Linear Equations Ax = b Due to the large, sparse nature of the matrix, iterative methods are generally used to solve the system of equations (Saad 2003 Greenbaum 1997
Obtaining the Numerical Solution: Based on the difference scheme, map the solution of the sys- tem of linear equations from the discrete mnite-dimensional space the function space of the PDE to obtain the numerical solution. this example, Ws, approximates u(x; Uj), so a interpolation step can be used finalize the solution.
B ALGORITHMIC DETAILS
B. COMPUTES MATRICES U₁ AND C
The following computational procedure is adapted from De Sturler (1999): Parks et al. (2006).
GCRO-DR can modified solve by carrying over Us from the i-th system the (i + 1)-th system. We have the relation A Ux Cs. We modify Us and Cs satisfy
A (0) U₁ = Ck:
C# C = I
with respect to (++1) follows:
Q R = reduced QR decomposition of A (lt1)
CRew Q,
meew Uold R
15

---

Published as a conference paper at ICLR 2024
B.2 GCRO-DR
The following computational procedure adapted from Parks et al. (2006).
Algorithm 2: GCRO-DR
Choose m, the maximum size of the subspace, and k. the desired number of approximate eigenvectors. Let tol be the convergence tolerance. Choose an initial guess Zp- Compute =
ro Avo. and
if Yk is defined (from solving previous linear system) then
Let [Q. R| be the reduced QR-factorization of AYg-
C = Q
U Y
æ1 = zu + U k H
TO
r1 ro C Hro
else
= ro/l|roll /|oolll
Perform m steps of GMRES, solving min ||e Hm for and generating Vn+1 and 10
c= =|roll2e1
H m
12
21 In V
13
Vm+i(c- mm
14
Compute the eigenvectors ž, of (Hm h hm+1,m H H C. MZ 8, associated with
the smallest magnitude eigenvalues f, and store in Pk. hm+l,m m is the element in m + and column m of matrix H
15
YE Vm Pk
16
Let [Q, R] be the reduced QR-factorization of H Px.
C₁ Vm+1Q
18
U = Y.
19
while |r.1l 2|2 > tol do
i +1
Perform
Amoldi steps with the linear operator (I C₁ CHC A, letting Vi 1.=r--//|rr... 11: and generating Vm-k+1, Hm-k and Bm
22
Let D; be diagonal scaling matrix such that U U DE. where the columns of Ug have unit norm.
23
m
1Ux
Vn
24
W m++1
[C
Vm-k+1]
k
B Bm-k
25
G
Hm-k
Solve min [WH m+ Ti-1 Gn for y.
27
2; 1; + Vn
28
r, = - Wm-1Gm3
Compute the eigenvectors ε, of G G 2, 0, H. WH.. m+ Vm associated with smallest
magnitude eigenvalues 0, and store in Pk
30
Y, V.P.. Pt
Let (Q, R] be the reduced QR-factorization of G,,, Pt.
C₂ Wm+1Q
33
U = Y.R
Let Y, Up (for the next system)
16

---

Published as a conference paper at ICLR 2024
C SUPPLEMENTARY THEORETICAL ANALYSIS
C.I SUPERLINEAR CONVERGENCE PHENOMENON
The phenomenon of superlinear convergence is freguently observed in the iterative processes of algorithms like GMRES (Van der Vorst & Vuik 1993; Moret 1997 Simoncini & Szyld 2005). The introduction of a deflation space enables these iterative algorithms to swiftly transition into the superlinear convergence phase, thereby expediting their performance (Gaul et al., 2013; Nicolaides) 1987]- The principle of recycling essentially embodies this concept, leveraging deflation space for acceleration. This phenomenon has been discerned in algorithms associated with recycling as well (Gaul) 2014).
For the problem under consideration, when sequence of linear systems manifests inherent corre- lations in their sequential order, an apt deflation space can be curated. This ensures rapid entry into the superlinear convergence phase, bypassing the initial slower convergence stages, and thereby enhancing the overall convergence velocity. Empirical data from our experiments underscores this observation. Theoretically, we can further dissect the potential reductions in iterations offered by SKR. The specific experimental phenomenon is illustrated D.5. D.5.2
D DETAILS OF EXPERIMENTAL DATA
D.1 SPECIFIC PARAMETERS OF THE MAIN EXPERIMENT
Baseline:
To ensure optimal speed, we first generated PDE linear equation systems using either Python or Matlab. Then solve them in the programming environment. We utilized the latest version from PETSc 3. 1..4 for GMRES.
Three Perspectives:
When solving large matrices, effective preconditioning can greatly accelerate the resolution process and enhance the stability of the algorithm. Recognizing that different scenarios require distinct preconditioning approaches, we conducted tests with more than ten of the most commonly used methods.
1. Precondition:
The tolerance of a solution inherently determines the iteration count and, by ex- tension, the time required for solution. Distinct algorithms exhibit varied convergence rates, and specific NOs hold unique tolerance standards. To ensure an encompassing comparison, we evalu- ated computational durations diverse tolerances, zeroing on 5-8 optimal error precisions for our tests.
Tolerance:
Matrix Dimensionality: Algorithmic performance tends fluctuate based the matrix's di- mensionality. Moreover, varying NOs mandate matrices of disparate dimensions. Hence, our study incorporated five distinct matrix dimensions.
1. Computational Duration: This serves the most unambiguous metric to gauge an algorithm's is efficacy. Iteration Count: Algorithmic stability is emblematic of its numerical sensitivity, with the iteration count offering direct insights into said stability.
Two Performance Metrics:
D.2 DATA SET
1. Darcy Flow
We consider two-dimensional Darey flows, which can be described by the following equation Lii
et al. 2020 Rahman et al.) 2022. Kovachki et al. 2021 Lu et al. 2022
(K(x,)) v) f.
where K is the permeability field, h is the pressure, and f is a source term which can be either a constant space-dependent function. In experiment, K (z. is derived using the Gaussian Random Field (GRF) method. The parameters inherent to the GRF serve as the foundation for our sort scheme.
17

---

Published as a conference paper at ICLR 2024
Figure 4: Solutions of Darcy flow equations with close parameters
Figure 5: Solutions of Darcy flow equations with divergent parameters
As illustrated, with the matrix two-norm serving as the metric for distance, the first the solution for the Darey Flow problem with two NO parameters that are very close. The subsequent Figure [ depicts solution for the Darcy Flow problem where the two NO parameters differ significantly. Clearly, when the parameters are closer, there's strong correlation between the PDE and its solution, which underpins our sorting algorithm.
2. Thermal Problem
We consider a two-dimensional thermal steady state equation, which can described by the fol- lowing equation (Sharma et al. 2018 Koric & Abueidda 2023):
8TT
PT
0,
Ox?
0y²
where T is the temperature. We examine the steady-state thermal equation in thermodynamics. The temperatures on the left and right boundaries are determined by random values ranging from -100 to and 0 to 100, respectively. These temperature values fundamental to sort approach.
With Elemert Mesh Displaydd
Tempersiur,
0.6
0.6
02-
02
-24
0.4
88
45
488
44
02
56
0.2
Figure 6: Finite element mesh and solution for the thermal problem
18

---

Published as a conference paper at ICLR 2024
As shown in the Figure 6 for the 'Thermal Problem', chose irregular boundary test the performance of our algorithm. The left image displays the finite element mesh obtained using FEM, an essential Step|Ajin converting the PDE into system of linear equations, previously mentioned. The right image showcases the solution one of the PDEs.
3. Poisson Equation
We consider a wwo-dimensional Poisson equation, which can be described by the following equa-
tion (Hsieh et al. 2019 Zhang et al. 2022
u f.
Physical Contexts in which the Poisson Equation Appears: Electrostatics; 2. Gravitation; 3. Fluid Dynamics.
We address the Poisson equation within square domain. The boundary conditions all four sides, as well as the f value on the left side of the equation, are generated using truncated Chebyshev poly- nomials. The coefficients of these five Chebyshev polynomials the basis for sorting (Driscoll et al. 2014).
Pocssnn Equatioo Solution
Poeraon Ecustion So uilon
0.
08
06
400
0.4
45
02
-16
Figure 7: Solution two PDEs before sorting
Equation
0.1
05
05
4.6
0.4
02
02
Figure 8: Solution of two PDEs after sorting
As depicted in the figure, with the matrix two-norm serving as the metric for distance, we provide visual representations of the Poisson Equation solutions featured in this study. The first Figure displays the solutions of the two adjacent equations before sorting, while the second Figure shows the solutions of the adjacent equations after sorting. It is evident that the sorting algorithm has enhanced the correlation between the preceding and following PDEs.
4. Helmholtz Equation We consider two-dimensional Helmholtz equation, which can be described
by the following equation (Zhang et al. 2022):
2 u 12 0,
Physical Contexts in which the Helmholtz Equation Appears: 1. Acoustics; Electromagnetism; 3. Quantum Mechanics.
19

---

Published as a conference paper at ICLR 2024
In Helmholtz Equation, k is the wavenumber, related to the frequency of the and the properties of the medium in which the wave is propagating. In our experiment, derived using the GRF method. The parameters inherent to the GRF serve as the foundation for our sort scheme.
22
22
30
39
Figure Solutions of Helmholtz equations with close parameters
Figure 10: Solutions of Helmholtz equations with divergent parameters
As depicted the Figure 10 with the matrix two-norm serving as metric for distance, comparable pattern emerges with the Helmholtz equations. Specifically, when the parameters of the neural operator are closely matched Helmholtz equations, a pronounced correlation exists between the PDE and its solution. In contrast, this correlation diminishes for equations with significantly differing parameters.
D.3 PRECONDITION
The experiments in this paper involve the following preconditioning techniques.
1. None: When preconditioning is set to 'none', essentially, no preconditioning operation is per- formed. This implies that the linear system is iteratively solved as-is, without any prior transforma- tions or modifications.
Diagonal Preconditioning (Jacobi) (Saad) 2003): This is simple preconditioning approach that considers only the diagonal elements of the coefficient matrix. The preconditioner matrix is the inverse of the diagonal elements of the coefficient matrix.
3. Block Jacobi (BJacobi) (Benzi et al. 1996): An extension of the Jacobi preconditioning, where the coefficient matrix broken down into smaller blocks, each corresponding to a subdomain or subproblem. Each block is preconditioned independently using its diagonal part.
Successive Over-relaxation (SOR) Youngg 1954): The SOR method a variant of the Gauss- Seidel iteration, introducing relaxation factor accelerate convergence. The preconditioner trans- forms the original problem into a weighted new problem, typically aiding in accelerating conver- gence for some iterative methods.
Additive Schwarz Method (ASM) (Toselli & Widlund 2004): ASM a domain decomposition method, where the original problem is split into multiple subdomains subproblems. Problems of
20

---

Published as a conference paper at ICLR 2024
each subdomain are solved independently. and these local solutions then combined into global solution.
6. Incomplete Cholesky (ICC) (Lin & Moré 1999 ICC is preconditioning method based on the Cholesky dccomposition, but drops certain off diagonal elements during the decomposition, making it "incomplete". It's utilized for symmetric positive definite problems.
ILU is based on LU decomposition, but like ICC, drops certain off diagonal elements during the decomposition. ILU can be applied to nonsymmetric prob- lems.
Incomplete LU (ILU) [Sad 2003):
D.4 ENVIRONMENT
To ensure consistency our evaluations, all comparative experiments were conducted under uni- form computing environments. Specifically, the environments used detailed follows:
1.
Environment (Envl):
Platform: Docker version 20.10.0
Operating System: Ubuntu 22.04.3 LTS
Processor: Dual-socket Intel Xeon Gold 6154 CPU, clocked 3.00GHz
Environment (Env2):
Platform: Windows 11, version 21H2, WSL
Operating System: Ubuntu 22.04.3 LTS
Processor: 13th Gen Intel Core 17-13700KF, clocked at 3.40 GHz
D.5
ANALYSIS OF RELEVANT EXPERIMENTAL RESULTS
D.5.1
CONVERGENCE SPEED ANALYSIS (TIME)
1e-1
1e-2
tolerance (2 norm)
Solver
Precondition
Slope
1e-3
SKR(ours) None
SKR(oures)
None
-6.80 x 10
SKR(ours) Jacobi
(KR(ee)s
Jacohi
-1.25 x 10
1e-4
SKR(ours) Bjacobi
(KR(e
(KR(eers
BJacobi
SOR
-2.22 10
-2.13 x
SKR(ours) SOR
SKR
ASM
-2.85 x 10
1e-5
SKR(ours) ASM
GMRES
None
-6.27 x 10
GMRES None
GMRES
Jacobi
-4.17 x 10
GMRES Jacobl
GMRES
Blacobi
-2.27 x 10-b
1e-6
GMRES Bjacobi
GMRES
10-3
GMRES SOR
GMRES
ASM
-1.38 10
1e-7
GMRES ASM
2
4
time (s)
Figure 11: Illustrating the Helmholtz Equation with matrix size of 10t, the graphic delves into the relationship between computational accuracy and average computation time. The left plot displays convergence curves under various preconditions, with the x-axis as average computation time and the y-axis as computational accuracy. The right plot presents linear fits for the three points with the minimum tolerance for each precondition, providing slopes as metrics for their convergence trends.
It is evident that the computational speed of our SKR algorithm far exceeds that of the GMRES algorithm. GMRES and similar Krylov subspace convergence algorithms exhibit a phenomenon known as superlinear convergence phenomenon in Appendix C. which means in log-error versus time or iterations graph, the GMRES convergence curve can be interpreted as combination of two segments (the first with smaller absolute slope value. followed by larger one). Due to the presence of superlinear convergence, an equitable and objective assessment of our SKR and GMRES in the high-precisiio phase requires discarding the data from the low-precision phase. Therefore,
21

---

Published as a conference paper at ICLR 2024
selected three high-pue righ-precision points for each algorithm to perform a linear fit and obtain the convergence slopes for this phase. A comparative analysis revealed that our algorithm's s slope has greater absolute value, indicating that our SKR also converges more quickly than GMRES in the high-precision stage.
D.5.2
CONVERGENCE SPEED ANALYSIS (ITERATION)
1e-1
1e-2
Solver
Precondition
Slope
tolerance (2 norm)
1e-3
SKR(oures]
None
-7.27 x 10
(KR(eurss)
Jacobi
-1.44 x 10
SKR(eurs)
SKR (oRree
BJacohi
-2.68
1e-4
SKR(ours) Jacabi
SKR(oures
SOR
-2.60 10
SKR(ours) Bjocobi
SKR(ours) SOR
(KR(ourss)
ASM
-4.49 x 10
1e-5
SKR(ours) ASM
GMRES
None
-3.31 x 10
GMRES None
GMRES
Jacobi
-2.37 x 10
1e-6
GMRES Jacobi
GMRES
BJacobi
-1.18
GMRES Bjacabi
GMRES
SOR
10
GMRES SOR
GMRES
ASM
-1.29 x 10
1e-7
GMRES ASM
2000
4000
6000
8000
10000
Average iterations
Figure 12: Illustrating the Helmholtz Equation with matrix size of 10 the graphic delves into the relationship between computational accuracy and average iteration count. The left plot displays convergence curves under various preconditions, with the x-axis as average iteration count and the y-axis as computational accuracy. The right plot presents linear fits for the three points with the minimum tolerance for each precondition, providing slopes as metrics for their convergence trends.
It is evident that, at the same level of precision, SKR algorithm requires significantly fewer iter- ations than GMRES. This strongly supports the notion that SKR achieves acceleration by expediting the convergence of the Krylov subspace, thereby reducing the number of iterations. As mentioned in Appendix D to appraise the convergence rates of both algorithms more equitably during the high- precision phase, we conducted linear fits using three high- precision data points. The compar- ison reveals that the absolute value of the slope during the high-precision phase is greater for our SKR algorithm than for GMRES, indicating faster rate of convergence. This robustly validates the superior stability of the SKR algorithm compared to GMRES.
D.5.3 STABILITY ANALYSIS
we can draw the following conclusions: 1. The use of the SKR algorithm almost never reaches the maximum iteration count, indicating its excellent iterative stability. 2. GMRES, on the other hand, frequently reaches the maximum iteration count without converging, and the effectiveness different preconditions varies significantly.
Based on the Figure 13
Combining these Results D.5.2 with the previous experiment, it is evident that our SKR algorithm exhibits much greater stability compared to the GMRES algorithm.
D.6 DETAILED EXPERIMENTAL DATA
We conducted nearly 3,000 experiments. Each experiment employed dataset specifically designed to mimic genuine NO training data. Below are some actual data from these experiments.
The title of each table below specifies the corresponding experimental dataset, preprocessing method, experimental environment, and details MPI parallelization. For instance, 'Table Darcy Flow, None, Env1, MPI72' indicates that the experimental results for the Darcy Flow problem recorded, with no matrix preprocessing algorithm applied, all within computing environment Env1, and all experiments utilizing MPI with 72 parallel threads.
As another example. 'Table 17] Poisson Equation, None, En 7153: MPI10, 11237: MPI10, 20245: MP120, 45337: MPI5, 71313: MP15 details that the table records the experimental results
22

---

Published as a conference paper at ICLR 2024
SKR(ours) SOR
SKR(ours) None
SKR(ours) Bjacobi
80
SKR(ours) ASM
SKR(ours) Jacobi
GMRES SOR
GMRES None
GMRES Bjacobi
GMRES ASM
Max Count
GMRES Jacobi
40
20
10-7
10-6
10-5
10-4
10-3
10
10-1
Tol
Figure 13: Illustrating the proportion of instances where different algorithms reach the maximum iteration count, the graphic delves into performance under various precisions for the Darcy flow problem with matrix size of 10 and maximum iteration count of 10
for the Poisson Equation, with no matrix preprocessing selected and all within the computing envi- ronment Envl. Here, experiments with matrix sizes of 7153 and 11237 were run with MPI using 10 parallel threads, size 20245 with MPI using 20 threads, and sizes 45337 and 71313 with MPI using threads.
These tables are divided into two parts:
the upper half records the average time to solve each linear system for two different algorithms under various conditions, in seconds, with smaller numbers indicating faster computation.
The lower half records the average number of iterations needed to solve each linear system for the two algorithms under different conditions, with fewer iterations indicating better algorithm stability and faster convergence.
The first row indicates different required error precisions, and the second column shows the matrix sizes derived from different grid densities. Each row's data composed two parts, with upper half showing the data for the GMRES algorithm-our baseline- and the lower half showing the corresponding data for our SKR algorithm.
23

---

Published as a conference paper at ICLR 2024
Table 3: Darcy Flow, None, Envl, MP172
solver)
1.E-01
1.F-02
1.E-03
1.E004
1.F-05
1.E-06
1.F-07
1.F-08
GMRES
0.13
0.20
0.24
0.27
0.29
0.31
0.34
0.36
2500
SKR
0.08
0.09
0.10
0.11
0.11
0.12
0.13
0.13
GMRES
0.28
0.45
0.55
0.61
0.66
0.69
0.71
0.72
6400
SKR
0.14
0.17
0.19
0.21
0.23
0.24
0.25
0.27
GMRES
0.39
0.70
0.79
0.85
0.88
0.89
0.91
10000
SKR
0.21
0.28
0.32
0.36
0.41
0.45
0.53
0.52
time
GMRES
3.89
5.62
6.00
6.10
6.19
6.28
22500
SKR
1.10
2.06
2.48
2.97
3.32
3.62
GMRES
26.28
51,20
57.63
60.89
62.88
64.83
40000
SKR
15.19
35.80
32.60
41.01
44.22
42.63
GMRES
1,394
2,575
3,115
3,608
3,983
4,325
4641
4963
2500
SKR
96
131
151
167
183
199
213
227
GMRES
3,192
5,402
6,709
7,617
8,193
8,582
8847
9059
6400
SKR
218
281
321
353
387
415
443
GMRES
3,874
7,457
8,426
9,022
9,375
9,596
9713
9775
10000
SKR
313
446
520
588
682
758
898
881
iter
GMRES
6,225
8,971
9,614
9,818
9,891
9,956
22500
SKR
745
1,445
1,763
2,122
2,367
2,600
GMRES
38,896
76,287
85,637
90,043
93,298
95,869
40000
SKR
9,283
21,950
19,994
25,212
27,207
26,177
Table 4: Darcy Flow, Jacobi, Envl, MP172
solver) tol
I.E-01
1.E-02
1.E-03
I.E-04
1.E-05
1.E-06
1.E-07
1.E-08
GMRES
0.13
0.19
0.22
0.24
0.26
0.28
0.30
0.31
2500
SKR
0.07
0.09
0.09
0.10
0.10
0.11
0.11
0.12
GMRES
0.25
0.44
0.52
0.58
0.62
0.65
0.68
0.70
6400
SKR
0.12
0.15
0.17
0.18
0.19
0.20
0.22
0.23
GMRES
0.64
0.76
0.82
0.85
0.88
0.91
10000
SKR
0.18
0.23
0.26
0.29
0.31
0.33
0.34
0.38
time
GMRES
0.70
1.24
1.33
1.38
1.40
1.40
22500
SKR
0.67
0.92
1.26
1.58
1.75
1.95
GMRES
23.72
47.98
53,99
57.54
60.86
62.79
40000
SKR
17.09
25.65
25.81
28.98
29.86
29.45
GMRES
1370
2357
2786
3110
3407
3695
3953
4192
2500
SKR
81
110
126
140
152
164
176
187
GMRES
2741
5235
6271
7074
7621
8037
8355
8622
6400
SKR
173
231
261
286
311
332
356
376
GMRES
3370
6642
7949
8611
8987
9322
9503
9610
10000
SKR
252
347
400
440
477
522
$49
613
iter
GMRES
4770
8741
9417
9778
9900
9921
22500
SKR
812
1 1148
1606
2024
2256
2526
GMRES
35088
70796
80000
85202
89638
92399
40000
SKR
10523
15770
15762
17713
18395
18073
24

---

Published as a conference paper at ICLR 2024
Table 5: Darcy Flow, BJacobi, Envl, I, MPI72
solver)
1.F-01
1.F-02
1.F-03
1.E-04
1.E-05
1.F-06
1.F-07
1.F-08
GMRES
0.09
0.14
0.17
0.20
0.22
0.23
0.25
0.27
2500
SKR
0.07
0.08
0.08
0,08
0.09
0.09
0.10
0.10
GMRES
0.22
0.39
0.49
0.56
0.60
0.64
0.67
0.70
6400
SKR
0.10
0.12
0.13
0.14
0.15
0.16
0.16
0.17
GMRES
0.30
0.59
0.70
0.78
0.85
0.89
0.97
10000
SKR
0.13
0.16
0.18
0.19
0.21
0.22
0.23
0.24
time
GMRES
0.49
1.08
1.29
1.41
1.50
1.55
22500
SKR
0.32
0.35
0.55
0.53
0.50
0.54
GMRES
15.43
32.62
32.00
36.61
49.46
51.52
40000
SKR
4.70
9.14
8.35
11.76
12.55
14.58
GMRES
676
1289
1675
1934
2182
2411
2636
2814
2500
SKR
57
76
87
96
104
113
121
129
GMRES
3642
4779
5497
5934
6310
6614
6916
6400
SKR
117
154
175
191
207
222
237
251
GMRES
2380
5015
5926
6712
7266
7705
8030
8347
10000
SKR
145
205
233
256
277
297
317
339
iter
GMRES
2667
6276
7520
8352
8839
9135
22500
SKR
338
378
630
601
568
614
GMRES
21396
44977
44643
51229
68445
71642
40000
SKR
2739
5382
4980
7030
7450
8713
Table 6: Darcy Flow, SOR, Envl, MPI72
solver tol
I.E-01
1.E-02
1.E-03
I.E-04
1.E-05
1.E-06
1.E-07
1.E-08
GMRES
0.08
0.13
0.16
0.18
0.19
0.20
0.21
0.23
2500
SKR
0.06
0.08
0.08
0.08
0.09
0.09
0.09
0.10
GMRES
0.20
0.33
0.41
0.46
0.50
0.54
0.57
0.60
6400
SKR
0.10
0.12
0.13
0.14
0.15
0.15
0.16
0.17
GMRES
0.28
0.53
0.62
0.69
0.74
0.79
0.82
0.85
10000
SKR
0.13
0.16
0.18
0.19
0.20
0.22
0.22
0.24
time
GMRES
0.44
0.99
1.17
1.30
1.38
1.41
22500
SKR
0.25
0.36
0.41
0.54
0.51
0.63
GMRES
15.72
35.73
41.95
46.75
50.32
52.57
40000
SKR
8.63
8.68
11.27
16.21
11.08
12.65
GMRES
622
1393
1822
2084
2271
2483
2675
2852
2500
SKR
57
87
96
105
113
121
129
GMRES
2059
3644
4688
5315
5833
6238
6611
6956
6400
SKR
119
1.56
175
192
209
224
238
253
GMRES
2439
5120
6028
6718
7246
7778
8136
8467
10000
SKR
1.50
206
236
259
281
303
320
341
iter
GMRES
2602
6261
7536
8371
8882
9150
22500
SKR
257
391
460
621
592
739
GMRES
22214
50690
59841
66630
71761
74751
40000
SKR
5146
5136
6735
9743
7635
25

---

Published as a conference paper at ICLR 2024
Table 7: Darcy Flow, ASM, Envl, MP172
solver)
1.E-01
1.F-02
1.F-03
1.E-04
1.E-05
1.E-06
1.F-07
1.F-08
GMRES
0.05
0.07
0.08
0.09
0.10
0.11
0.12
0.13
2500
SKR
0.06
0.06
0.07
0.07
0.07
0.08
0.08
0.08
GMRES
0.14
0.21
0.30
0.37
0.42
0.48
0.52
0.57
6400
SKR
0.09
0.10
0.10
0.11
0.12
0.12
0.12
0.13
GMRES
0.29
0.56
0.67
0.75
0.82
0.89
1.01
10000
SKR
0.11
0.13
0.14
0.15
0.16
0.16
0.17
0.18
time
GMRES
0.51
1.14
1.46
1.69
1.84
1.95
22500
SKR
0.21
0.29
0.33
0.44
0.39
0.42
GMRES
19.54
42.58
50,62
57.46
61.99
66.30
40000
SKR
2.55
7.57
3.64
9.26
6.42
9.41
GMRES
100
233
320
408
500
570
638
726
2500
SKR
32
42
48
53
58
62
67
71
GMRES
679
1199
1761
2215
3003
3316
3631
6400
SKR
99
109
117
126
133
141
GMRES
1502
3144
3783
4284
4681
5101
5461
5820
10000
SKR
92
125
141
153
166
178
189
201
iter
GMRES
1982
4748
6098
7088
7763
8236
22500
SKR
177
265
306
441
378
414
GMRES
18042
39800
47228
53725
57977
61608
40000
SKR
3664
1735
4526
3112
4611
Table 8: Darcy Flow, ICC, Env2, No MPI
solver tol
I.E-01
1.E-02
1.E-03
1.E-04
1.E-05
1.E-06
1.E-07
1.E-08
GMRES
0.01
0.03
0.03
0.04
0.04
0.04
0.05
0.05
2500
SKR
0.02
0.03
0.04
0.04
0.04
0.05
0.05
0.05
GMRES
0.11
0.29
0.38
0.47
0.60
0.71
0.81
0.91
6400
SKR
0.11
0.15
0.16
0.18
0.19
0.20
0.21
0.22
GMRES
0.63
1.27
1,65
1.97
2.22
2.43
2.85
10000
SKR
0.24
0.31
0.35
0.37
0.40
0.43
0.45
0.48
time
GMRES
2.15
3,49
4.46
5.43
6.17
6.76
22500
SKR
0.56
0.77
0.87
0.95
1.03
1.10
GMRES
36.38
53.80
58.95
63.86
69.77
75.75
40000
SKR
1.72
2,70
3.07
14.62
3.65
3.97
GMRES
32
149
172
196
221
245
270
298
2500
SKR
25
32
37
41
44
48
54
GMRES
236
629
825
1032
1320
1551
1787
2000
6400
SKR
49
66
74
82
88
94
100
106
GMRES
878
1769
2311
2752
3117
3405
3684
3985
10000
SKR
71
96
108
117
126
135
143
151
iter
GMRES
1208
1964
2506
3051
3460
3797
22500
SKR
72
103
118
130
141
150
GMRES
11479
16958
18618
20189
22018
23811
40000
SKR
130
210
240
1184
288
313
26

---

Published as a conference paper at ICLR 2024
Table 9: Darcy Flow, ILU, Env2, No MPI
solver tol
1.-01
1.E-02
1.E-03
1.E-04
1.E-05
1.E-06
1.E-07
1.E-08
GMRES
0.01
0.03
0.03
0.03
0.04
0.04
0.04
0.05
2.000
SKR
0.02
0.03
0.03
0.04
0.04
0.04
0.05
0.05
GMRES
0.10
0.27
0.37
0.47
0.56
0.64
0.72
0.80
6400
SKR
0.11
0.16
0.17
0.18
0.19
0.20
0.22
GMRES
0.63
1.19
1.52
1.88
2.13
2.36
2.55
2.75
10000
SKR
0.23
0.30
0.34
0.37
0.39
0.42
0.44
0.47
time
GMRES
2.38
3.96
4.86
5.90
6.69
7.42
22500
SKR
0.66
0.89
1.01
1.10
1.19
1.28
GMRES
30.39
53.52
61.25
68.44
72.64
76.58
40000
SKR
1.71
2.65
2.98
3.35
3.59
3.82
GMRES
30
141
157
179
199
219
242
261
2500
SKR
24
30
34
38
42
45
51
GMRES
222
599
823
1034
1237
1409
1594
1760
6400
SKR
49
63
72
79
85
92
97
103
GMRES
889
1670
2133
2653
2995
3328
3587
3836
10000
SKR
71
94
105
115
123
132
140
148
iter
GMRES
1414
2361
2896
3521
3997
4419
22500
SKR
89
122
140
154
167
180
GMRES
9768
17195
19676
21936
23358
24658
40000
SKR
127
204
231
256
277
300
Table 10: Thermal Problem, None, Envl, No MPI
solveri
IE-04
I.E-06
I.E-07
1.E-08
I.E-09
I.B-11
GMRES
0.21
0.29
0.37
0.45
0.54
0.62
0.70
2755
SKR
0.11
0.14
0.16
0.19
0.22
0.23
0.27
GMRES
1.74
2.53
3.32
4.12
4.93
5.73
6.53
7821
SKR
0.44
0.58
0.71
0.85
1.01
1.13
1.26
GMRES
3.44
5.08
6.71
8.34
9.99
11.60
13.21
11063
SKR
0.76
1.01
1.30
1.56
1.83
2.05
2.30
time
GMRES
7.48
11.29
15.11
18.92
22.57
26.04
26.50
17593
311.57
SKR
5.12
6.97
GMRES
27.23
85.47
105.81
106.87
70031
SKR
5.54
13.85
20.37
25.65
GMRES
503
693
882
1075
1278
1476
1674
2755
SKR
50
66
76
91
106
118
133
GMRES
1468
2140
2819
3500
4181
4864
5547
7821
SKR
104
130
162
216
243
GMRES
2057
4015
4996
5977
6946
7908
11063
iter
98.30
17593
SKR
130
174
217
278
326
375
439
GMRES
1909
4592
7295
9816
31157
SKR
85
178
249
345
GMRES
2531
7967
9888
10000
70031
SKR
115
300
442
S58
27

---

Published as a conference paper at ICLR 2024
Table 11: Thermal Problem, Jacobi, Envl, No MPI
solver'(td
1.E0
LE-00
1.E-U6
L.E-U7
LE-08
1.E-09
LE:10
LE-11
GMRES
0.20
0.26
0.31
0.37
0.43
0.50
0.58
2755
SKR
0.11
0.15
0.16
0.19
0.20
0.23
0.26
GMRES
1.37
1,78
2.14
2.46
2.73
3.04
3.36
SKR
0.44
0.56
0.71
0.89
1.01
1.13
1.24
GMRES
2.49
3.35
4.12
4.80
5.42
6.04
6.69
11063
SKR
0.75
1.01
1,30
1.57
1.77
2.02
2.28
time
GMRES
17593
17.77
22.71
31157
SKR
1.70
3.77
5.18
6.99
GMRES
15.67
34.74
68. 17
101.00
70031
SKR
5.88
13.55
19.76
24.83
GMRFS
455
606
732
857
994
1169
1347
SKR
101
114
129
GMRES
1142
1490
1788
2049
2283
2530
2801
7821
SKR
104
131
164
214
GMRES
11063
iter
GMRES
1334
2111
2941
3552
4090
4646
5466
17593
SKR
132
172
214
262
321
379
428
GMRES
1186
2078
3713
4737
31157
SKR
77
184
251
343
GMRES
1450
3216
5767
8316
70031
SKR
122
291
431
541
Table 12: Thermal Problem, BJacobi, Env1, No MPI
solver'<eo [tol
LE-04
1.E-05
1.E-06
1.E-07
LE-08
1.E-09
LE-10
LE-11
GMRES
0.07
0.09
0.11
0.13
0.15
0.17
0.19
2755
SKR
0.05
0.08
0.08
0.09
0.10
0.12
0.12
GMRES
0.45
0.58
0.72
0.84
0.97
1.10
1.23
SKR
0.23
0.31
0.33
0.36
0.43
0.46
0.53
GMRES
0.85
1.13
1.41
1.69
1.95
2.21
2.46
SKR
0.36
0.47
0.51
0.62
0.66
0.77
0.81
time
GMRES
1.55
2.10
2.65
3.18
3.72
4.27
4.81
17593
1.23
1.50
31157
SKR
1.24
1.49
1.99
2.47
GMRES
13.00
21.62
27.71
34.00
70033
SKR
2.97
5.09
7.03
8.86
GMRES
114
154
188
228
266
306
344
SKR
21
35
40
-55
.49
GMRES
289
376
466
549
633
718
802
7821
S SKR
74
83
GMRES
1142
11063
iter
GMRES
454
616
778
936
1094
1256
1418
17593
107
118
132
GMRES
522
902
1206
1498
31157
SKR
46
62
83
103
GMRES
832
1383
1793
2132
SKR
50
95
136
173
28

---

Published as a conference paper at ICLR 2024
Table 13: Thermal Problem, SOR, Env1, No MPI
solver [tol
LEE044
1.E-U5
1.E-06
1.E-07
LE-U8
LE-00
L.E-10
LE-11
GMRES
0.06
0.09
0.10
0.13
0.14
0.17
0.19
2755
SKR
0.07
0.08
0.09
0.10
0.12
0.12
0.13
GMRES
0.58
0.79
1.00
1.20
1.41
1.61
1.79
SKR
0.25
0.33
0.36
0.44
0.47
0.55
0.58
GMRES
0.86
1,16
1.45
1.77
2,08
2.41
2.75
11063
SKR
0.44
0.51
0.62
0.75
0.80
0.92
1.04
time
GMRES
4.78
7.79
8.70
17893
1.74
6.96
9.89
12.66
31157
SKR
1.16
1.82
2.36
2.99
GMRES
12.99
24.30
31.12
37.73
70031
SKR
3.91
6.36
7.21
8.83
GMRES
106
1.1
183
228
255
304
335
2755
SKR
50
55
GMRES
371
505
639
775
907
1035
1158
7821
SKR
GMRES
527
1256
11063
iter
GMRES
751
1064
1370
1670
1962
2238
2500
17893
SKR
107
122
138
154
GMRES
619
1123
1594
2036
31157
SKR
45
77
99
129
GMRES
807
1600
2182
26.37
70031
SKR
49
113
145
200
Table 14: Thermal Problem, ASM, Env1, No MPI
solver'<ed too
LE-04
1.E-0
1.E-06
1.E-07
L.E-08
L.E009
LE-10
LE-11
GMRES
0.08
0.10
0.12
0.14
0.16
0.18
0.21
2755
SKR
0.05
0.08
0.09
0.09
0.10
0.12
0.13
GMRES
0.48
0.62
0.77
0.90
1.03
1.17
1.31
SKR
0.24
0.32
0.34
0.37
0.44
0.47
0.54
GMRES
0.92
1.21
1.51
1.81
2.08
2.37
2.61
SKR
0.37
0.48
0.53
0.64
0.68
0.79
0.83
time
GMRES
1.66
2.24
2.82
3.38
3.96
4.54
5.12
17593
1.26
31157
SKR
1.28
1.53
2.05
2.54
GMRES
12.67
21.09
29.90
35.49
70033
SKR
3.06
5.23
7.20
9.06
GMRES
114
154
188
228
266
306
344
SKR
21
35
40
455
.49
GMRES
289
376
466
549
633
718
802
7821
SKR
74
83
GMRES
1142
11063
iter
GMRES
454
616
778
936
1094
1256
1418
17593
107
118
132
GMRES
522
902
1206
1498
31157
SKR
46
83
103
GMRES
832
1383
1793
2132
SKR
50
95
136
173
29

---

Published as a conference paper at ICLR 2024
Table 15: Thermal Problem, ICC, Env1, No MPI
solver [tol
LEE044
1.E-U5
1.E-06
1.E-07
LE-U8
LE-00
LE--10
LE-11
GMRES
0.07
0.09
0.11
0.13
0.15
0.17
0.19
2755
SKR
0.05
0.08
0.09
0.09
0.30
0.12
0.12
GMRES
0.46
0.59
0.73
0.86
0.99
1.12
1.25
SKR
0.24
0.31
0.33
0.36
0.43
0.46
0.53
GMRES
0.87
1.16
1.44
1.72
2,00
2.26
2.50
11063
SKR
0.36
0.47
0.52
0.62
0.66
0.77
0.81
time
GMRES
3.26
4.93
17893
1.23
5.62
7.51
9.33
31157
SKR
1.24
1.51
2.02
2.50
GMRES
14.18
22.57
28.99
34.54
70031
SKR
3.00
3.10
7.04
8.92
GMRES
114
1.5
188
228
266
306
344
2755
SKR
40
49
GMRES
289
376
466
549
633
718
802
7821
SKR
GMRES
1142
11063
iter
GMRES
454
616
778
936
1094
1256
1418
17893
SKR
79
107
118
132
GMRES
522
902
1206
1498
31157
SKR
46
62
83
103
GMRES
832
1383
1793
2132
70031
SKR
50
95
136
173
Table 16: Thermal Problem, ILU, Env1, No MPI
solver'<eo tool
LE-04
1.E00
1.E-06
1.E-07
LE-08
1.E-09
LE-10
LE-11
GMRES
0.07
0.09
0.11
0.13
0.15
0.17
0.19
2755
SKR
0.05
0.08
0.08
0.09
0.10
0.12
0.12
GMRES
0.45
0.57
0.71
0.83
0.96
1.09
1.21
SKR
0.23
0.31
0.33
0.36
0.43
0.45
0.53
GMRES
0.85
1.12
1.40
1.67
1.94
2.18
2.43
SKR
0.36
0.47
0.51
0.62
0.66
0.76
0.81
time
GMRES
1.54
2.09
2.63
3.16
3,70
4.24
4.79
17593
1.50
31157
SKR
1.23
1.49
2.00
2.47
GMRES
12.99
21.46
29.25
33.26
70033
SKR
2.95
5.03
6.96
8.78
GMRES
114
1.4
188
228
266
306
344
SKR
21
35
40
-55
.49
GMRES
289
376
466
549
633
718
802
7821
SKR
74
83
GMRES
1142
11063
iter
GMRES
454
616
778
936
1094
1256
1418
17593
107
118
GMRES
522
902
1206
1498
31157
SKR
46
62
83
103
GMRES
832
1383
1793
2132
SKR
50
95
136
173
30

---

Published as a conference paper at ICLR 2024
Table 17: Possion Equation, None, Envl, 7153: MPI10, 11237:
MPI10, 20245: MPI20, 45337:
MPI5, 71313: MPI5
n
solver)tol
1.E-05
1.E-06
1.E-07
1.E-08
1.E-09
1.E-10
1.E-11
GMRES
0.07
0.09
0.11
0.13
0.15
0.17
0.19
7153
SKR
0,07
0.11
0.10
0.11
0.12
0.13
0.15
GMRES
0.24
0.33
0.41
0.50
0.66
0.75
11237
SKR
0.16
0.19
0.22
0.26
0.29
0.34
0.35
GMRES
0.41
0.57
0.74
0.90
1.06
1.24
1.39
20245
SKR
0.24
0.29
0.33
0.39
0.44
0.50
0.56
time
GMRES
1.81
2.74
3.69
4.64
5.60
6.56
7.51
45337
SKR
1,46
1.91
2.36
2.79
3.25
3.69
4.15
GMRES
3.89
6.13
8.39
10.74
13.02
15.33
17.68
71313
SKR
3.05
4.10
5.12
6.16
7.20
8.25
9.28
GMRES
274
374
475
576
676
776
878
7153
SKR
122
139
158
GMRES
362
506
651
796
942
1088
1234
11237
SKR
89
112
136
160
185
210
235
GMRES
554
804
1054
1306
1558
1810
2062
20245
SKR
128
165
201
239
278
316
355
iter
GMRES
958
1475
1998
2526
3054
3583
4112
45337
SKR
209
279
349
419
489
560
631
GMRES
1359
2158
3810
4641
5472
6303
71313
SKR
289
397
500
606
710
816
922
Table 18: Possion Equation, Jacobi, Envl. 7153: MPI10, 11237: MPI10, 20245: MPI20, 45337:
MPI5, 71313: MPI5
n
solver)
I.E-05
I.E-06
I_E-07
1.E-08
1.E-09
1.E-10
L.E-1
GMRES
0.07
0.09
0.12
0.13
0.16
0.17
0.19
7153
SKR
0.08
0.09
0.10
0.12
0.12
0.14
0.15
GMRES
0.25
0.33
0.42
0.50
0.68
0.78
11237
SKR
0.18
0.21
0.24
0.27
0.31
0.34
0.37
GMRES
0.41
0.58
0.75
0.91
1.06
1.23
1.40
20245
SKR
0.25
0.30
0.35
0.40
0.46
0.51
0.56
time
GMRES
1.85
2.78
3.73
4.69
5.65
6.61
7.69
45337
SKR
1.48
1.91
2.37
2.80
3.25
3.69
4.13
GMRES
3.94
6.16
8.41
10.69
12,92
15.18
17.41
71313
SSKR
3.07
4.08
5.09
8.10
9.12
GMRES
279
482
581
782
71.53
SKR
88
124
142
160
178
GMRES
368
512
658
802
947
1092
1236
11237
SKR
91
114
138
163
188
212
237
GMRES
558
804
1051
1298
1546
1794
2042
20245
SKR
129
165
201
238
276
314
352
iter
GMRES
970
1484
2003
2526
3049
3571
4094
45337
SKR
349
418
486
GMRES
1362
2148
2954
3762
4567
5369
6167
71313
SKR
291
394
495
596
699
800
904
31

---

Published as a conference paper at ICLR 2024
Table 19: Possion Equation, BJacobi, Envl, 7153: MPI10, 11237: MPI10, 20245: MPI20, 45337:
MPI5, 71313: MPI5
solver)tol
I.E-05
1.E-06
1.E-07
1.E-08
1.E-09
1.E-10
1.E-11
GMRES
0.07
0.08
0.10
0.11
0.13
0.14
0.15
7153
SKR
0,07
0.08
0.09
0.10
0.11
0.12
0.13
GMRES
0.20
0.26
0.31
0.36
0.46
0.51
11237
SKR
0.14
0.17
0.20
0.22
0.25
0.28
0.30
GMRES
0.34
0.47
0.58
0.71
0.84
0.94
1.06
20245
SKR
0.21
0.26
0.29
0.33
0.38
0.42
0.46
time
GMRES
1.29
1.77
2.22
2.63
3.02
3.40
3.76
45337
SKR
1.15
1,45
1.78
2.09
2.39
2.70
3.01
GMRES
2.57
3.53
4.34
5.09
5.86
6.62
7.45
71313
SKR
2.26
2.97
3.63
4.27
5.56
6.24
GMRES
198
262
327
393
458
585
7153
SKR
71
85
99
114
129
GMRES
262
349
432
513
591
668
743
11237
SKR
70
88
106
125
144
163
182
GMRES
419
595
768
940
1111
1279
1446
20245
SKR
102
132
160
191
220
251
280
iter
GMRES
588
823
1041
1241
1434
1614
1795
45337
SKR
152
199
246
294
340
386
434
GMRES
1081
1337
1577
1818
2061
2320
71313
SKR
200
269
334
396
462
523
588
Table 20: Possion Equation, SOR, Envl, 7153: MPI10, 11237: MPI5, 71313: MPI5
MPI10, 20245:
MPI20, 45337:
n
solver)
I.E-05
I.E-06
I.E-07
1.E-08
1.E-09
1.E-10
L.E-1
GMRES
0.07
0.08
0.09
0.11
0.12
0.14
0.15
7153
S SKR
0.07
0.08
0.09
0.10
0.11
0.12
0.13
GMRES
0.20
0.26
0.32
0.37
0.47
0.51
11237
SKR
0.15
0.18
0.21
0.23
0.26
0.28
0.31
GMRES
0.33
0.46
0.58
0.70
0.82
0.94
1.08
20245
SKR
0.22
0.26
0.32
0.35
0.39
0.43
0.47
time
GMRES
1.43
2.04
2.62
3.17
3.68
4.15
4.60
45337
SKR
1.16
1.47
1.80
2.11
2.42
2.73
3.06
GMRES
2,76
3.94
4.98
5.90
6.78
7.64
8.49
71313
SKR
2.31
3.03
3.69
4.37
5.72
6.37
GMRES
203
339
408
477
545
71.53
SKR
72
86
115
130
144
GMRES
270
362
450
534
614
692
766
11237
SKR
71
88
107
126
145
164
183
GMRES
421
599
775
953
1129
13055
1480
20245
SKR
102
132
160
191
220
251
280
iter
GMRES
647
942
1221
1484
1729
1960
2174
45337
SKR
298
344
390
439
GMRES
822
1192
1516
1804
2076
2342
2610
71313
SKR
203
272
339
405
470
535
599
32

---

Published as a conference paper at ICLR 2024
Table 21: Possion Equation, ASM, Env1, 7153:
MPHO0 11237:
MPHO,, 20245: MPI20, 45337:
MPI5, 71313: MPI5
n
solver tol
I.E-05
1.E-06
1.E-07
1.E-08
1.E-09
1.E-10
1.E-1
GMRES
0.07
0.09
0.10
0.12
0.14
0.15
0.17
7153
SKR
0.08
0.10
0.10
0.13
0.14
GMRES
0.27
0.35
0.42
0.49
0.56
0.63
0.69
11237
SKR
0.17
0.21
0.23
0.27
0.30
0.33
0.36
GMRES
0.48
0.66
0.83
1.01
1.18
1.35
1.53
20245
SKR
0,25
0.31
0.36
0.41
0.46
0.52
0.57
time
GMRES
1,42
1,93
2.42
2.87
3.32
3.71
4,10
45337
SKR
1,20
1.51
1.85
2.16
2,47
2.79
3.12
GMRES
2.80
3.86
4.71
5.52
6.35
7.18
8.08
71313
SKR
2.35
3.09
3.76
4.42
5.75
6.44
GMRES
198
327
393
458
7153
SKR
57
71
85
99
114
129
144
GMRES
262
349
432
513
591
668
743
11237
SKR
70
88
106
125
144
163
182
GMRES
419
595
768
940
1111
1279
1446
20245
SKR
102
132
160
191
220
251
280
iter
GMRES
1041
1434
1614
1795
45337
SKR
152
246
294
340
434
GMRES
775
1081
1337
1577
1818
2061
2320
71313
SKR
200
269
334
396
462
523
588
Table 22: Possion Equation, ICC, Envl, No MPI
solver
1.E-05
1.E-06
I.E-07
1.E-08
1.E-09
1.E-10
1.E-11
GMRES
0.09
0.10
0.10
0.11
0.12
0.13
0.13
7153
SKR
0.12
0.13
0.14
0.13
0.13
0.14
0.14
GMRES
0.15
0.17
0.19
0.20
0.21
0.23
0.25
11237
SKR
0.21
0.20
0.21
0.21
0.22
0.38
0.39
GMRES
0.32
0.37
0.40
0.45
0.50
0.55
0.59
20245
SKR
0.37
0.38
0.39
0.69
0.72
0.74
0.77
time
GGMRES
0.91
1.12
1.26
1.41
1.56
1.71
1.90
45337
SKR
0.88
1.58
1.66
1.75
1.83
1.93
2.27
GMRES
2.01
2.19
2.45
2.91
3.35
3.66
3.88
71313
SKR
2,50
2.67
2.86
3.11
3.62
3.81
3.97
GMRES
22
30
33
36
42
7153
SKR
16
18
GMRES
26
36
40
45
52
59
11237
SKR
20
17
19
18
20
22
24
GMRES
35
42
50
59
67
74
80
20245
SKR
18
20
19
22
25
28
31
iter
GMRES
49
62
75
84
97
110
120
45337
SKR
19
23
27
31
35
39
43
GMRES
98
139
154
167
71313
SKR
39
44
33

---

Published as a conference paper at ICLR 2024
Table 23: Possion Equation. ILU, Envl, No MPI
solver
1.E-05
1.E-06
1.E-07
1.E-08
1.B-09
1.E-10
1.E-1
GMRES
0.26
0.26
0.27
0.28
0.28
0.29
0.29
7153
SKR
0.29
0.30
0.32
0.31
0.32
0.32
0.33
GMRES
0.41
0.43
0.44
0.45
0.46
0.49
0.49
11237
SKR
0,48
0,50
0.51
0.51
0.53
0.52
0,53
GMRES
0.84
0.85
0.87
0.91
0.93
0.97
0.99
20245
SKR
0.94
0.96
0.94
0.95
0.98
1.00
time
GMRES
2.02
2.18
2.26
2.36
2.70
2.87
45337
SKR
2.15
2.24
2.22
2.29
2.98
3.02
3.13
MRES
3,46
3.67
3.96
4.27
4.52
4.79
4.93
71313
SKR
3.54
3.63
4.68
4.79
4.94
5.11
5.23
GMRES
15
18
20
22
24
26
27
7153
SKR
10
11
12
13
14
15
16
GMRES
22
24
27
29
31
33
11237
SKR
14
16
19
GMRES
27
31
34
44
20245
SKR
18
15
17
18
20
20
iter
GMRES
32
39
46
53
62
69
76
45337
SKR
17
20
20
21
24
26
29
GMRES
40
49
60
71
78
86
95
71313
SKR
20
23
26
29
33
36
Table 24: Helmholtz. Equation, None, Env1, MP172
solver tol
1.E-01
1.E-02
1.E-03
1.E-04
1.E-05
1.E-06
1.E-07
GMRES
0.81
1.37
1.75
2.07
2.36
2.61
2.87
2500
SKR
0.17
0.21
0.22
0.23
0.25
0.27
0.29
GMRES
1.74
3.09
3.72
4.14
4.49
4.77
4.93
6400
SKR
0.29
0.36
0.41
0.45
0.48
0.51
0.54
GMRES
2.44
4.06
4.88
5.22
5.41
5.51
5.57
10000
SKR
0.40
0.52
0.61
0.67
0.71
0.80
0.86
time
GMRES
16.22
33.05
39.39
43.37
46.67
49.21
22500
SKR
3.60
6.73
9.31
6.13
9.04
8.54
GMRES
1,305
2,279
2,936
3,520
3,999
4,461
4882
2500
SKR
98
130
151
168
185
199
214
GMRES
3,036
5,436
6,604
7,357
7,978
8,421
8749
6400
SKR
212
277
319
354
388
414
441
GMRES
4,187
7,052
8,478
9,118
9,407
9,574
9713
10000
SKR
303
407
473
538
570
648
710
iter
GMRES
26,565
53,855
64,542
71,353
76,793
80,991
22500
SKR
2,647
5,014
6,937
4,557
6,743
6,394
34

---

Published as a conference paper at ICLR 2024
Table 25; Helmholtz Equation, Jacobi, Envl, MPI72
n
solver tol
1.E-01
1.E-02
1.E-03
1.E-04
1.E-05
1.E-06
1.E-07
GMRES
0.72
1.20
1.45
1.65
1.87
2.07
2.28
2500
SKR
0.16
0.19
0.21
0.23
0.24
0.25
0.27
GMRES
1.45
2.64
3,30
3,71
4.08
4.32
4.54
6400
SKR
0.27
0.33
0.37
0.39
0.42
0.45
0.47
GMRES
2.18
3.82
4.57
5.00
5.15
5.33
5,40
10000
SKR
0.37
0,45
0.51
0.56
0.60
0.64
0.68
time
GMRES
16.31
29.04
37.38
41.38
43.97
45.58
22500
SKR
4.92
4.24
6.01
6.46
5,49
5.66
GMRES
1,155
1,988
2,405
2,762
3,136
3,509
3849
2500
SKR
82
110
126
140
153
164
176
GMRES
2,507
4,597
5,791
6,562
7,220
7,637
8028
6400
SKR
176
228
261
289
315
335
356
GMRES
3,726
6,621
7,959
8,630
9,004
9,237
9430
10000
SKR
249
327
373
415
451
489
521
iter
GMRES
26,707
47,340
60,878
67, 67,205
71,418
74,351
22500
SKR
3,611
3,113
4,442
4,786
4,054
4,175
Table 26: Helmholtz Equation, BJacobi, Envl, MP172
n
solver' tol
1.E-01
1.E-02
1.E-03
1.E-04
1.E-05
1.E-06
1.E-07
GMRES
0.50
0.75
0.96
1.12
1.23
1.35
1.47
2500
SKR
0.12
0.14
0.15
0.16
0.18
0.18
0.19
GMRES
1.09
2.01
2.48
2.89
3.18
3.48
3.71
6400
SKR
0.19
0.23
0.26
0.28
0.29
0.31
0.33
GMRES
1.49
2.60
3.36
3.86
4.27
4.55
4.72
10000
SKR
0.24
0.29
0.32
0.35
0.37
0.39
0.42
time
GMRES
13.13
19.82
25.01
28.58
30.87
32.45
22500
SKR
2.32
3.51
0.71
0.77
0.85
3.41
GMRES
737
1,154
1,513
1,772
1,979
2,168
2352
2500
SKR
58
76
87
96
105
113
122
GMRES
1,766
3,332
4,181
4,853
5,359
5,821
6272
6400
SKR
116
154
175
193
208
223
238
GMRES
2,405
4,227
5,541
6,376
7,019
7,511
7883
10000
SKR
147
197
222
244
265
283
301
iter
GMRES
20,335
30,564
38,699
44,223
47,387
50,197
22500
SKR
1,639
2,505
2,427
35

---

Published as a conference paper at ICLR 2024
Table 27: Helmholtz Equation, SOR, Env I, MPI72
n
solver tol
1.E-01
1.E-22
1.E-03
1.E-04
1.E-05
1.E-06
1.E-07
GMRES
0.48
0.77
0.93
1.10
1.19
1.30
1.38
2500
SKR
0.14
0.16
0.17
0.18
0.19
0.20
0.21
GMRES
1.17
1.98
2.40
2.77
3.04
3.30
3.59
6400
SKR
0.21
0.25
0.27
0.29
0.31
0.33
0.34
GMRES
1.46
2.60
3.22
3.77
4.16
4.43
4.66
10000
SKR
0.25
0.31
0.34
0.37
0.40
0.42
0.44
time
GMRES
12.28
18.00
23.99
27.31
30.17
32.18
22500
SKR
2.18
0.65
0.75
1.26
3.76
2.93
GMRES
719
1,231
1,507
1,782
1,966
2,134
2290
2500
SKR
57
76
88
96
105
114
122
GMRES
1,987
3,414
4,215
4,823
5,339
5,783
6262
6400
SKR
118
156
176
194
210
225
239
GMRES
2,441
4,405
5,503
6,432
7,074
7,595
7976
10000
SKR
149
197
226
248
270
288
308
iter
GMRES
19,540
28,939
38,419
43,545
48,285
51,382
22500
SKR
1,545
2,712
2,094
Table 28: Helmholtz Equation, ASM, Env1, MP172
n
solver's tol
1.E-01
1.E-02
1.E-03
1.E-04
1.E-05
1.E-06
1.E-07
GMRES
0.12
0.34
0.42
0.49
0.55
0.63
0.69
2500
SKR
0.13
0.14
0.15
0.16
0.17
0.18
0.18
GMRES
1.01
1.61
1.95
2.23
2.59
2.84
3.06
6400
SKR
0.19
0.22
0.23
0.25
0.26
0.27
0.28
GMRES
1.55
2.83
3.44
3.98
4.35
4.72
5.07
10000
SKR
0.22
0.27
0.29
0.31
0.33
0.34
0.36
time
GMRES
15.36
25.21
32.21
35.62
38.32
41.75
22500
SKR
0.95
0.56
2.33
0.75
0.77
3.08
GMRES
64
306
384
459
531
601
668
2500
SKR
33
42
48
53
58
63
67
GMRES
1,054
1,713
2,118
2,412
2,803
3,081
3328
6400
SKR
68
89
100
110
118
126
134
GMRES
1,611
3,034
3,659
4,264
4,676
5,068
5448
10000
SKR
91
119
136
149
160
170
182
iter
GMRES
15,683
25,566
32,785
36,475
38,992
42,431
22500
SKR
498
1,310
1,738
36

---

Published as a conference paper at ICLR 2024
Table 29: Helmholtz Equation, ICC, Env2, NoMPI
solver' tol
1.E-01
1.E022
1.E-03
1.E-04
1.E-05
1.E-06
1.E-07
GMRES
0.01
0.02
0.03
0.04
0.05
0.06
0.07
2500
SKR
0.03
0.04
0.04
0.05
0.05
0.05
0.06
GMRES
0.22
0.71
0.93
1.15
1.29
1.43
1.55
6400
SKR
0.13
0.17
0.18
0.20
0.21
0.23
0.24
GMRES
0.65
1.58
2.12
2.50
2.80
3.10
3.32
10000
SKR
0.25
0.33
0.37
0.41
0.44
0.47
0.49
time
GMRES
2.03
4.77
6.08
6.93
7.66
8.34
22500
SKR
0.62
0.90
1.02
1.12
1.21
1.30
GMRES
48
97
160
232
294
362
424
2500
SKR
28
38
43
48
52
55
59
GMRES
492
1,578
2,048
2,499
2,866
3,160
3431
6400
SKR
57
78
88
96
103
109
116
GMRES
915
2,229
2,983
3,518
3,942
4,356
4665
10000
SKR
104
117
129
139
149
159
iter
GMRES
1,149
2,704
3,442
3,926
4,338
4,733
22500
SKR
81
122
140
155
168
181
Table 30: Helmholtz Equation, ILU, Env2, NoMPI
solver' tol
1.E-01
1.E-02
1.E-03
1.E-04
1.E-05
1.E-06
1.E-07
GMRES
0.00
0.01
0.01
0.01
0.01
0.01
0.01
2500
SKR
0.01
0.02
0.03
0.03
0.03
0.03
0.04
GMRES
0.20
0.64
0.85
1.02
1.15
1.27
1.37
6400
SKR
0.13
0.16
0.18
0.19
0.21
0.22
0.23
GMRES
0.62
1.29
1.75
2.30
2.63
2.88
3.11
10000
SKR
0.25
0.33
0.37
0.40
0.43
0.46
0.49
time
GMRES
2.94
5.60
6.74
7.73
8.62
9.28
22500
SKR
0.74
1.04
1.19
1.31
1.43
1.52
GMRES
20
28
32
36
40
46
56
2500
SKR
21
25
28
30
32
34
GMRES
437
1,420
1,875
2,257
2,548
2,822
3056
6400
SKR
57
77
86
93
100
107
114
GMRES
869
1,818
2,483
3,254
3,726
4,085
4401
10000
SKR
76
104
115
127
137
147
157
iter
GMRES
1,768
3,373
4,058
4,649
5,181
5,587
22500
SKR
101
147
169
187
203
219
37

---

Published as a conference paper at ICLR 2024
E ADDITIONAL EXPERIMENTS
E. OTHER APPLICATION SCENARIOS
Neural operators is not a prerequisite for the application of our SKR algorithm. Indeed, our initial use of neural operators revealed the time-consuming issue of dataset generation. After an extensive study of Krylov methods, we developed the SKR algorithm, which has proven to be highly success- ful. Neural operators PDE problems, and the time-intensive nature of dataset generation has been a sigmificant hurdle in the field, which is why our paper is set in this context. However, the SKR algorithm' 's utility is not limited to neural operator problems and can be applied to other areas, such as:
of the most renowned and efficient data-driven methods for solving
1. Data-Driven PDE Algorithms: There are various data-driven PDE solving or parameter optimiza- tion algorithms in both traditional heuristics and AI for PDE domains. These algorithms often require the frequent generation of large volumes of PDE data for datasets, a process that our SKR algorithm can accelerate. For instance, Greenfeld et al. (2019) Hsieh et al. (2019).
2. Specific Physical Problem Solving: In certain scenarios, modeling the physical problems at hand translates into solving multiple related linear systems. Our SKR algorithm is adept at efficiently handling such issues. For instance, this is evident the modeling of fatigue and fracture through finite element analysis (Gullerud & Dodds Jr 2001), which employs dynamic loading across numerous steps, resulting a substantial number of interrelated linear systems. Similarly, in the resolution process of lattice quantum field theory (Muroya et al., 2003), the end goal involves solving large collection of correlated linear systems.
E.2
PARALLELIZATION STRATEGIES
The parallelism of our SKR algorithm is indeed a critical component. Typically, matrix algorithms employ three parallelization strategies.
E.2.1 MPI-BASED PARALLELIZATION
Using MPI to parallelize certain computational processes within the matrix algorithm, such as the Krylov subspace iteration in GMRES, which often relies on specific matrix preconditioning meth- ods. Regarding this MPI-based parallel approach, our experiments have conducted numerous com- parative tests in Appendix D.6, and our SKR algorithm has consistently achieved excellent results.
E.2.2 DECOMPOSE COMPUTATIONAL TASKS
Utilizing MPI to decompose the computational task into several independent parts, each solved separately. As you rightly mentioned in your question, this is indeed excellent approach. We will illustrate the parallelism strategy of our SKR algorithm with a specific example, showcasing its sustained and notable acceleration impact.
Imagine have server CPU with 72 parallel threads. Our parallel SKR algorithm comprises the following parts:
Sorting: As with the original SKR, we order these data points using a suitable sorting algorithm into sequence where consecutive parameters are strongly correlated, meaning the distance between parameter matrices is minimal (we usually use the 1, 2, or infinity norms of matrices the metric in this Banach space). The choice of sorting algorithm varies with the size of the dataset:
Greedy Sorting: Simple implement and negligible in computational cost for smaller datasets, such as our 7200 data points. However, for larger datasets, say 10" data points, the computational load of the greedy algorithm increases substantially and cannot be paral- lelized across multiple CPUs.
FFT Dimension Reduction + Fractal Division + Greedy Sorting: For very large datasets, propose parallelizable method that first reduces dimensionality via FFT manage the high-dimensional coordinates, then applies a fractal division algorithm based on the Hilbert curve, and finally, performs greedy sorting within each divided section.
38

---

Published as a conference paper at ICLR 2024
2. Solving Linear Systems: Assuming the task hand is to generate dataset comprising 7200 data points, upon sequencing the data to establish a strong correlation, we subsequently parti- tion these 7200 data points into 72 batches, each containing 100 points, distributing the tasks solving the corresponding linear systems to different CPU threads. Each thread employs our SKR algorithm for solving, achieving parallelism in the process.
To demonstrate the feasibility of our parallel approach and its effective acceleration, we have designed the following experiments.
An experimental simulation conducted to generate dataset for the Helmholtz equation, comprising 7200 data points. Each corresponding matrix has a dimension of 10000, with Suc- cessive Over-Relaxation (SOR) used uniformly for matrix preconditioning. Both our SKR al- gorithm and the baseline GMRES algorithm were implemented to distribute the solving of 7200 tasks across 72 threads, with each thread solving 100 linear systems. Due the potential variance in computation completion times across threads, the reported computation times and iteration counts have been averaged for consistency.
The upper half of the table presents the average computation time taken by both algorithms to solve each linear system, measured in seconds. The lower half details the average number of iterations required to solve each linear system. The first row of the table specifies the precision requirements for solving linear systems.
Table 31: Comparative Experimets of Parallel Versions
1E-03
IE-05
1E-07
time(s)
Parallel GMRES
0.08
0.105
0.122
time(s)
Parallel SKR(ours)
0.011
0.013
0.015
iter
Parallel GMRES
3734
4906
5715
iter
Parallel SKR(ours)
126
148
167
The experimental results reveal that our SKR algorithm achieves 6.7-8.0 fold acceleration in computation time and 30-34 fold reduction in the number of iterations compared to the baseline.
E.2.3 BLOCK PARALLE
The block concept is a strategy for parallelizing large matrices and reducing memory overhead. Drawing from the idea of Krylov blocks, we have redesigned block version of our SKR algorithm to facilitate parallel processing. The fundamental principle involves transforming the conventional matrix algorithm into block matrix variant, where each block is computed on distinet CPU, thereby achieving parallel execution of the matrix algorithm. This parallel approach significantly reduces memory usage, making the algorithm suitable for extremely large matrix computations. We will now demonstrate the remarkable acceleration achieved by our block version of the SKR algorithm through experimental results:
An experimental simulation was conducted to generate a dataset for the Helmholtz equation. Each corresponding matrix has dimension of 10000, with Successive Over-Relaxation (SOR) used uniformly for matrix preconditioning. running on 72 parallel MPI threads.
The upper half of the table presents the average computation time taken by both algorithms to solve each linear system, measured in seconds. The lower half details the average number of iterations required to solve each linear system. The first row of the table specifies the precision requirements for solving linear systems.
The experimental results reveal that our SKR algorithm achieves a 150-200 fold acceleration in computation time and 24-26 fold reduction in the number of iterations compared to the base- line. The significant improvement computation time mainly attributed to the effects MPI parallelization.
39

---

Published as a conference paper at ICLR 2024
Table 32: Comparative Experimets of Block Versions
1E-03
1E-05
1E-07
time(s)
GMRES
3.22
4.16
4.66
time(s)
Block SKR(ours)
0.015
0.025
0.03
iter
GMRES
5503
7074
7976
iter
Block SKR(ours)
224
269
305
E.3 DATASET VALIDITY
will explain from both theoretical and experimental perspectives that the data set generated using our SKR algorithm will not affect the training of neural operators.
E.3.1 THEORETICAL PERSPECTIVE
From the theoretical perspective, under specific eтror tolerances, our SKR algorithm does not alter the outcome of the generated dataset. As detailed in Section 5.1] our SKR algorithm is backed by rigorous convergence analysis, ensuring accurate convergence to the true solution. Especially for linear systems with high correlation, our SKR algorithm demonstrates superior theoretical convergence rate over GMRES.
In the field of numerical solutions for large matrices, there are generally true solutions without error, only numerical solutions with errors smaller than an acceptable threshold. We consider the numerical solution satisfactory when the error of the computed linear system's 'ss solution is below the acceptable threshold, which we believe accurately represents the true solution.
E.3.2 EXPERIMENTAL PERSPECTIVE
Experimentally, we found that our SKR algorithm does not affect the effectiveness of the generated dataset when training neural operators.
For the Darcy flow dataset, we simultancously generated 5256 data points with a precision 10-8 using both our SKR algorithm and GMRES, of which 5000 were designated for the training set and 256 for the test set, with the matrix size being 2500. The data in the table represent the relative error under the two-norm during the training process, where the first row indicates the number of training epochs and the rightmost column shows the final convergence error.
Table 33: Comparative Efficacy of Neural Operators Trained on Datasets Generated by Different Algorithms
0
100
200
300
400
final
GMRES 1
0.912
0.127
0.096
0.055
0.025
0.02
GMRES 2
0.927
0. 1.4
0.094
0.054
0.026
0.02
SKR(ours)
0.886
0.132
0.085
0.055
0.026
0.02
SKR(ours) 2
0.912
0.145
0.088
0.051
0.028
0.02
To mitigate the impact of randomness during the training process, we trained the neural network models twice using datasets generated by GMRES and SKR respectively to compare their train- ing dynamics. The results of this experiment demonstrate that the datasets generated by our SKR algorithm and GMRES yield identical outcomes when used to train neural operators.
Based on the aforementioned theoretical analysis and experimental validation, we can confidently state that the datascts generated by our SKR algorithm and the baseline GMRES algorithm have no impact on the training of neural operators. That is, the performance of the datasets they generate is equivalent.
40


Reviews:
[{"review_id": "VhJW47QhTR", "reviewer": "Reviewer_BjLB", "summary": "In this paper, authors proposed an iterative method called Sorting Krylov Recycling (SKR) to boost the efficiency of solving systems of linear equations, thus accelerating data generation for neural operators training. The proposed method leverages the power of Krylov subspace recycling and a sorting algorithm to arrange interrelated systems in a sequence, where adjacent systems exhibit high similarities. \nThe effectiveness of the SKR method was demonstrated through both theoretical analysis and extensive experiments. The SKR performance was compared with other state-of-the-art methods, such as GMRES, and by considering the impact of different preconditioning techniques on the convergence rate and computational efficiency. The results indicate that SKR can deliver impressive acceleration, reducing the wall clock time by factor of up to 13.9 and requiring up to 30 times fewer iterations.", "questions": "- One general question for me was whether this approach can be used beyond neural operators, and can be used for any data driven PDE solvers? I did not quite understand the necessity of using NO to generate PDE, and it would be great if authors could elaborate on this.\n- For many problems such as Poisson equation the resulting matrix, A, is symmetric and using iterative solvers for symmetric case (such as CG solver) could significantly help speed up the process. with that being said, could authors modify the algorithm by using symmetry property to make it perform better for symmetric matrices?", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "- Paper is well-written and overall easy to follow, with additional useful materials in appendices. \n- Extensive experiments and nice ablation studies to illustrate the impact of different components of the proposed algorithm\n- The code and the reproducibility details are useful for other researchers to leverage the developments for their data driven PDE solvers, or benchmark their ongoing/future studies on related topics.", "weaknesses": "- It would have been nice if authors provided more elaboration on the results for different types of PDEs, in particular for Poisson equation, for which the SKR and GMRES perform comparably to some extent. \n- In addition to time and iteration count, what about the numerical accuracy of their solvers w.r.t the grid resolution? Also, how can one interpret the slopes of plots indicating tolerance vs time and tolerance vs interations?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "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": "1VTUGZEF1L", "reviewer": "Reviewer_ZpMg", "summary": "The paper addresses a fundamental challenge in the training of Neural Operators (NOs) for solving Partial Differential Equations (PDEs) by introducing a novel algorithm called Sorting Krylov Recycling (SKR). The primary challenge is the time-consuming and computationally expensive process of generating labeled data for training NOs, which involves solving numerous systems of linear equations. The proposed SKR algorithm is designed to improve the efficiency of solving these systems and significantly accelerate data generation for NOs training.", "questions": "Could you provide more details about the specific types of PDEs or problems where the SKR algorithm is expected to have the most significant impact in terms of efficiency improvement?\n\nThe paper mentions that SKR can achieve a remarkable speedup. Are there specific parameters or settings that are critical for achieving this speedup, and are there scenarios where the speedup might be less pronounced?\n\nHave you considered potential challenges or limitations in implementing the SKR algorithm in practical applications, and are there strategies to address these challenges?\n\nBeyond the context of generating training data for NOs, are there other domains or problems where SKR's approach of optimizing linear system solutions could be applied effectively?\n\nHow does the SKR algorithm handle variations in the complexity of PDEs, and does it exhibit consistent efficiency improvements across different levels of complexity?", "limitations": "", "rating": 6, "confidence": 2, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "Innovative Algorithm: The SKR algorithm is a novel and innovative approach to accelerate the generation of training data for NOs, focusing on the efficient resolution of linear systems. It addresses a fundamental challenge in the field of data-driven PDE solvers.\n\nEfficiency Improvement: The paper demonstrates that SKR can achieve a remarkable speedup in the generation of training data, with a potential acceleration of up to 13.9 times. This significant efficiency improvement has the potential to reduce computational costs in the training of NOs.\n\nReal-World Relevance: The need for efficient data generation for NOs is of practical importance, especially in scientific domains where PDEs play a crucial role, such as climate modeling, fluid dynamics, and electromagnetism.", "weaknesses": "Complexity: While SKR is a promising algorithm, the paper does not discuss its potential complexities or challenges in practical implementation. It is essential to evaluate the algorithm's feasibility and usability in real-world scenarios.\n\nApplication Scope: The paper primarily focuses on the problem of generating training data for NOs. It would be beneficial to discuss broader applications or scenarios where SKR could be applied beyond this specific context.\n\nGeneralization: The paper does not extensively discuss the generalizability of the SKR algorithm to different types of PDEs or scenarios. It is important to assess whether SKR is applicable in a wide range of PDE-related problems.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "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": "IB5dtpFCJB", "reviewer": "Reviewer_7U2z", "summary": "This paper proposes a new method, SKR, to speed up the process of repeatedly solving a given PDE with different parameters (e.g. initial conditions, boundary conditions, right-hand side, etc.). This scenario arises when generating data to train neural operators, as this process effectively entails solving the same PDE multiple times with different input parameters. Solving a PDE often involves a potentially very large linear system, which can be solved using Krylov methods. SKR leverages the fact that the same PDE is solved across samples with different parameters to accelerate the data generation process, thanks to Krylov subspace recycling and a specific sorting algorithm.", "questions": "Is there a plan to add support for widely used PDE software such as Firedrake, FEniCS, etc.? This would enable data generation for neural operator training with just a few lines of code using SKR in high-level languages like FEniCS or Firedrake, enhancing the method's dissemination.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "This work addresses the important problem of data generation for NO training. The idea relies on the well-known Krylov subspace recycling method and, to the best of my knowledge, has not been previously proposed.", "weaknesses": "The paper does not discuss the aspect of parallelism, which raises concerns about the fairness of the benchmark. More specifically, in the traditional setting, PDE solves are independent of each other. Consequently, data generation can be achieved by splitting the PDE systems to solve across different MPI processes to reduce the cost and achieve faster data generation. However, in the SKR setting, the PDE solves are no longer independent, which undermines the claimed advantage of SKR over the classical approach.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "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": "e5giKRt9MB", "reviewer": "Reviewer_A7dr", "summary": "The paper addresses lengthy workload required for generating data for neural operators, which are one class of typical machine learning models used as surrogate models for solving partial different equations. The key idea of the paper is the use of Krylov subspace recycling method, a technique for solving a series of interrelated systems. Together with the sorting algorithm based on inherent similarities, the proposed method greatly accelerates the data generation process.", "questions": "See the weeknesses.", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "The main advantage of the paper is an oriented way of generating dataset based on Krylov subspace recycling method, whose validity is supported by rigorous theoretical analysis. The paper reveals common features of (generating process of) NO datasets. The proposed method is evaluated with a wide range of its variants and the reduction in the computational cost is shown significant for most of the cases.", "weaknesses": "Although the reduction in the computational workload of SKR is significant, it is still not clear to me that SKR is a better choice as a generator for NO dataset, since the paper lacks evaluation of the dataset generated by SKR in the context of training NOs. One of my questions in this regard is:\nHow does the performance of NO instances (e.g. FNO, DeepONet) trained on datasets generated by SKR compare to those trained on their original datasets? \nSpeaking from the practical aspect, it might be also beneficial for readers if the authors could give run-time comparison of SKR to parallelized solvers.\n\nThe presentation of the paper does not look clear to me. The followings are a couple of examples that made it difficult for me to understand the paper\n* In Section 3.1 ”The data points for these PDEs in the parameter space of the NO inputs are relatively dense” : How did the authors withdraw this conclusion? Compared to what PDEs are the inputs dense?\n* Table 5 (and its caption) is entirely confusing and hard to understand.\n* Caption in Figure 9 “Specifically, when the parameters of the neural operator are closely matched in Helmholtz equations” : how do you compare the parameters of NO to Helmonltz equations?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "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:

Decision: Accept (spotlight)