SymMaP: Improving Computational Efficiency in Linear Solvers through Symbolic Preconditioning
Authors: Hong Wang
Jie Wang
Minghao Ma
Haoran Shao
Haoyang Liu
Venue: NeurIPS 2025 poster ()
Keywords: Matrix Preconditioning
Symbolic Discovery
Linear System Solver
Primary area: machine_learning_for_sciences

Abstract: Matrix preconditioning is a critical technique to accelerate the solution of linear systems, where performance heavily depends on the selection of preconditioning parameters.
Traditional parameter selection approaches often define fixed constants for specific scenarios.
However, they rely on domain expertise and fail to consider the instance-wise features for individual problems, limiting their performance. 
In contrast, machine learning (ML) approaches, though promising, are hindered by high inference costs and limited interpretability. 
To combine the strengths of both approaches, we propose a symbolic discovery framework—namely, **Sym**bolic **Ma**trix **P**reconditioning (**SymMaP**)—to learn efficient symbolic expressions for preconditioning parameters. 
Specifically, we employ a neural network to search the high-dimensional discrete space for expressions that can accurately predict the optimal parameters. 
The learned expression allows for high inference efficiency and excellent interpretability (expressed in concise symbolic formulas), making it simple and reliable for deployment. 
Experimental results show that SymMaP consistently outperforms traditional strategies across various benchmarks.

SymMaP: Improving Computational Efficiency in Linear Solvers through Symbolic Preconditioning
Hong Wangl223 Jie Wang Minghao Ma', Haoran Shao', Haoyang Liu
University of Science and Technology of China
2 CAS Key Laboratory of Technology in GIPAS, University of Science and Technology of China 3 MoE Key Laboratory of Brain-inspired Intelligent Perception and Cognition, University of Sci- ence and Technology of China
wanghong17700mail. watt::0mmmaauuaa....mmmmmmmm www.uutt..... settc.edd.ccmm cn, jiewangzeustc.edi
Abstract
Matrix preconditioning is critical technique to accelerate the solution of linear systems, where performance heavily depends on the selection of preconditioning parameters. Traditional parameter selection approaches often define fixed constants for specific scemarios. However, they rely on domain expertise and fail to consider the instance-wise features for individual problems, limiting their performance. In contrast, machine learning (ML) approaches, though promising, are hindered by high inference costs and limited interpretability. To combine the strengths of both approaches, we propose a symbolic discovery framework-namell, Symbolic Matrix Preconditioning (Sym-MaP)-to learn efficient symbolic expressions for preconditioning parameters. Specifically, we employ a neural network to search the high- nigh-dimmesional discrete space for expressions that can accurately predict the optimal parameters. The learned expression allows for high inference efficiency and excellent interpretability (expressed in concise symbolic formulas), making it simple and reliable for deployment. Experimental results show that SymMaP consistently outperforms traditional strategies across various benchmarks
Introduction
Linear systems are foundational in machine learning, physics, engineering, and other scientific fields B1
15 16 Since analytical solutions are often unavailable, efficient numerical algorithms become essential M.. Matrix preconditioning, key technique in this domain, accelerates iterative solvers and improves computational stability [59 77. For instance, the successive over-relaxation (SOR) method optimizes convergence by integrating Gauss-Seidel iterations with weighted update scheme governed by the over-relaxation factor w [18].
32
The effectiveness of matrix preconditioning depends on key parameters, such as in SOR. Selecting can accelerate convergence, while < 1 may stabilize the process. This trade-off makes w critical parameter, directly influencing the performance of preconditioning.
Traditional parameter selection strategies often rely domain expertise to define fixed constants for specific scenarios. However, (challenge 1) different problem parameters often require distinct optimal preconditioning parameters. Traditional strategies ignore instance-wise features-specifie characteristics of individual problems, such as equation coefficients. This limits their adaptability to varying problem instances and tasks.
"Corresponding author. Our code is available at
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

In contrast, machine learning (ML) approaches hold great promise but come with other challenges. First, (challenge 2) ML inference, while efficient on GPUs, performs poorly CPU-only envi- ronments due to limited parallel processing capabilities. This is particularly problematic in linear system solver deployments, where GPU resources are often unavailable. Second, (challenge 3) the "black-box" nature of many ML techniques hinders a deeper understanding of the learned policies, raising concerns about their reliability.
In light of this, a natural solution is to combine the reliability and superior performance of these two paradigms. We propose symbolic discovery framework naamely.--ammlly...... Symbolic Matrix Preconditioning (SymMaP) to learn efficient symbolic expressions for preconditioning parameters. The framework consists of three main steps. SymMaP first begins by applying grid search identify the optimal pre- conditioning parameters based on task-specific performance metrics. Next, the framework performs a risk-seeking search the high-dimensional discrete space of symbolic expressions, evaluating the best-found symbolic expression using a risk-seeking strategy. Finally, these symbolic expres- sions can be directly integrated into the modern solvers for linear systems, significantly improving computational efficiency.
The key contributions and advantages of SymMaP are summarized follows:
We propose symbolic discovery framework, SymMal, to learn efficient symbolic expres- sions for preconditioning parameters.
SymMaP exhibits excellent generalization, making adaptable a wide range of precondi- tioning methods and optimization objectives.
The symbolic expressions derived by SymMaP are both interpretable and easy to integrate into solver environments, offering practical and transparent approach to enhancing the performance of linear system solvers.
2 Preliminaries
2.1 Matrix Preconditioning Technique
Matrix preconditioning is a technique employed to accelerate the convergence of iterative solvers and enhance the stability of algorithms. It is generally employed in solving linear systems [8]. which are typically expressed in the form: Ax The fundamental idea of preconditioning is to transform the original problem into an equivalent one with better numerical properties. Specifically, this technique involves finding preconditioner M that approximates either the inverse of A or some form conducive to iterative solutions 77]. Consequently, the original equation is transformed into M Ax Mb. There are generally two optimization objectives: 1. To accelerate the convergence of iterations by altering the spectral distribution of the matrix A. 2. To reduce the condition number of the matrix A, thereby lessening its ill-conditioning and enhancing the stability of the iteration. Some common preconditioning techniques include the Jacobi, Gauss-Seidel, SOR 168), algebraic multigrid (AMG) [61], etc.
2.2 Prefix Notation
Prefix notation is mathematical format where every operator precedes its operands, eliminating the need for parentheses required in conventional infix notation and simplifying symbolic manipulation. This representation is particularly advantageous in symbolic regression, it allows mathematical expressions to be expressed as sequences of tokens that can be easily processed by neural networks.
In this notation, operators can be unary (e.g., sin, cos) binary (e.g.. x, +),, while operands can be constants or variables B0). Each prefix expression uniquely corresponds to a symbolic tree structure, facilitating the conversion back to the original mathematical expression 29
The sequential nature of prefix notation aligns well with the architecture of recurrent neural net- works (RNNs), which process information step by step. Unlike infix notation which may require variable-length look-ahead to determine the next valid token, prefix notation allows RNN to generate expressions through an auto-regressive process where each decision is well-defined based on previous tokens, and by removing the need for parentheses, it reduces the vocabulary size of possible tokens, which greatly enhances model training efficiency.

---

12
10
Rerations
Time
100%
-104
9.%
Time Ratio
67%
Iterations
Time (s)
03%
-109
100
8.%
optimal paramister Mgp
10a
70%
Optimal Parameter tip
SymMaP
Optimal Fixe Constant
Defoult Parameter
0.00
0.25
SOR preconditioning parameter
0.50
5.75
1.00
1.25
1:30
1.75
2.00
60%
le5
le6
Tolerance (2 Norm)
Figure I: Left. Variation in iteration counts and computation times under different SOR precondition- ing parameters applied linear system from second-order elliptic PDE. Right. Ratio average computation times at various tolerances to default parameter times under different SOR parameter selection schemes, evaluated on the second-order elliptic PDE dataset.
3
Motivation
The selection of matrix preconditioning parameters significantly affects their effectiveness 77 To design appropriate algorithmic prediction parameters, we first analyze the optimization space for preconditioning parameter selection and investigate the existence of optimal parameters. Next, analyze the unique challenges present in this scenario. Finally, to address these challenges, we design symbolic discovery framework select the parameters.
3.1 Motivation for Optimizing Preconditioning Parameters
As illustrated in Figure L1 the choice of relaxation factors w significantly impacts the iteration count and computation time, when solving a second-order elliptic partial differential equation (PDE) [14] using SOR preconditioning [18]. There exists optimal parameter Wop that minimizes the computation time, with specific details available in the Appendix B.2
To further analyze the optimization space of preconditioning parameters, we evaluate the impact of various parameter selection strategies on preconditioning performance. As shown Figure I the Optimal Parameter Wap' represents the parameter that minimizes computation time in each experi- ment, serving as the theoretical upper limit of our optimization. The "Optimal Fixed Constant' refers to fixed constant that minimizes average computation time, and 'Default Parameter' corresponds to the default setting of w = 1 in the portable extensible toolkit for scientific computation (PETSc) (4]. The gap between the optimal fixed constant and the optimal parameter highlights significant potential for optimizing preconditioning parameter selection, motivating this paper. The performance of our SymMaP algorithm approaches the optimal parameter, demonstrating its accuracy in learning the optimal parameter expression.
3.2 Challenges in Predicting Efficient Preconditioning Parameters
We aim to develop universal framework for predicting efficient parameters. However, the context of solving linear systems imposes specific challenges:
(C1) Strong Generalization Capability: Real-world scientific computing scenarios vary signif- icantly. For instance, the choice of PDE grid form can lead to significant variations matrix structure (25). resulting in distinct optimal parameters. Moreover, preconditioning addresses various optimization goals, such as reducing computational time, the number of iterations, and lowering condition numbers [2]. This necessitates that parameter prediction algorithms possess robust general- ization capabilities: they should take problem scenarios and features as inputs while applying them to different preconditioning methods and optimization goals.
(C2) Computational Efficiency: linear system solver typically relies on Krylov subspace meth- ods implemented in low-level libraries optimized for CPU architectures, such PETSc 44, LA-

---

1. Training Data Generation
2. Evaluate Expressions and Train the Deep Symbolic Model datasct (+)),,)) initial parameter 0g
3. Sequential Model
PDB(a)
liess syskets A(ojs - Ko)
& Sibling
Parent
RNN
Sampled Taken
paramater
M(+)
kg
Aljaa (o)M)(
0 0n
update
prefix expression
generating expression
sequence
adsptive
epeend prrsece Aop
gnd scarch
constant optimization
RL ttttrinee
parmmeter of
calculate R(r)
(1),)) (a.
NO
4. Deployment in linear solvers
R(F) 1 or iter
max
Token Library
Symbolic log
Expression
integrated to linear solver
compiled library
Yes
expression
Figure 2: Illustration of how SymMaP discovers efficient symbolic expressions for preconditioning parameters. Part I demonstrates the acquisition of optimal parameters and dataset generation; Part illustrates the training process of the RL-based deep symbolic discovery framework; Part shows how the sequential model generates symbolic policies; Part 4 presents the deployment of symbolic expressions.
PACK 3]. Algorithms like generalized minimal residual method (GMRES) 252 and conjugate gradient (CG) 200 iteratively compute the matrix's invariant subspace, favoring single-threaded or execution modes. Preconditioning techniques aim to accelerate these solvers without significant additional computational overhead, often adopting implicit iterative formats (e.g., SOR 77) or utilizing low-cost matrix decompositions (e.g.. AMG [61). Therefore, any parameter prediction algorithms must be compatible with CPU environments and seamlessly integrated into existing algorithm libraries. At the same time, it must maintain low computational costs to preserve the performance benefits of preconditioning.
limited
(C3) Algorithmic Transparency: Algorithms in scientific computing often require rigorous analysis under mathematical theories. Opaque prediction algorithms could confuse researchers. For instance, the relaxation factor w in SOR needs to avoid being too close to or in some scenarios [. This an issue that opaque algorithms cannot avoid in advance. Moreover, interpretable algorithms can guide researchers to conduct further studies and reveal the underlying mathematical structures of problems. Therefore, these pose challenges to the transparency and interpretability of the parameter prediction algorithms.
3.3 Symbolic Discovery to Preconditioning Parameter Selection
Symbolic discovery extracts mathematical expressions from data, establishing relationships between problem features and optimal preconditioning parameters. Its integration into matrix preconditioning overcomes parameter selection challenges through a generalizable, efficient, and transparent approach.
Firstly, symbolic discovery can accommodate various types of input parameters and can specifically tailor symbolic expression learning for different preconditioning methods and optimization goals [10), thereby meeting the requirement for broad applicability in scientific computing tasks (C1). Secondly, the explicit expressions derived are computationally lightweight and can be quickly evaluated at runtime. They integrate seamlessly into existing CPU-based algorithm libraries like PETSc 4 with almost no overhead (C2). Thirdly, the symbolic discovery provides transparent and interpretable expressions [50]. allowing researchers to understand the influence of parameters within existing theoretical frameworks and identify potential numerical stability issues. This interpretability fosters trust in the algorithm's predictions and supports further theoretical exploration (C3).

---

4 Method
This study focuses enhancing the performance of parameterized preconditioners in solving linear systems derived from parameterized PDEs. Specifically, we investigate preconditioners with continuous parameters, such as the relaxation factor in SOR, while excluding those with discrete parameters like the level of fill-in ICC or ILU factorization.
We introduce novel framework, SymMal, for symbolic discovery in matrix preconditioning. As shown in Figure 2 we first obtain the optimal preconditioning parameters through a grid search to construct training dataset. Then employ RNN to generate symbolic expressions in prefix notation, which are then evaluated for their fitness. The RNN is trained using a reward function based on the performance of the generated expressions. By optimizing the RNN parameters maximize this reward function, we generate symbolic expressions that approximate the relationship between the problem's 's feature parameters (PDE parameters) and the optimal preconditioning parameters. Finally, we deploy the learned symbolic expressions into linear system solvers. The detailed steps are as follows and pseudocode is provided in the Appendix
4.1 Input Features and Training Data Generation
Input Features. In the context of solving parameterized PDEs, which frequently arise in linear systems. we consider feature parameters that characterize the equations. For instance, second-order clliptic PDE can be expressed as: Q11 Ugz + 0124zy + Uy + uz + a2 + apu f. where the coefficients 011, d12, ag2:... represent the feature parameters of PDE (see Appendix D.I for details). These feature parameters, denoted as Z4, serve as input features for the symbolic discovery process in SymMaP.
Training Data Generation. We assume that the preconditioning performance and parameters are continuous. For each linear system, we determine optimal preconditioning parameters through an adaptive grid search. Using SOR preconditioning as an example, we optimize the relaxation factor w within [0. 2] minimize computation time (or condition number). The search process involves: 1. Conducting an initial coarse grid search (step size: 0.01) to evaluate computation time (or condition number) for each w. 2. Identifying candidate regions with optimal performance. Performing refiedd grid search (step size: 0.001) within these regions.
This process yields the required training dataset, where each data point contains: 1. the problem feature parameters I;. 2. the optimal preconditioning parameters ys. = 1, 2,...... n,, and n is the number of data, typically set test set.
1200, with 1000 allocated for the training set and 200 for the
4.2 The Generation of Symbolic Expressions
Token Library. For SymMaP, we define the library L of mathematical operators and operands as {+, x, + -qrr,, exp, log, pow, 1.0). Although other operators such as poly, sin and cos are frequently used (62), decided exclude them because they offer limited explanatory power in matrix preconditioning and significantly increase the time and memory consumption during training.
After converting the mathematical expressions into prefix notation, use this tokenized represen- tation as pre-order traversal of the expression tree [69]. In each iteration, the RNN receives a pair consisting of parent node and a sibling node as inputs. Then the RNN outputs categorical distribution over all possible next tokens. The parent node refers to the last incomplete operator that requires additional operands to form a complete expression. The sibling node, in the context of binary operator, represents the operand that has already been processed and incorporated into the expression. In cases where parent sibling node is applicable, they are designated as empty nodes. This structured input method enables the RNN to maintain contextual awareness and effectively predict the sequence of tokens that form valid mathematical expressions.
The Sequential Model. During the generation categorical distribution for each "next token" each step. This distribution is represented as vector (1))) where denotes the i-th step and 0 represents the parameters of the RNN. The elements of the vector correspond to the probabilities of each token, conditioned on the previously selected tokens in
single symbolic expression, the RNN emits a
5

---

the traversal [6]
(I)
Here, Tr denotes the index of the token selected at the i-th step. The probability of generating the entire symbolic expression T is then the product of the conditional probabilities of all tokens [46 30]
p = II 00 (T)).
Optimization of Constants: The library L incorporates a 'constant token,' which allows for the inclusion of various constant placeholders within sampled expressions. These placeholders serve as the parameters & in the symbolic expression. We seek to find the optimal values of these parameters by maximizing the reward function: E arg max; R(T: {), utilizing nonlinear optimization method. This optimization is executed within each sampled expression as an integral part of computing the reward, prior to each training iteration.
4.3
The Reward Function
Once a symbolic expression is fully generated (i.e., the symbolic tree reaches all its leaf nodes), we evaluate its fitness by calculating the normalized root-mean-square error (NRMSE), metric frequently used in genetic programming symbolic discovery [54]. is defined as NRMSE
oy -V - -0:0 where ý = -(;;) is the predicted value for the i-th sample, x; is the problem feature parameter, Yi is the optimal preconditioning parameter, Øy is the standard deviation of the target values y, and the number of data. To bound this fitness value between 0 and 1, we apply squashing function: R(T) 1/(1 + NRMSE). Our objective is to maximize R(T)) thereby minimizing the NRMSE and improving the accuracy of the generated expressions.
4.4 The Training Algorithm
Although the objective function is well-defined, is important to note that R(T) is not deterministic value but random variable dependent the RNN's parameters O. Therefore, the key challenge is to establish an appropriate criterion for evaluating this random variable and then apply gradient-baser optimization methods accordingly.
Risk Policy, is natural to consider the expectation of the reward function, i.e., (r=)p(r:0)) [((()) as the objective function to optimize. We can apply the "log-integral" trick [66 and obtain
Ve (o [[([)]]][[]]]]]]]]]][]]] (E- [R(T) V p(r; 0)].
(3)
Thus, even though the expectation of the reward function is not directly differentiable with respect to 0, we can approximate the gradient using the sample mean.
In the context of symbolic regression, model performance is often driven by a few exceptional results that outperform others by a significant margin 446
58
With this in mind, we adopt a risk-seeking
policy, which aims to maximize:
J(0, 0 [((T) [ R Q(0, E)]-
Here, E is the risk factor, typically = 0.05. ((0,=) ) is the (1 E) -quantile of the reward distribution under parameter i.e.
@ e) = inf E RCCDF(((T): 2 1 11--) E),
(5)
where (DF(R(T):: 0) refers to the cumulative distribution function. From this, the gradient of J(0,s) can be derived as [46]:
(0,(0, = ())p() (A(()) Q(0,e)) V Vologp(T: R(T)) Q(0,e)
(6)
This gradient be estimated using Monte Carlo sampling:
(e) ~
Σ - logp wow.olg.....oo................oo [R((())=(10.5)
EN
1=1
@(0, e) is the empirical (1 Ei--an-uanttle of the reward function. By concentrating on the top per- centile of samples, SymMaP emphasizes optimizing the best-performing solutions in preconditioning thereby obtaining the optimal symbolic expressions for preconditioning parameters.

---

4.5 Deployment in Linear System Solver
After the training process, we obtained a symbolic expression for predicting the preconditioning parameter. The learned formula is exceptionally concise and incurs minimal computational cost. Therefore, we directly compile the learned policy into lightweight shared object using a simple script and then integrate it into the linear system solver package (e.g.. PETSc).
Experiments
We conducted comprehensive experiments to evaluate the SymMaP framework, organized into three primary sections: 1. Assessment of three different preconditioners and optimization goals across various datasets to determine the effectiveness of SymMaP, 2. Analysis of associated computational cost and the interpretability of the learned symbolic expressions, Ablation studies.
Preconditioners:
We considered three different preconditioners and various optimization metrics: 1. SOR preconditioner with the relaxation factor w [88] SSOR preconditioner with the relaxation factor w [88; 3. AMG preconditioner with the threshold parameters 0r [6]].
Datasets: We investigated linear systems derived from five distinct PDE classes: 1. Darey Flow Prob-
lems 35 2. Second-order Elliptic PDEs 14] 3. Biharmonic Equations (5). 4. Thermal Problems [65] 5. Poisson Equations (65 All cases except biharmonic equations yield symmetric matrices. No- tably, the non-symmetric matrices from biharmonic equations are incompatible with SSOR and AMG preconditioning techniques. Additionally, we conducted partial tests the Lippnann-Schwinger equation (a numerical integration problem) and Markov chains (an optimization problem). Details of these experiments be found in Appendix [..
Baselines: We compared SymMaP against various parameter selection methods for preconditioning Specifically, the comparison involved the following scenarios: 1. No matrix preconditioning, 2. Default parameters in PETSc 4 Fixed constants, Optimized fixed constants.
Experiment Settings: All preconditioning procedures were uniformly implemented using the C- based PETSc library 4 to maintain evaluation consistency. The experiments were conducted within PETSc's linear solver framework (GMRES and CG) [18], with condition numbers computed through built-in function ISPComputeEExtreeeSiiinn arValues.
Details on preconditioners, the mathematical forms of datasets, tolerance metric (relative residual), and the runtime environment are available Appendices D.I D.3 and D.2] respectively. Infor- mation on the generation of training datasets for the following experiments and parameters of the SymMaP algorithm are outlined in Appendices D.44 DD4 and D.5 The generated dataset and training time are available in Appendix D.6 For introduction to related work, see Appendix
5.1 Main Experiments
In these experiments, as shown in Tables we optimized relaxation factors w in both SOR and SSOR preconditioning, and threshold parameters 0, AMG preconditioning. For SOR and SSOR. we identified w values that minimize computation time, forming the training dataset for SymMaP to learn symbolic expressions that optimize computational times for solutions. Similarly, for AMG, we
Table I: Comparison of average computation times (seconds) for SOR with different w selections, and tolerance is 1e-7. SymMaP and 2 are the two learned expressions that achieved the highest reward function scores, with the best-performing method highlighted in bold.
No peecondition
PETSe default
Fixed constant
Fixed coustant
Dataset
Optimal constant
SymMaP
SymMaP
0.2
1.8
Biharmonic
4.2 x 10a
7.67
2.04
4.86
1.60
1.31
1.24
1.26
Dancy Flow
1.0 x 10*
33.1
13.5
17.5
9.91
9.54
8.50
8.60
Elliptic PDB
4.0 x 10a
31.3
21.0
21.4
17.5
16.6
15.8
16.3
Poisson
10a
4.12x10
1.95xx10
1.95x10 2
1.38x10
1.35x10
1.36x10
Thermal
2.8 102
2
5.98x10-2
2.07x10-7)
1.18x10
5.94x10-2
5.76x1102
5.91x10-2

---

Table 2: Comparison of average computation times (seconds) for SSOR with different selections, and tolerance is 1e-7. SymMaP and 2 are the two learned expressions that achieved the highest reward function scores, with the best-performing method highlighted in bold.
No precondition
PETSc default
Fixed constant
Fixed constant
Dataset
Matrix
Optimal constant
SymMaP
SymMaP
wm1
m 0.2
1.8
Dancy Flow
4.9 x 10*
4.18
0.488
0.757
1.09
0.448
0.412
0.523
Elliptic PDE
4.0 x 104
23.9
10.5
14.7
8.72
8.68
7,70
7.74
Poisson
23 x 10a
1.12x10
1,02x
10 2
1.93x10 2
x,91x10 4
x0
9.10
9.92x10
Thermal
x 10a
2.34x10-1
2.69x 10
5.08x10-7
9.87x10-2
2.24x10-2
2.13x1102
Table 3: Comparison of average condition numbers for preconditioned matrices using different threshold parameter Or selections in AMG. SymMaP I and 2 are the two learned expressions that achieved the highest reward function scores, with the best-performing method highlighted in bold.
No precondition
PETSc default
Fixed constant
Fixed constant
Dataset
Marrix size
Optimal constant
SymMaP
SymMaP 2
By
9p 0.8
Darcy Flow
1.0 x 10*
752862
8204
19146
11426
7184
4824
S786
Elliptic PDE
4.0 x 10
6792
184.6
205.4
212.5
182.8
168.8
170.3
Poisson
1.0 x 10
1242
4.55
68.85
68.85
4.55
3.72
3.72
Thermal
288 x 10
7325
11.9
627.2
627.2
9.91
9.71
9.71
selected Or values that minimize the condition number of preconditioned matrices, Partial symbolic expressions can be found in Appendix [E.]]
Experimental results indicate that SymMaP consistently outperforms all others across all experimental tasks. For SOR, Table I shows that SymMaP reduces computation times by up 40% compared to PETSc's default settings and by 10% against the optimal constants. In SSOR, Table shows that cuts computation time and iteration counts by up 27%, over PETSc's defaults, and achieves reductions of 11% in time compared to optimal constants. For AMG, Table 3 shows that SymMaP lowers the condition number by up to 40% relative PETSc's defaults and 32% against the optimal constants. These results highlight SymMaP's ability to effectively derive high- high-performance symbolic expressions for various preconditioning parameters, showcasing its broad applicability and strong generalization across different preconditioning tasks.
Additionally, we analyze SymMaP's ties, generalization performance (e.g.. performance outside the training parameter range and under varying matrix sizes, geometries), error reduction trajectories, and statistical characteristics of some experiments (c.g., medians, quartiles, etc.), which can be found in Appendices E.2 E.6 One of the symbolic expressions derived by SymMaP exhibits significant ad advantages over existing learning-bast optimize only a single preconditioner), interpretabilityy and compatibility with CPU-based scenarios. Even without considering these advantages, we compared SymMaP with existing learning- based algorithms. SymMaP consistently demonstrated superior performance across all experiments, with detailed results provided in Appendix E.7
prediction capabilities, custom metric capabili-
E.3 E.4 E.5
algorithms in terms of applicability (as current tearning-basee methods often
5.2 Comparison with Neural Network Performance
To evaluate the deployment overhead and prediction per- formance of SymMaP, we compared it with a basic multi- layer perceptron (MLP) architecture. The MIP implemen- tation consists of three fully connected layers, taking PDE parameters input and generating preconditioning param- eters as output. We employed ReLU activation functions and trained the model using mean squared error (MSE) between predicted and optimal parameters as the loss func- tion. Both symbolic expression and MLP were executed in a CPU environment to simulate a modern solver envi- ronment.
Table Comparison of the runtime re- quired for symbolic expression and MLP predict the SOR relaxation factor and the subsequent average solution time, us- ing the Darcy flow dataset with a matrix size of 102 and tolerance is le-5,
Runtime (s)
Solution time (s)
MLP Symbol
x 10
7.1x 10
1.x10
7.1x10
As shown in Table 4 the runtime of symbolic expressions learned by SymMaP was only 20% of that of the MLP,
8

---

primarily due to the poor performance of neural networks in pure CPU environment, highlighting SymMaP's computational efficiency. Furthermore, the average solution times for parameters predicted by both symbolic expressions and MLP were closely matched. This demonstrates that symbolic expressions possess equivalent expressive capabilities to neural networks this scenario, effectively approximating the optimal parameter expressions.
5.3 Interpretable Analysis
In Table we report subset of the learned symbolic expressions, with the mathematical significance of the related symbols detailed in Appendix E.8 More symbolic expressions can be found in Appendix E.T These symbolic ex-
Table 5: Partial symbolic expressions learned from some of the main experiments.
Dutase
Symbolic expression
pressions are notably more concise and selective, not utilizing all candidate parameters and sym- bols, which aids researchers in analyzing their underlying relationships.
SOR
Ritarmonic
1.0/(4.0 1.0/2g) 1.0/z1 1.0 + 1.0/2gg + 1.0 1.0/(zz 4.0))
SOR
Elliciic PDB
SOR
Darcy Flow
1.0 1.0/(=4 1.0)
Ellipiic PUE
1.0 1.0//a/ 1.2)
AMG
Ellipiic PDE
(21 1)/7
For instance. in the context of SOR and SSOR preconditioning, empirical evidence suggests that smaller relaxation factors should be chosen when diagonal components are relatively small. Our experimental findings corroborate this: for the second-order elliptical PDE dataset, the symbolic expressions derived for SOR and SSOR preconditioning depend solely on with larger 22 values leading to smaller predicted relaxation factors, exemplified by 1.0 (22+1.2) Here, 22 represents the coupling coefficient of the elliptical PDE, which directly influences the relative size of the non- diagonal components of the generated matrix, whereas other coefficients have minimal impact. As the coupling coefficient increases, the relative numerical value of the non-diagonal components increases, and the diagonal components reduce correspondingly, aligning with empirical observations.
These experimental outcomes demonstrate that SymMaP can derive interpretable and efficient symbolic expressions for parameters, further aiding researchers in understanding and exploring the underlying mathematical principles.
5.4 Ablation Experiments
We conducted an ablation study using SymMaP to evaluate the impact of dif- ferent mathematical operator selections, the main
Table 6: Comparison of the impact of the choice of math- ematical operators on preconditioning and training time. The first column lists the selected operators, the second column shows the condition numbers of preconditioned matrices derived from AMG parameter predictions on the Darcy flow dataset (lower is better), and the third column displays SymMaP training times.
as described in Table experiments, We utilized the operator set (,,-,,,+ +, sqrt, exp, log, pow, 1.0] listed in the second row.
The results indicate that this selection of operators achieves balance between pre- dictive performance and training time ef- ficiency, meeting our expectations. A de- tailed hyperparameter sensitivity analysis concerning learning rate, batch size, and dataset size provided in Appendix E.9 Furthermore, comparative experiments against alternative parameter search meth- ods and other interpretable algorithms presented in Appendices E. and E. Il
Functionset
Condition number
timels)
thy X6, poly
15351
sin, esp,
7086.9
703.17
M. sggt, exp. log, cos, pow, 1.0
7172.6
635.82 32
th -i 1.0, pow
7241.5
703.26
+yy diy pow, +hx Xi, th, 1.0
746.30
7301.4
702.46

---

Limitations and Conclusions
In this paper, we propose SymMaP, deep symbolic discovery framework designed for predicting efficient matrix preconditioning parameters. Experiments show that SymMaP can predict high- performance parameters and is applicable across a variety of preconditioning and optimization objectives. Additionally, SymMaP is easy to deploy with virtually no computational cost.
Future work will focus on optimizing preconditioning for specific matrix structures, such as symmetric and upper triangular matrices. We also aim to analyze the mathematical significance of the learned symbolic expressions from theoretical perspective, such exploring the impact of problem features on the solution process through pseudospectral analysis. Furthermore, plans to extend SymMaP to additional preconditioning methods (c.g., ILU, ICC) are underway. We are confident in the symbolic model's 's immense potential for broad real-world applications, especially in matrix preconditioning.
Acknowledgements
The authors would like thank all the anonymous reviewers for their insightful comments and valuable suggestions. This work supported by the National Key R&D Program of China under contract 20222D0119801, and the National Nature Science Foundations China grants U23A20388, 62021001.
10

---

References
[1] Ravi Agarwal. Difference equations and inequalities: theory, methods, and applications. CRC Press, 2000.
[2] Sivaram Ambikasaran, Carlos Borges, Lise- Marsie Imbert-Gerard and Leslie Greengard. Fast, adaptive. high-order accurate discretization of the lippmann -enmaann-schwngee equation in two dimen- sions. SIAM Journal Scientific Computing, 3173)A1770-A1187,, 2016.
[3] E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Green- baum, S. Hammarling, A. McKenney, and D. Sorensen. LAPACK Users' Guide. Society for Industrial and Applied Mathematics, Philadelphia, PA, third edition, 1999.
[4] Satish Balay, Shrirang Abhyankar, Mark F. Adams, Steven Benson, Jed Brown, Peter Brune, Kris Buschelman, Emil M. Constantinescu, Lisandro Dalcin, Alp Dener, Victor Eijkhout, Jacob Faibussowitsch, William D. Gropp, Václav Hapla, Tobin Isaac, Pierre Jolivet, Dmitry Karpeev, Dinesh Kaushik, Matthew G. Knepley, Fande Kong, Scott Kruger, Dave A. May, Lois Curfman McInnes, Richard Tran Mills, Lawrence Mitchell, Todd Munson, Jose E. Roman, Karl Rupp, Patrick Sanan, Jason Sarich, Barry F. Smith, Stefano Zampini, Hong Zhang, Hong Zhang, and Junchao Zhang. PETSc Web page. https://petsc. orgp 2024.
[5] Igor Barrata, Joseph P Dean, Jørgen S Dokken, Michal Habera, Jack HALE, Chris Richardson, Marie E Rognes, Matthew W Scroggs, Nathan Sime, and Garth N Wells. Dolfinx: The next generation fenies problem solving environment. 2023.
[6] Lipman Bers, Fritz John, and Martin Schechter. Partial differential equations. Mathematical Soc., 1964.
American
[7] Ke Chen. Matrix preconditioning techniques and applications. Number 19, Cambridge University Press, 2005.
[8] Xiangning Chen, Chen Liang. Da Huang, Esteban Real, Kaiyuan Wang, Hieu Pham, Xuanyi Dong, Thang Luong, Cho-Jui Hsich, Yifeng Lu, et al. Symbolic discovery of optimization algorithms. Advances neural information processing systems, 36, 2024.
[9] Philippe G Ciarlet and Peerre-Araud Raviart. A mixed finite element method for the biharmonic equation. In Mathematical aspects of finite elements in partial differential equations, pages 125 145. Elsevier, 1974.
[10] Miles Cranmer, Alvaro Sanchez Gonzalez, Peter Battaglia, Rui Xu, Kyle Cranmer, David Spergel, and Shirley Ho. Discovering symbolic models from deep learning with inductive biases. Advances neural information processing systems, 33:17429 17442, 2020.
[11] James Demmel. Applied numerical linear algebra. SIAM, 1997.
[12] Huanshuo Dong, Hong Wang, Haoyang Liu, Jian Luo, and Jie Wang. Accelerating pde data generation via differential operator action in solution space. arXiv preprint 202502240..05957, 2024.
[13] Tobin A Driscoll, Nicholas Hale, and Lloyd N Trefethen. Chebfun guide, 2014.
[14] Lawrence C Evans. Partial differential equations, volume 19. American Mathematical Society, 2022.
[15) Wenhan Gao and Chunmei Wang. Active learning based sampling for high-dimensional nonlinear partial differential equations. Comput. Phys., 475:111848 11848, February 2023.
[16] Wenhan Gao, Ruichen Xu, Yuefan Deng, and Liu. Discretization-inaaianeee? on the discretization mismatch emrors in neural operators. The Thirleenth International Conference Learning Representations, 2025.
[17] Roland Glowinski and Olivier Pironneau. Numerical methods for the first biharmonic equation and for the two-dimensional stokes problem. SIAM review, 2172):1 1979.
[18] Gene H Golub and Charles F Van Loan. Matrix computations. JHU press, 2013.
11

---

[19] Markus Götz and Hartwig Anzt. Machine learning coaning--ided aided numerical linear algebra: Convo- lutional neural networks for the efficient preconditioner generation. In 2018 IEEE/ACM 9th Workshop on Latest Advances in Scalable Algorithms for Large-Scale Systems (scalA), pages 49-56. IEEE, 2018.
[20] Anne Greenbaum. Iterative methods for solving linear systems. SIAM, 1997.
[21] Daniel Greenfeld, Meirav Galun, Ronen Basri, Irad Yavneh, and Ron Kimmel. Learning to optimize multigrid pde solvers. In International Conference on Machine Learning, pages PMLR, 2019.
[22] Liu Haoyang, Wang Jie, Cai Yuyang, Han Xiongwei, Kuang Yufei, and HAO Jianye. Optitree: Hierarchical thoughts generation with tree search for LLM optimization modeling. In The Thirty-ninth Annual Conference on Neural Information Processing Systems, 2025.
[23] Paul Häusner, Ozan Öktem, and Jens Sjölund. Neural incomplete factorization: learning preconditioners for the conjugate gradient method. Transactions Machine Learning Research, 2024.
[24] Jun-Ting Hsieh, Shengjia Zhao, Stephan Eismann, Lucia Mirabella, and Stefano Ermon. Learn- ing neural pde solvers with convergence guarantees. arXiv preprint (90i:0:996..1220. 2019.
[25] Claes Johnson. Numerical solution of partial differential equations by the finite element method. Courier Corporation, 2009.
[26] Seid Koric and Diab W Abueidda. Data-driven and physics- informedd deep learning operators for solution of heat conduction equation with parametric heat source. International Journal of Heat and Mass Transfer, 203: 123809, 2023.
[27] Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Learning maps between function spaces. arXiv preprint (0.iv::2/8.0844:: 2021.
(28) Yufei Kuang, Jie Wang, Haoyang Liu, Fangzhou Zhu, Xijun Li, Jia Zeng, HAO Jianye, Bin Li and Feng Wu. Rethinking branching on exact combinatorial optimization solver: The first deep symbolic discovery framework. In The Twelfth International Conference on Learning Representations, 2024.
[29] Guillaume Lample and François Charton. Deep learning for symbolic mathematics. arXiv preprint arXiv:1912.01141 2019.
[30] Mikel Landajuela, Brenden K Petersen, Sookyung Kim, Claudio P Santiago, Ruben Glatt, Nathan Mundhenk, Jacob Pettit, and Daniel Faissol. Discovering symbolic policies with deep reinforcement learning. In International Conference on Machine Learning, pages 5979-5988 PMLR, 2021.
[31] Steven Leon, Lisette G De Pillis, and Lisette G De Pillis. Linear algebra with applications. Pearson Prentice Hall Upper Saddle River, NJ, 2006.
[32] Randall J Vequee Finite difference methods for ordinary and partial differential equations: steady-state and time-dependent problems. SIAM, 2007.
[33] Xijun Li, Fangzhou Zhu, Hui-Ling Zhen, Weilin Luo, Meng Lu, Yimin Huang, Zhenan Fan, Zirui Zhou, Yufci Kuang, Zhihai Wang, et al. Machine learning insides optverse ai solver: Design principles and applications. arXiv preprint arXiv:2401.05960 2024.
[34] Yichen Li, Peter Yichen Chen, Tao Du, and Wojeiech Matusik. Learning preconditioners for conjugate gradient pde solvers. In International Conference Machine Learning, pages 19425-19439. PMLR, 2023.
[35] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differen- tial equations. arXiv preprint arXiv:2010.08899, 2020.
12

---

[36] Haoyang Liu, Yufei Kuang, Jie Wang, Xijun Li, Yongdong Zhang. and Feng Wu. Promoting generalization for exact solvers via adversarial instance augmentation, 2023.
[37] Haoyang Liu, Jie Wang, Zijie Geng, Xijun Li, Yuxuan Zong, Fangzhou Zhu, Jianye HAO, and Feng Wu, Apollo-MILP: An alternating prediction-correction neural solving framework for mixed-integer linear programming. In The Thirteenth International Conference on Learning Representations, 2025.
[38] Haoyang Liu, Jie Wang, Wanbo Zhang, Zijie Geng, Yufei Kuang. Xijun Li, Bin Li, Yongdong Zhang, and Feng Wu. MILP-studio: MILP instance generation via block structure decompo- sition. In The ehirtt-eighth Annual Conference on Neural Information Processing Systems, 2024.
[39] Hongyu Liu, Haoyang Liu, Yufci Kuang, Jie Wang, and Bin Li. Deep symbolic optimization for combinatorial optimization: Accelerating node selection by discovering potential heuristics. In Proceedings of the Genetic and Evolutionary Computation Conference Companion, pages 2067-2075, 2024.
[40] 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 Engineering, 393:114778, 2022.
[41] Jian Luo, Jie Wang. Hong Wang, Zijie Geng, Hanzhu Chen, Yufei Kuang, et al. Neural krylov iteration for accelerating linear system solving. Advances in Neural Information Processing Systems, 37: 1886661288667,, 2024.
[42] Ilay Luz, Meirav Galun, Haggai Maron, Ronen Basri, and Irad Yavneh. Learning algebraic multigrid using graph neural networks. In International Conference on Machine Learning, pages 6489-6499. PMLR, 2020.
[43] Qitan Lv, Tianyu Liu, and Hong Wang. Exploiting edited large language models as general scientific optimizers. arXiv preprint arXiv:2503.0962) 2025.
[44] Daniel J Mankowitz, Andrea Michi, Anton Zhernov, Marco Gelmi, Marco Selvi, Cosmin Paduraru, Edouard Leurent, Shariq Iqbal, Jean- nannBBatistte Lespiau, Alex Ahern, et al. Faster sorting algorithms discovered using deep reinforcement learning. Nature, 618(7964):277-263 2023.
[45] RF Pawula. Approximation of the linear boltzmann equation by the fokker-planck equation. Physical review, 162(1):186, 1967.
[46] Brenden K Petersen, Mikel Landajuela, T' Nathan Mundhenk, Claudio P Santiago, Soo K Kim, and Joanne Kim. Deep symbolic regression: Recovering mathematical expressions from data via risk-seeking policy gradients. arXiv preprint arXiv: 191220487, 2019.
[47] Riccardo Poli, WB Langdon, and NF McPhee. A field guide to genetic programming (with contributions by jr k2z8))2008). Published via http://ull.. com, 2008
[48] Tianle Pu, Zijie Geng, Haoyang Liu, Shixuan Liu, Jie Wang. Li Zeng, Chao Chen, and Changjun Fan. RoME: Domain-robust mixture-of-experts for MILP solution prediction across domains. In The Thirty-ninth Annual Conference on Neural Information Processing Systems, 2025.
[49] Md Ashiqur Rahman, Zachary E Ross, and Kamyar Azizzadenesheli. U-no: U-shaped neural operators. arXiv preprint arXiv:2204. /27, 2022.
[50] Cynthia Rudin. Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead. Nature machine intelligence, 1(5):206-215, 2019.
[51] John W Ruge and Klaus Stüben. Algebraic multigrid. Multigrid methods, pages 73-130. SIAM, 1987.
[52] 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.
13

---

[53] Yousef Saad. Iterative methods for sparse linear systems. SIAM, 2003.
[54] Michael Schmidt and Hod Lipson. Distilling frree-form natural laws from experimental data. science, 324(5923)881-85, 2009.
[55] 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 ar0iv 1374, 2018.
[56] Rita Stanaityte. ILU and Machine Learning Based Preconditioning for the Discretized Incom- pressible Navier-Stokes Equations. PhD thesis, University of Houston, 2020.
[57] Ali Taghibakhshi, Scott MacLachlan, Luke Olson, and Matthew West. Optimization -aasad -bsed algebraic multigrid coarsening using reinforcement learning. Advances in neural information processing systems, 14:12129-12140 2021.
[58] Aviv Tamar, Yonatan Glassner, and Shie Mannor. Policy gradients beyond expectations: Conditional value-at-risk. Citeseer, 2015.
[59] Lloyd N Trefethen and David Bau. Numerical linear algebra. SIAM, 2022.
[60] Vladislav Trifonov, Alexander Rudikov, Oleg Iliev, Yuri M Laevsky, Ivan Oseledets, and Ekaterina Muravleva. Learning from linear algebra: A graph neural network approach to preconditioner design for conjugate gradient solvers. arXiv preprint arXiv:2405.1555) 2024.
[61] Ulrich Trottenberg. Cornelius W Oosterlee, and Anton Schuller. Multigrid. Elsevier, 2000.
[62] Silviu-Marian Udrescu and Max Tegmark. Ai feymman: A physics-inspired method for symbolic regression. Science Advances, 6 6(16),eaay2631, 2020.
[63] Gennadi Vainikko. Fast solvers of the lippmann-schwinge equation. In Direct and inverse problems of mathematical physics, pages 423-440. Springer, 2000.
[64] Franz. J Vesely et al. Computational Physics. Springer, 1994.
[65] Hong Wang. Zhongkai Hao, Jie Wang, Zijie Geng, Zhen Wang, Bin Li, and Feng Wu. Ac- celerating data generation for neural operators via krylov subspace recycling, arXiv preprint arXiv:2401.09516 2024.
[66] Ronald J Williams. Simple statistical gradient-foloowing algorithms for connectionist reinforce- ment learning. Machine learning, 8:229-256, 1992.
[67] Zherui Yang, Zhehao Li, Kangbo Lyu, Yixuan Li, Tao Du, and Ligang Liu. Learning sparse approximate inverse preconditioners for conjugate gradient solvers gpus. In Advances in Neural Information Processing Systems, 2025.
[68] David Young. Iterative methods for solving partial difference equations of clliptic type. Trans- actions of the American Mathematical Society, 76(1):92-111, 1954.
[69] Wojeiech Zaremba and Ilya Sutskever. Learning execute. arXiv preprint 0rXiv:1410.4615. 2014.
[70] Enrui Zhang, Adar Kahana, Eli Turkel, Rishikesh Ranade, Jay Pathak, and George Em Karni- adakis. hybrid iterative numerical transferable solver (hints) for pdes based on deep operator network and relaxation methods. arXiv preprint arXiv:2208. 13273, 2022.
14

---

NeurIPS Paper Checklist
The checklist is designed encourage best practices for responsible machine learning research, addressing issues of reproducibility, transparency, research ethics, and societal impact. Do not remove the checklist: The papers not including the checklist will be desk rejected. The checklist should follow the references and follow the (optional) supplemental material. The checklist does NOT count towards the page limit.
Please read the checklist guidelines carefully for information on how to answer these questions. For each question in the checklist:
You should answer [Yes], [No] or [NA].
[NA] means either that the question is Not Applicable for that particular paper or the relevant information is Not Available.
Please provide short (1-2 sentence) justificatin right after your answer (even for NA).
The checklist answers are an integral part of your paper submission. They are visible to the reviewers, area chairs, senior area chairs, and ethics reviewers. You will be asked to also include it (after eventual revisions) with the final version of your paper, and its final version will be published with the paper.
The reviewers of your paper will be asked to use the checklist as one of the factors their evaluation. While "[Yes] is generally preferable to "[No]" ", is perfectly acceptable to answer "[No] pprovided proper justification is given (e.g., "error bars are not reported because would be too computationally expensive" or "we were unable to find the license for the dataset we used"). In general, answering "[No] or "[NA] is not grounds for rejection. While the questions are phrased in binary way, we acknowledge that the true answer is often more nuanced, so please just use your best judgment and write justification to elaborate. All supporting evidence can appear either in the main paper or the supplemental material, provided in appendix. If you answer [Yes] to question, in the justification please point to the section(s) where related material for the question can found.
IMPORTANT, please:
Delete this instruction block, but keep the section heading "NeurIPS Paper Checklist",
Keep the checklist subsection headings, questions/answers and guidelines below.
Do not modify the questions and only use the provided macros for your answers.
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification:
The claims are put in the abstract and Section |
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize other settings.
It is fine include aspirational goals as motivation long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors? [Yes]
Answer:
Justification: For further details, please refer to the "Limitations and Conclusions" section, found in Section
15

---

Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g.. independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or speech-to-text system might not be used reliably provide closed captions for online lectures because it fails handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, authors should discuss possible limitations of their approach to address problems of privacy and faireess
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community, Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full of assumptions and a complete (and correct) proof?
Answer: [NA]
Justification:
The paper does not include theoretical results.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors encouraged provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
Justification: For detailed information, please refer to the "Method" section Section 4 Complete pseudocode can be found in Appendix C Specific experimental parameters
discussed Appendix D.2 Appendix D.4 and Appendix D.5
Guidelines:
The answer NA means that the paper does not include experiments.
16

---

If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. w
For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurlPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how to reproduce that algorithm. (b) the contribution is primarily a new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or a way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification:
We provide the instruction and code in supplemental material.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (nttps://nips.co cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
/ 7/nips. /ttppp.cc/pbliccguideeee http://oeeuuiiissooooooooooo
for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
17

---

Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Question:
Answer: [Yes]
Justification: Specific experimental parameters are discussed in Appendix D.2 Ap- pendix D.4 and Appendix D.5
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: For a detailed comparison, please refer to the experimental sections in Sec- tion 5.1 and Appendix E.6
Guidelines:
The answer NA means that the paper does include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
It should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of exccution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: Specific experimental environments can be found in Appendix D.2
Guidelines:
The answer NA means that the paper does not include experiments.
18

---

The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics http://pwwpppaae......o/ppp//baa
Answer: [Yes]
Justification:
We have checked the NeurIPS code of ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
the authors answer No, they should explain the special cireumstances that require deviation from the Code of Ethics.
The authors should make to preserve anonymity (e.g., if there is special consid- eration due to laws or regulations their jurisdiction..
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
The work presented in this paper aims to develop a deep symbolic discovery framework designed for predicting efficient matrix preconditioning parameters, focused solely on the field of scientific computing. It does have any direct positive negative social impacts.
Justification:
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious unintended (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML..
11. Safeguards
Question: Does the paper describe safeguards that have been put in place for responsible release of data models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
19

---

Answer:
[NA]
Justification:
We have checked this and confirmed the paper poses no such risks.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging. and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g.. code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: The original papers or URLs of the codes, models, and data sets used this article have been cited in the paper.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website). the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the bataasess.ppprsstttt paperswi theode. com/datasets package should be provided. For popular datasets,
has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. If this information is available online, the authors are encouraged reach out to the asset's 's creators.
13. New assets
Question: Are new assets introduced in the paper well documented and the documentation provided alongside the assets?
Answer: [Yes]
Justification:
We provide the instruction in supplemental material about our code and data.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL or include anonymized zip file.
14. Crowdsourcing and research with human subjects
20

---

Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification: The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurIPS Code of Ethies, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country or institution) were obtained?
Question:
Answer: [NA]
Justification: The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, formatting purposes and does not impact the core methodology, scientific rigorousness, or originality of the research, declaration is not required.
Answer: [NA]
Justification: This paper only utilizes LLMs for writing and editing.
Guidelines:
The answer NA means that the method development in this research does not involve LLMs as any important, original, or non-standard components.
Please refer to our L.MM policy fbttps://neuria ce / for what should or should be described.
21

---

A Related work
A.1 Machine Learning for Algorithm Discovery
Machine learning has the potential to uncover implicit rules beyond human intuition from training data, enabling the construction of algorithms that outperform handcrafted programs. Approaches to algorithm discovery in machine learning encompass symbolic discovery, program search, and more. Specifically, program search focuses on optimizing the computational processes of algorithms. For example, [44] explores the discovery of faster sorting algorithms, while [ investigates efficient optimization algorithms.
In contrast, symbolic discovery aims to search within the space of small mathematical expressions rather than computational streams [46] 30). This approach is analogous to an extreme form of model distillation, where knowledge extracted from black-box neural networks is distilled into explicit mathematical expressions. Traditional methods for symbolic discovery have relied on evolutionary algorithms, including genetic programming [47.. Recently, deep learning has emerged powerful tool in this domain, offering enhanced representational capacity and new avenues for solving symbolic discovery problems [54] [0].
A.2 Neural Networks for Matrix Preconditioning
Recent studies have explored the use of neural networks to improve matrix preconditioning tech- niques. [21 42
demonstrate the effectiveness of neural networks in refining multigrid preconditioning algorithms, thus streamlining the computa- tional process. [9] utilized Convolutional Neural Networks (CNNs) for the optimization of block Jacobi preconditioning algorithms, while (56) developed corresponding Incomplete Lower-Upper Decomposition (ILU) preconditioning algorithms leveraging machine learning insights. Although these algorithms achieved impressive results, they still face challenges such as limited interpretability and reduced computational efficiency when deployed in pure CPU environments. This paper attempts to address these issues by incorporating symbolic discovery into the framework.
57
67
65 12 41 36 37 38, 48
28
39
33
B Detailed introduction of matrix preconditioning
B.1 Overview of Matrix Preconditioning Methods
Jacobi Method:
The Jacobi preconditioner utilizes only the diagonal elements of matrix to precondition linear system. By approximating the inverse of the diagonal matrix, this method is computationally simple and effective for systems with strong diagonal dominance. However, its convergence rate can be slow, and its performance diminishes for poorly conditioned or weakly diagonally dominant matrices. The Jacobi method is typically used as a baseline for comparison with more sophisticated preconditioners 53
Gauss-Seidel (GS) Method: The Gauss-Seidel preconditioner improves upon the Jacobi method by considering both the lower triangular and diagonal parts of the matrix in sequential manner. Unlike the Jacobi method, which updates all variables simultaneously. the GS method updates each variable in sequence using the most recent values. This leads to faster convergence, especially for diagonally dominant matrices. However, the GS method can still struggle with poorly conditioned systems, and its forward-only approach can limit performance in some applications 53
Successive -ela-Relaxttion (SOR): The SOR method builds on the Gauss-Seidel method by introducing relaxation factor w to accelerate convergence. This factor allows for over- relaxation (w >> 1) or under-relaxation (w < 1), tuning the method for faster performance on certain types of problems. SOR can significantly reduce the number of iterations needed for convergence compared to both the Jacobi and GS methods, but choosing the optimal relaxation factor is
(68)
Symmetric Successive Over-Relaxatie (SSOR): SSOR a symmetric version of the SOR method, where relaxation is applied in both forward and backward sweeps of the matrix. This bidirectional process improves stability and well-suited for use with iterative solvers like the conjugate gradient method, which requires symmetric preconditioners. SSOR's s
22

---

symmetry ensures that the preconditioner maintains the properties needed for efficient and stable convergence, making a popular choice for symmetric positive-definite systems [T8].
Algebraic Multigrid (AMG): AMG is advanced preconditioning technique designed to handle large, sparse systems of linear equations, especially those arising from the discretiza- tion of partial differential equations. Unlike traditional methods, AMG operates multiple levels the matrix structure, coarsening the matrix to form a hicrarchy of smaller systems that casier to solve. Solutions the coarser grids are then interpolated back to the finer grids. This multilevel approach makes AMG highly efficient for large-scale problems, as dramatically reduce the number of iterations needed to achieve convergence. AMG is often used in combination with methods like SSOR or Gauss-Seidel as a smoother on each grid level, and is particularly effective in cases where the problem exhibits multiscale nature
Relationship Among Jacobi, GS, and SOR Methods: The Jacobi method is the simplest of the three, using only diagonal information. The GS method improves upon the Jacobi method by using both diagonal and lower triangular matrix elements achieve faster convergence. SOR further refines the GS method by introducing relaxation factor to optimize the update process. Both the GS and SOR methods can be seen iterative improvements on the Jacobi method, with SOR offering more flexible and potentially faster alternative by adjusting the relaxation factor. SSOR extends SOR symmetrically, making suitable for use in more advanced iterative solvers like the conjugate gradient method [33 18]
B.2 Parameters in Matrix Preconditioning
The choice of preconditioning parameters significantly influences the effectiveness of the precon- ditioning process, especially in the iterative solving of linear systems [7]. Below, we discuss three specific preconditioning techniques SOR,, SSOR, and key parameters affect the preconditioning results.
particularly how their
B.2.1 Relaxation Factor w in SOR and SSOR Methods
In the SOR preconditioning method, the relaxation factor is a critical parameter that determines the acceleration of iteration. SOR evolves from the Gauss-Seice method by introducing w to speed up convergence. The SOR iteration formula is given by:
[++1) (D+ WL [1 w) (M) + wb wUx(()]
(8)
where D, L, and U are the diagonal, strictly lower triangular, and strietly upper triangular parts of the matrix A, respectively [18]
The SSOR preconditioning method can be represented by the following formula:
MSSOR
D - WU)D-' wL),
w2
where MSSOR constitutes the preconditioner, and D, L, U, and defined similarly to their roles in the SOR method. This symmetrical formulation enhances the stability and effectiveness of the preconditioning, particularly benefiting symmetric positive-definite matrices by optimizing the convergence properties of the iterative solver [18]
The choice of w directly impacts the speed of convergence and the condition number of the matrix. Different problems and scenarios often require different choices w, which typically need to de- termined based on the specific properties of the problem and through numerical experimentation [18). In the PETSc library, the default relaxation factor w for both SOR and SSOR is set to 1, at which point SOR degenerates to GS preconditioning.
B.2.2 Threshold Parameters fr in AMG
In the AMG method, the threshold parameter 0r determines whether the non-zero elements of the matrix are "strong" enough to be considered in the construction of coarse grid during the multigrid
23

---

process. This parameter is crucial for establishing the connectivity between coarse and fine grids in the hierarchical multilevel structure [ST]
The AMG method solves the equation system through multiple levels of grids, each corresponding version of the original problem. During this process, the threshold parameter is used to determine whether a given non-zero matrix element is strong enough to keep the corresponding grid points connected during coarsening.
A lower threshold often leads more elements being considered strong connections, which might increase the complexity of the coarse grid but can help preserve the essential features of the original problem, thus improving the efficiency and convergence of the multigrid method.
A higher threshold might result in fewer strong connections, thereby reducing the complexity of the coarse grid. However, this weaken the effectiveness of the AMG method, especially in maintaining the features of the original problem.
Different values of Or directly influence the condition number the preconditioned matrix. Selecting the appropriate threshold parameter typically involves considering the specific structure and features of the problem, and adjustments are made through experimental fine-tuning to achieve the optimal balance [61]. In the PETSc library, the default threshold parameter 0z is set to O.
24

---

C
Algorithm Pseudocode
Algorithm 1: RNN-based Symbolic Discovery Process
Input: RNN with parameter 0, the library of tokens
parent (i + empty node
20 f- (O) sibling(() (ibblngg))))))))
// x is the concatenation of parent and sibling // Initialize hidden state of RNN
modes
ho 0
for t 1............................1 do
(4h, h4) next token the ApplyConstraint (h,, F)
1) // de is the categorical distribution of the Regularize the distribution // Arity(T,) denotes the number of operands of т
Sample token Ti ~ th
Arity(T:) > 0 then
(arent(t) f Tt sibling () + empty node
else
// When Arity(T:) =0, go back the last incomplete operator node count +- 0 for i 1,
do
// Backward iteration
count + count + Arity(T,) -1
if count = 0 then
parent(t) < Ti
sibling(t) (f) t Ti+1
break
count then
// The expression sequence is complete
break
parent(t)]] (sibling())
St
Output:
Prefix expression sequence
Algorithm 2: Deep Symbolic Optimization for Matrix Preconditioning Parameter
Input :RNN with initial parameter Oo, the library of tokens L, batch size N, iteration number
J, risk factor e, and learning rate
< 0 j + repeat
for
i 1,2, ..., N
N do (0) + SymbolicDiscover(.
f- argmax {g in T as constant placeholder: R(T:E))
// Constant
optimization
(0) f ReplaceConstant
Compute タ1 using (1) and 0
See Eq.
Compute 92 as entropy gradient
- + agğ + 92)
// Update the parameter
Train model: update pe via PPO by optimizing J(0;e)
until or convergence
Output: : The best symbolic expression T*
25

---

D Experiment Settings
D.1 Datasets
1. Darcy Flow Problem
We consider two-dimensional Darcy flows, which can be described by the following equation 35 49
27 40]:
(10) where K is the permeability field, h is the pressure, and is a source term which can be either constant or a space-dependent function.
(K(z.))Vh(z.g))
In our experiment, K (x, 3) is generated using truncated Chebyshev polynomials. We convert the Darcy flow problem into a system of linear equations using the central difference scheme of Finite Difference Methods (FDM) (32). The coefficients of the Chebyshev polynomials serve as input features for our symbolic discovery framework.
2. Second-order Elliptic Partial Differential Equation
We consider general two-dimensional second-order elliptic partial differential equations, which are frequently described by the following generic form [4
Lu 01 Ugz 0122 + 02222yy + Ur 02lly + agu I
(11)
where do, 01, 022 011, d12, 022 are constants, and f represents the source term, depending æ, y- The variables u, Ug Uy are the dependent variables and their partial derivatives. The equation is classified elliptic 4011022 > a/2-
In our experiments, 411, 022, 01, 02, are uniformly sampled within the range (-1, while the coupling term 012 sampled within (-0.01, 0.01). 0.00).. We then select equations that satisfy the elliptic condition to form our dataset. Similar to the approach with the Darcy flow problem, we convert the PDE into a system of linear equations using the central difference scheme of FDM. The coeffi- cients do, a1, a2, @11, 412, serve as input features for our symbolic discovery framework. When discussing SSOR preconditioning, we set aj and a2 zero to ensure the resulting matrix remains symmetric.
3. Biharmonic Equation
We consider the biharmonic equation, fourth-order elliptic equation, defined on a domain Ω C R.2. The equation is expressed as follows 9 17
5)
4d f in Ω = [0, a] x [0, b],
(12)
where 7
V2V2 represents the biharmonic operator and f = 4.0m sin(Tx)) sin(ry) is the prescribed source term.
In our experiments, we construct the dataset by varying the solution domain f2 [0, a] [0, bj. We utilize the discontinuous Galerkin finite element method from the FEniCS library to transform this problem into system of linear equations [5]. The parameters a, of the domain serve as input features for our symbolic discovery framework.
4. Poisson Equation
We consider two-dimensional Poisson equation, which can be described by the following equa-
tion 65 24 70]:
(13) Physical Contexts in which the Poisson Equation Appears: Electrostatics: 2. Gravitation; Fluid Dynamics.
V2u = in Ω = [o, 1J2.
In experiments, we address the Poisson equation within square domain where both the boundary conditions on all four sides and the source term f on the equation's left-hand side are generated using third-order truncated Chebyshev polynomials. The finite difference method with central difference scheme is employed to discretize the equation into a linear system. The Chebyshev coefficients serve as parameters for our symbolic discovery framework [3]
26

---

5. Thermal Problem
We consider two-dimensional thermal steady state equation, which be described by the following
equation 65 55 26
8T
8T
0y² - 0,
(14)
Oxp2
where T is the temperature. We examine the steady-state thermal equation in thermodynamies. As with the previous equation, we still solve this equation in the square domain. The boundary temperatures on the left and right boundaries are determined by random values ranging from -100 to and to 100, respectively. The top and bottom boundary temperature functions are generated by third-order truncated Chebyshev polynomials. The boundary temperature and the coefficients of the Chebyshev polynomials serve as parameters for our symbolic discovery framework.
6. Lippmann-Schwingg Equation
We consider the 2D Lippmann-Schwinge integral equation describing quantum scattering pro- cesses 63
(xr) o(r) / )
(15)
where ((r) is the incident wavefunction and Go denotes the free-particle Green's function. For symmetric potentials, we parameterize (r) using Chebyshev polynomials on radial grid r [0, Rmaz]:
N
2r 2r_
V(r) Σ Ca k=0
Rmozz
(16)
Through Galerkin discretization with Chebyshev basis functions, we convert the integral equation into a sparse linear system:
[ - Go v
(17)
where Go contains integrated Green's 's function matrix elements and V is the diagonal potential matrix constructed from Chebyshev coefficients. This formulation enables efficient computation of scattering cross-sections across different potential configurations. The Chebyshev coefficients {cx. serve parameters for our symbolic discovery framework.
Markov chain Cootzmmann-Distributee Brownian Motion)
We consider a Markov chain optimization problem governed by the steady-state ID Fokker-Planck equation with Meebyshev-parameteriie potential energy field 45 64
(DVp+ + pV =
(18)
where (xx) = E (xT is the potential energy defined by Chebyshev polynomials Ts. D is the diffusion coefficient, and f is source term. The equilibrium distribution follows Peg e (=
In our experiments, the potential field (xx is generated using truncated Chebyshev polynomials with coefficients {cr} sampled uniformly from [-1, 1]. We discretize the Pokker-Planck equation into via spectral collocation methods with Cheby basis functions. The Chebyshev coefficients {cx+ serve as parameters for our symbolic discovery framework.
sparse, nonsymmetric linear system A((c))
D.2 Environment
To ensure consistency in our evaluations, all comparative experiments were conducted under uniform computing environments. Specifically, the environments used detailed follows:
1. Environment (Env1):
Platform: Windowsl1 version 22631.4169, WSL
Operating System: Ubuntu 22.04.3
CPU Processor: AMD Ryzen 9 5900HX with Radeon Graphics CPU, clocked at 3.30GHz
Environment (Env2):
Platform & Operating System: Ubuntu 18.04.4 LTS
27

---

CPU Processor: Intel(R) Xeon(R) Gold 6246R CPU at 3.40GHz
GPU Processor: GeForce RTX 3090 24GB Library: CUDA Version 11.3
Speed tests for solving linear systems were performed in Env 1, while all training related to symbolic discovery was conducted in Env 2.
D.3 Tolerance Metric for Linear Systems
Relative Residual:
To estimate the relative accuracy of the solution i for a linear system Ar = b, we employ the relative residual follows:
||A bl12
Relative Residual
(19)
|
Here, å represents the solution predicted by the model, and is the right-hand side vector of the linear system. When ë is the exact solution, the Relative Residual equals O.
D.4 Training Data Generation
We employed an adaptive grid search to generate the training dataset. Initially, we traversed coarse grid, sampling every 0.05, and from this dataset, we selected the three points with the smallest values. Subsequently. we conducted finer grid search around these points, sampling every 0.001, identify the point with the minimum value, which we designated as our optimal parameter. Particularly, after experimental validation confirmed the dataset's convexity, we utilized binary search sampling method for a dataset derived from the second-order elliptic equation's 's SOR preconditioning, Starting with points at 0.0. 1.0, and 2.0, we compared these values. If the value at 0.0 was lowest, we computed at 0.5; if at 2.0, then at 1.5; and if at 1.0, then at both 0.5 and 1.5. This process was repeated until minimum point was achieved with a precision of 0.001.
For SOR preconditioning, we evaluated second-order elliptic equations, Darcy flow equations, and biharmonic equations, with solution time as the metric for optimal preprocessing parameters, achieved by minimizing solution time using the previously described grid method. In SSOR preconditioning. applied to second-order elliptic and Darcy flow equations. we utilized hybrid metric that combined normalized computation time and iteration counts, aiming to simultancously optimize both iteration counts and solution times. For AMG preconditioning, also examined with second-order elliptic and Darcy flow equations, we used the condition number of the preconditioned matrix as the metric, where lower value indicates better performance.
D.5 Parameters of SymMAP
Experimental Setup. SymMAP is implemented using the LSTM architecture with one layer and 32 units. More details about the hyperparameters provided in Table 7
Table 7: Hyperparameters of SymMAP (Default Model)
Hyperparameter
Value
Number of LSTT layers
Number LSTM units
32
Number of training samples
2,000,000
Batch size
1,000
Risk factor
0.05
Minimal expression length
Maximal expression length
64
Learning rate
0.0005
Weight of entropy regularization
0.03
Restricting searching space. We employ specific constraints within our framework streamline the exploration of expression spaces effectively and ensure they remain within practical and manageable bounds:
28

---

1. Bounds on expression length. To strike a balance between complexity and manageability. we set boundaries for expression lengths: a minimum of 4 and maximum of 64 characters. This ensures that expressions are neither overly trivial nor excessively complicated.
2. Constant combination. We restrict expressions such that the operands of any binary operator are not both constants. This is out of the simple intuition that, if both operands are constants, the combination of the two can be precomputed and replaced with a single constant.
3. Inverse operator exclusion. We preclude unary operators from having their inverses as children to avoid redundant computations and meaningless expressions, such as in log(exp(2:)).
4. Trigonometric Constraints. Expressions involving trigonometric operators should not include descendants within their formulation. For instance, sin + cos(z)) is restricted because it combines trigonometric operators in a way that is uncommon in scientific contexts.
D.6 Computational Time for Related Algorithms
Dataset Generation Time:
Darcy Flow Problem: 40 hours
Second-order Elliptic Partial Differential Equation: 40 hours
Biharmonic Equation: 100 hours
Poisson Equation: hours
Thermal Problem: 6 hours
SymMAP Training Time: For each run, 1000 iterations are performed.
Without polynomials the Token Library: approximately 800 seconds.
With polynomials the Token Library: approximately 2600 seconds.
MLP Training Time (in Experiment 5.2)::
Training until full convergence: 1.5 hours.
29

---

E Experimental Data and Supplementary Experiments
E.1 Symbolic Expressions from Main Experiments
This section documents some of the learned expressions from the main experiments, corresponding to "SymMaPI" in Tables and 3
Second-orde elliptic PDE problem, AMG preconditioning:
21 13 + 1.0
+ 7.0
Parameter meanings: 21-16 represent the coefficients 611, 012 a22, 01, 02, and ag in the second-order elliptic equation.
Second-order elliptic PDE problem, SSOR preconditioning:
(2z- 4z 2zg
2- -.2356677.244
0-0.1178522
+ 0.22994224 -
- 44g
-7.534392; + 61
0.389893 0.3144x:4 0.0111835)
Parameter meanings: 21-24 represent the coefficients d11, 012, 022, and ag in the second- order elliptic equation.
Darcy flow problem, SSOR preconditioning:
1.0
2 (I14 .4) + 1.0
Parameter meanings: T1-Ti6 represent the 16 coefficients of second-order truncated
Chebyshev polynomial in two dimensions, ordered as follows: 1, z, x², x², x²y, x²y y², xy? x² y x² y², y zyy 2zy², and z²y².
Darcy flow problem, AMG preconditioning
1.0
1.0 +
1.0
210+23+05+2.0
Parameter meanings: 116 represent the 16 coefficients as described above.
Biharmonic Equation, SOR preconditioning:
1.0ga
1.0
4x²
1.0
10
22
Parameter meanings: and represent the length and width of the equation's boundary, respectively.
Biharmonic Equation, AMG preconditioning:
1.0
1.0
3.0 + 1.z. + 1.0 32
x2
Parameter meanings: 21 and 22 represent the length and width of the equation's 's boundary, respectively.
30

---

Poisson Equation, SOR preconditioning:
1.0
cxp
33 + exp(2 exp(x))
Parameter meanings: 21-28 represent the coefficients of two second-order truncated Cheby- shev polynomials for the boundary functions.
Poisson Equation, SSOR preconditioning:
1.15024107160485 (0. (0.106506997999920..; 1)
1/16
Parameter meanings: I1 -ig represent the coefficients of two second-order truncated Cheby- shev polynomials for the boundary functions.
Poisson Equation, AMG preconditioning:
1.0
z + 7.0
Parameter meanings: 2g represent the coefficients of two second-order truncated Cheby- shev polynomials for the boundary functions.
Thermal problem, SOR preconditioning:
0. 0.778800783071405 1/4
exp
(1-x6)
Parameter meanings: -1-27 represent the coefficients of the Chebyshev polynomial for the boundary temperature function on the upper and lower boundaries, while .rş and represent the coefficients for the boundary temperature function on the left and right boundaries.
Thermal problem, SSOR preconditioning:
1.0
1.0
log (4.0(1 - 0.5zz)²)
Parameter meanings: 41-24 represent the coefficients the Chebyshev polynomial for
boundary temperature function on the upper and lower boundaries, while Jzş and represent the coefficients for the boundary temperature function on the left and right boundaries.
Thermal problem, AMG preconditioning:
1.0
2.71 0.718281822 2.71828182845900 (3355555000000001133333333333333 + 8.15484548837714
Parameter meanings: I -84 represent the coefficients of the Chebyshev polynomial for the boundary temperature function on the upper and lower boundaries, while ay and ze represent the coefficients for the boundary temperature function on the left and right boundaries.
31

---

E.2
Experiments
To further validate the ability of SymMaP to predict multiple parameters, we conducted experiments on AMG preconditioning with an SOR smoother, optimizing two parameters simultaneously: AMG's threshold parameter (default PETSc: 0) and the SOR relaxation factor w.
All other settings match the AMG preconditioning described in the main experiment. We jointly optimize both parameters, using the condition number as the evaluation metric.
Table 8: Condition number comparison across different PDE problems and parameter settings.
Problem
Optnall
SymMaP
0:2
1,3
Coustant
Duncy Flow
10,000
752862
8012
7623
8321
21464
17326
18354
12496
9341
9852
6432
4021
Elligic PDE
40,000
6792
108.1
161.2
163.8
191.5
183.2
187.3
213.4
194.2
198.4
159.3
156.1
Poisson
10,000
1242
4.31
4.71
3.56
67.43
70.34
54.26
67.43
70.34
54.26
3.45
2.94
Thermal
2,800
7525
11.2
10.3
12.6
463.2
432.1
48932
463.2
432.1
489.2
8.64
7.53
Bolizmann
50,000
9548
246.4
195.7
197.3
263.9
213.4
219.3
353.5
246.3 3
283.9
189.2
165.3
Lippoaan--Se Sctwinggr
50,000
4832
4937
12842
6383
3981
4372
2395
This ously.
SymMaP's ability to handle multiple preconditioner parameters simultane-
E.3 Custom Metric Experiments
In the main paper, experimental focus was primarily SOR, SSOR, and AMG preconditioners. We selected these methods because they exhibit a non-monotonic relationship between their key parameters (e.g.. the relaxation factor the threshold parameter 0r) and performance metrics such as solution time or condition number. This characteristic creates a non-trivial optimization landscape with distinct optimal parameter, making these preconditioners ideal for demonstrating the core capabilities of our symbolic discovery framework, SymMaP.
We acknowledge that other important preconditioners, such as Incomplete Cholesky (ICC) and Incomplete LU (ILU) factorization, present a different type of optimization challenge. For these methods, parameters like the fill-in level typically have monotonic relationship with individual performance metrics. For instance, increasing the fill-in level generally decreases the condition number of the preconditioned matrix but monotonically increases the computational cost of the factorization. A practical application, therefore, requires balancing this trade-off by defining custom, often heuristic, objective function. To maintain clarity and avoid introducing ambiguity from such custom metrics in our main experiments, we initially prioritized the more straightforward optimization cases offered by SOR, SSOR, and AMG.
To demonstrate the broader applicability and generalizability of SymMaP, we conducted an additional experiment on ICC preconditioning. In this scenario, the goal is to find an optimal fill-in level that balances the trade-off between matrix conditioning and computational cost. We defined a hybrid objective function to capture this balance:
Objective = 0.03 Condition Number + Time(s)
(20)
The experiment was performed on the symmetric second-order elliptic PDE problem (with matrix size of 40,000), maintaining consistency with the experimental settings described in the main paper. The task for SymMaP was to discover a symbolic expression for the optimal fill-in level that minimizes the objective function defined in Equation 20
Table 9: Performance comparison for ICC preconditioning on the second-order elliptic PDE dataset. The objective value calculated using Equation 20 where lower value is better.
Method
None
Default
Optimal Constant SymMaP
Objective Value
227.66
24.01
22.74
22.27
The results, presented in Table 9 show that SymMaP successfully identifies parameterization that outperforms the default and optimal constant strategies, achieving the lowest objective value. This experiment confirms that SymMaP is not limited to preconditioners with non-monotonie optimization
32

---

landscapes but can also effectively handle trade-off-brase off--bsedd optimization problems by learning symbolic policies for custom, obbti-obeective metries. This demonstrates the framework's 's flexibility and generalizability to wider range of preconditioning scenarios.
E,4
Generalizability Analysis
E,4.1
Generalization Across PDE Parameters
If the matrix properties remain stable under parameter variation, the learned symbolic expressions can generalize beyond the training range. Otherwise, as noted, performance may degrade. To validate this, conducted additional experiments:
Second-order elliptic equation, grid size: 40,000, test range: a E [-2, coupling term [-0.5, 0.5] (vs. training range: E [-1, 1], coupling term E [-0.01, 0.01)), tolerance: 10- preconditioner: SOR, linear solver: GMRES.
Table 10: Comparison of average computation times (seconds) for SOR with different selections, and tolerance Ic-3,
None
PETSc default
Fixed constant 0.1
Fixed constant 1.9
Optimal constant
SymMaP
Time (s)
16.98
4.04
1.26
15.2
0.94
0.86
The results in Table T0 demonstrate that the expressions learned by SymMaP perform well even side the range of training parameters, indicating that SymMaP has strong generalization capabilities.
E.4.2 Generalization Across Matrix Size
In preliminary tests, we observed that for certain preconditioners (e.g., SOR), the optimal precondi- tioning parameters are largely independent of the matrix size (resolution) for given problem. To ensure consistency, we fixed the matrix size in our experiments.
To validate the reasonableness of this observation, conducted additional experiments. Specifically, we trained a model on a dataset composed of matrices with mixed side lengths (uniformly x 10 to 5 x 104) and tested it test with similarly mixed side lengths.
11.03
13.00
6.4%
Figure 3: SOR optimal parameter distribution for mixed edge-length second-order elliptic PDE dataset.
Table 11: SOR experiment mixed-size datasets (uniformly 10 to 5 10*). Other settings consistent with the main experiment. Average time is used as the metric (in seconds).
PETSc default
Fixed constant w
Dataset
None
Optimal constant
SymMaP
0.1
0.2
1.8
1.9
Biharmonic
745
258
542
512
231
242
193
175
Darcy Flow
164
65.3
102
84.1
51.7
61.8
47.4
39,5
Elliptic PDE
27.2
18.6
19.2
20.8
15.2
16.7
12.4
10.8
Poisson
16.7
10.3
11.5
10.8
10.2
12.9
8.53
8.12
Thermal
87.3
26.4
83.4
80.2
46.2
48.5
24.5
22.3
33

---

The distribution of the optimal parameters across the second-ores order elliptic PDE dataset is shown in Figure The results in Table IT confirm that the matrix size has minimal impact on the optimal parameters for SOR. Furthermore, for preconditioners where the matrix size does influence the parameters, SymMaP can seamlessly incorporate the size as an additional input feature.
E.4.3 Generalization Across Variable Geometries
We conducted experiment where the problem's complexity arises from varying geometric domain rather than from PDE coefficients. We chose a second-order elliptic PDE defined on a variable triangular domain. For each problem instance, the triangle's 's three vertex coordinates were randomly sampled. These six coordinates were then provided as input features to SymMaP, which was tasked with predicting the optimal relaxation factor, for the SOR preconditioner. The performance was measured by the total solution time.
Table 12: Comparison of average computation times (in seconds) for SOR preconditioning on second-order elliptic PDE with a variable triangular domain.
Method
None
Default (w 1)
w 0.2 w 1.8
Optimal Constant
Theoretical Optimal
SymMaP
Time (s)
16.7
8.3
9.04 9.01
7.98
7.65
7.72
The results, shown in Table 12 indicate that SymMaP performs robustly, achieving solution time very close to that of the empirically determined optimal constant. This experiment validates that SymMaP can effectively learn policies for problems with geometric variations by parameterizing the domain shape and using those parameters as input features.
E,5 Error Reduction Trajectory
We evaluated the performance of different parameter settings the biharmonic equation dataset under varying tolerance settings. Specifically, we measured the solution time and plotted the error reduction curves. Other settings are consistent with the main experiment. As shown in Figure the experimental results demonstrate that SymMaP consistently delivers strong performance across all accuracy levels.
IL'S drficl PAAAEE
Lntinal conerr
Thme
Figure 4: Error reduction curves for different parameter settings on the biharmonic equation dataset.
34

---

E.6 Statistical Insights from the Second-Order Ellipticc Equation Experiments
Here, we use the second-order elliptic equation as an example present more details of the statistical characteristics of experimental data, as shown in Table
Table 13: Performance comparison for solving second-ord elliptic equation with an edge length of
40,000 and a tolerance of 10-7 using SOR preconditioning. The coefficients range from (-0.1, 0.1) for the main parameters and (0.005, 0.01) for the coupling terms. Average time is used as the metric (in seconds).
Statistic
None
PETSc default
Fixed constant 0.1
Fixed constant 1.9
Fixed constant
SymMaP
Mean 48.10
33.02
36.54
28.09
26.28
22.67
Median
11.79
9.87
9.99
8.53
8.42
8.04
Upper quartile
63.67
42.81
49.54
39.78
35.79
29.67
Lower quartile
3.38
2.68
2.65
2.31
2.10
2.10
Variance
1365.1
822.47
864.45
784.92
749.35
705.27
E.7
Comparative Experiments with Learning-Basel Approaches
SymMaP enhances traditional preconditioning algorithms by discovering symbolic expressions for optimal parameters. In our experimental design, we deliberately avoid direct comparisons with other learning-based approaches for the following reasons:
1. Focus on Preconditioning Enhancement: SymMaP is specifically designed to improve traditional preconditioners without altering their core algorithms. fair comparison would require other learning-based methods that also preserve the original preconditioning process. However, no such methods currently exist.
2. CPU Compatibility; Linear solver environments are predominantly CPU-based. SymMaP's symbolic expressions integrate seamlessly into CPU workflows, whereas neural network- based approaches often lack efficient CPU deployment.
3. Generality: SymMaP is highly adaptable diverse preconditioners and linear solvers. In contrast, existing learning-basee methods either: (a) combine with specific preconditioners for ICC-preconditioned CG;
(e.g., (11), (b) target narrow cases (e.g., 60 2
exclusively for CG).
4. Interpretability & Safety: Numerical algorithms demand rigorous analysis. Opaque predic- tors, such as neural networks, risk violating theoretical constraints (e.g., avoiding w 0 or in SOR), SymMaP's symbolic expressions enable proactive avoidance of such issues through analytical guarantees.
While these advantages are significant, we also conducted comparative experiments with other learning-basse preconditioning algorithms evaluate their performance. The experiments are divided into two groups: 1. Table [4 focuses symmetric matrices, using the CG algorithm with tolerance of 10 2. Table 15 targets nonsymmetric matrices, employing the GMRES algorithm with tolerance of 10" 10-7
All other settings are consistent with those described in the paper. The performance metric the average solve time, measured in seconds, where lower values indicate better performance.
35

---

B Syl-SSSSHSOOO
MD SD--SyMMPPPP.
0.:3
0.44
25
1.23
3.23
52
Table 15: Comparison of average computation times (seconds) with different methods, linear sol ver GMRES, and tolerance is Ie-7.
Table 14: Comparison average computation times (seconds) with different methods, linear solver is CG, and tolerance is 1c-7.
0.2
6.3
53
1.00
1.05
72
74
641
1.30
1.26
302
302
23
0.3
9.1
36.7
(00
0,43
0.89
23
30
90
272
33
032
83
254
x
0.33
+
0.52
B
0:42
94
22.
6.47
0:01
143
412
172
173
0.32
158
333
4:34
244
241
192
195
0.1
31
315
CD
1:94
5.99
133
60) -000
0.:11
0.33
PM
102
BDD
2.00
2:14
11.5
0:2
4:31
4.35
182
155
2011
201
132
0.45
645
14
JMGG
AMGH50R
155
155
182
AMG
1423
157
362
4.21
4:1
TKU
120
6.80
103
2.90
200
13A
124
636
200
Oaliuraon
(P0 & DPC
180 CPU
CPU
780
321
CTU GIU
Obee CPU
DPU GPU
Oet CHU
Stxt
630
4000
Elosk FIE
Daadd
36

---

E.7.1 Notes on These Methods
Classification of Existing Algorithms:
1. Methods Enhancing Linear Solvers: These approaches add optimization modules or modify the iteration process (e.g., [4T. 21 These methods target different components than SymMaP and could potentially be combined with our approach for further improvements.
2. Methods Predicting Preconditioning Matrices: These approaches use neural networks to predict preconditioning matrices (e.g.. [60)
34
Code Availability:
34
The official implementation is available on GitHub (NeuralPCG).
23
The official implementation is available on GitHub (NeuralIF).
60
The official implementation is available on OpenReview (PreCorrector).
The official implementation is available GitHub (NeurKItt).
21
The official implementation is available on GitHub (Learning to Optimize Multigrid
Solvers).
[24]
The official implementation is available GitHub (Leamnngg Neural PDE Solvers
with Convergence Guarantees).
Table details explanation:
Combinations: Sequential application of methods (e.g., 'AMG+?' or "141++?).
x: Method incompatible (c.g.. [34], 22 [24]. [60] only work for symmetric matrices/CG)
excessively slow (e.g.. [24 100 for Elliptic PDE).
None: No preconditioning.
SOR/SSOR/AMG///C/I... PETSC defaults.
SymMaP(SOR/SSOR) Optimized via SymMaP.
41: Accelerates Krylov solvers via recycling (distinet from preconditioning): combinable with other methods.
[2I): An AMG variant: allows nested preconditioning but conflicts with learning- based methods (B44, [23], 660).
244: A fixed-point iteration (non-Krylov); computationally prohibitive for large sparse systems.
Experimental Setup:
Hardware:
-PUU: AMD Ryzen 9 5900HX, MPI 8-core parallel.
-PUU: RTX 3090 (CPU-only uses PyTorch CPU inference).
Averaging: 1000 runs per problem.
Linear Solver: For faimess,, use PETSc's CG/GMRES.
E.7.2 Results
As shown in Tables 14 and 15 across all experiments, the combination of "y4Ma+ymMaP' consis- tently achieved the shortest computation times. When evaluating standalone algorithms (without combinations). SymMaP alone demonstrated the fastest performance in all test cases. The perfor- mance advantage of SymMaP particularly pronounced in CPU-only environments. SymMaP generalizes to any parameterized preconditioner, unlike other methods (e.g.. (23) for ICC, [21] for AMG). These results conclusively demonstrate the superior performance characteristies of the SymMaP approach.
37

---

E,8
Interpretable Analysis Details
Table 16: Symbolic expressions learned from the main experiments
Precondition
Dataset
Symbolic expression
SOR
Biharmonic
1.0 + +
1.0/21
SOR
Elliptic PDE
1.0 + 1.0/(22 + 1.0 + 1.0(112 + +.0))
SOR
Darey Flow
1.0 1.0/(zz + 1.0)
SSOR
Elliptic PDE
1.0 +.0 10/(220 + 1.2)
AMG
Elliptic PDE
(I123 + +1/7
As shown in Table 16 the variables are defined as follows: in the first row, and represent the size of the boundary for PDE solutions; in the second row, 22 represents the coefficient of second-order coupling term; in the third row, is the coefficient of the fourth x-term multiplied by the first y-ferm in two-dimensional Chebyshev polynomial; in the fourth row, 22 again denotes the coefficient of second-on -onder coupling term; in the fifth row, 21 signifies the coefficient of second-order non-coupling term.
E.9 Analysis of Hyperparameters
The performance of SymMaP is primarily influenced by the learning rate of the RNN, batch size, and dataset size. We conducted experiments to study the impact of these hyperparameters.
Symbolic Learning RNN Parameters:
Table 17: Performance comparison of SymMaP under various symbolic learning RNN parameters (lower condition numbers are preferable). The experiment focuses on optimizing AMG precondition- ing coefficients in the Darcy Flow dataset.
Learning Rate
Batch Size
Condition number
Training time(s)
500
6780
1173.09
0.01
1000
5168
863.51
2000
6898
522.80
500
5935
110..66
0.001
1000
11774
676.40
2000
5935
505.85
500
4718
1026.45
0.0005
1000
5935
703.17
2000
5935
549.36
500
12228
1324.00
0.0001
1000
7508
837.18
2000
6884
603.62
Results in Table 17 indicate that an appropriate combination of RNN learning rate and batch size can enhance performance.
38

---

Dataset size:
Table 18: Performance comparison of SymMaP across varying dataset sizes (lower condition numbers indicate better performance). The experiment evaluates the optimization of AMG preconditioning coefficients for the Darcy Flow dataset.
Dataset size
Condition number
Training time (s)
10
7032
669.68
50
6980
737.80
100
4892
812.02
500
3811
699.54
1000
5345
703.17
Table T8 demonstrates that increasing the dataset size enhances the performance of symbolic expres- sions learned by SymMaP, as expected.
E.10 Analysis of Data Generation Method
The SymMaP framework is modular, allowing the data generation step to be performed by any suitable search algorithm. We conducted this comparison on the second-order elliptic PDE problem with the SOR preconditioner (matrix size 40,000, 1,000 instances), evaluating standard grid search, random search, and Bayesian optimization against our adaptive approach. Table 19 summarizes the results, showing each method's data generation time, the Mean Absolute Error (MAF) of the discovered optimal parameters against our high-precision reference, and the final solution time achieved by SymMaP when trained the resulting dataset.
Table 19: Comparison of different data generation methods for finding the optimal SOR parameter.
Data Generation Method
Avg. Generation Time
MAE (vs. SymMaP)
Final Avg. Solution Time (s)
Standard Grid Search
100hh
7.80 10-4
16.0
Random Search (2000 samples)
~ 100hh
1.74 x 10-3
16.0
Bayesian Optimization
~ 1h
9.80 x 10-4
15.9
Adaptive Search (SymMaP)
25h
0
15.9
The analysis reveals two key findings. First, while methods like Bayesian optimization are signif- icantly more time-efficient for data generation, all tested algorithms converge to nearly identical optimal parameters (MAE < 1.74 x 10 10-3). Second, because the resulting training data consistent, SymMaP discovers the same symbolic expression and achieves the same final performance regardless of the data source. This demonstrates that SymMaP's effectiveness robust and independent the specific search algorithm used, provided it accurately identifies the optimal values.
E.11 Comparison with Alternative Interpretable and Symbolic Methods
To clarify the advantages of SymMaP, conducted additional experiments comparing it against other lightweight, interpretable, and symbolic modeling approaches. We categorize these alternatives into two groups: (1) standard interpretable models like linear and polynomial regression, and (2) alternative symbolic regression paradigms, such as Genetic Programming (GP), that do not use an RNN-based search.
E.11.1
Comparison with Standard Interpretable Models
key advantage of symbolic regression over fixed-form models its flexibility in both feature selection and functional form. While models like linear or polynomial regression are interpretable, their expressive power limited. contrast, SymMaP explores a much richer space of mathematical expressions. In fact, by constraining the library of operators (e.g., to only '+' and ***), our framework can effectively subsume these simpler models, which represent special cases within the broader search space.
39

---

To provide direct performance comparison, we conducted an experiment on the second-order elliptic PDE problem (using the SOR preconditioner, with matrix size of 40,000 and 1,000 data points). We tasked I Linear Regression, Quadratic Polynomial Regression, and SymMaP with predicting the optimal SOR relaxation factor.
Table 20: Performance comparison against standard interpretable models.
Fitting Method
Solution Time
Linear Regression
20.9
Quadratic Polynomial Regression SymMaP
17.0
15.8
As shown in Table 20) SymMaP discovers a policy that yields significantly faster solution times. This demonstrates that its expressive power is crucial for accurately capturing the complex, non-lineat landscape of optimal preconditioning parameters.
E.11.2 Comparison with Genetic Programming- Based Symbolic Regression
A more direct comparison is against other symbolic regression paradigms. Our choice of an RNN+RL approach was primarily motivated by its superior scalability to the high-dimensional feature spaces often encountered in preconditioning problems. For instance, the Darcy Flow dataset in our paper has 16 input features. Genetic Programming (GP) algorithms, which form the basis of popular tools like PySR, are known to scale poorly with a large number of features. PySR itself warns users that its performance can degrade high-dimensional settings, recommending dimensionality reduction as prerequisite.
To quantify this, we conducted head-to-head comparison between SymMaP and PySR on finding the optimal SOR relaxation factor for the Darcy Flow and Biharmonic datasets. For PySR, we report the result from the perfforrroninn expression found across five runs.
Table 21: Performance comparison of SymMaP against PySR.
Dataset
None
Default
Optimal Constant
PySR
SymMaP
Darcy Flow
33.1
13.5
9.54
8.87
8.50
Biharmonic
7.67
2.04
1.31
1.41
1.24
The results in Table 21 21 show that while both methods find effective expressions, SymMaP consistently outperforms PySR. Notably, for the Biharmonic problem, the expression found by PySR failed im- prove upon the optimal fixed constant, likely due to the performance limitations of GP this context. Furthermore, the RNN+RL approach significantly more efficient, requiring approximately 15 minutes to discover its expressions, whereas PySR took over hour.
In summary, for the specific challenge of preconditioning parameter prediction, which often involves larger number of input features, the RNN+RL search strategy in SymMaP offers superior performance and greater efficiency compared to both simpler interpretable models and alternative GP-based symbolic regression tools.
40


Reviews:
[{"review_id": "dhDMEVUz4X", "reviewer": "Reviewer_BuLC", "summary": "This paper presents a symbolic-based machine learning method for identifying optimal parameters for matrix preconditioning. The method is described in detail and offered as an end-to-end solution. Experimental results demonstrate that it can identify parameter values close to the optimum, thereby enhancing the adaptability of linear system solvers across different tasks.", "questions": "Please find my questions in the \"Strengths and Weaknesses\" section. Thanks.", "limitations": "The authors briefly discuss the limitations of their approach, but there are several important aspects that warrant further attention—such as the manual effort currently required to integrate the symbolic method's output into solvers, and the limited range of supported solvers. Expanding the limitations section to address these points in more depth would significantly strengthen the paper.", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n\n1. The paper addresses a fundamental problem in optimization and machine learning. While the topic may not be currently trending, it holds significant value in the domains of High-Performance Computing (HPC) and AI for Science.\n\n2. The proposed symbolic approach offers a degree of explainability compared to traditional black-box ML methods for searching optimal parameters in matrix preconditioning. The method is well implemented, featuring a complete pipeline: data generation, parameter search using reinforcement learning (RL), and a thoughtful state representation via RNN. The output is integrated into a linear system solver to demonstrate practical benefits.\n\n3. The paper includes sufficient background material, making it accessible even to non-expert readers. The writing is clear and approachable, which is commendable.\n\n\nWeaknesses:\n\n1. I am somewhat concerned about the limited description of how the symbolic method integrates into real-world linear system solvers. The paper currently lacks detail in this area. It would be highly beneficial to include a more thorough explanation—ideally, demonstrating integration with at least two widely used solvers. This would help the audience better understand the practical applicability and impact of the proposed method.\n\n2. The evaluation scope is relatively narrow, a limitation that the authors acknowledge in the conclusion. Given that the symbolic method should theoretically generalize better to different matrix preconditioning setups—due to its reduced reliance on training data—it's disappointing that this potential is only discussed rather than thoroughly evaluated. A more comprehensive evaluation on diverse preconditioning scenarios would significantly enhance the paper's contribution, especially in terms of its claimed generalizability and explainability.\n\n3. Several design choices, such as the use of adaptive grid search and RNNs, are mentioned but not sufficiently justified or compared against plausible alternatives. It would strengthen the paper to include discussion and experimental results demonstrating why these choices are effective or superior to other options.\n\n4. Did you consider including conventional hyperparameter search methods (e.g., Bayesian optimization, random search, or grid search) in your evaluation? These methods are common baselines in this area and should be included to contextualize the performance of your symbolic-RL approach.", "was_revised": true, "final_justification": "I remain positive about this paper after the response period.", "rebuttal": "**Dear Reviewer,**\n\nWe are sincerely grateful for your thorough review and insightful feedback on our paper. Your comments are invaluable, and we appreciate the opportunity to clarify our contributions and address your concerns. Below, we address each of your points in detail.\n\n----\n\n### **On Practical Integration with Real-World Linear Solvers  (Weakness 1)**\n\nThe process of integrating our algorithm is a very simple, two-step procedure:\n\n1. Parameter Computation via the Symbolic Expression:\n\nOnce SymMaP discovers a symbolic formula, applying it to a new problem is computationally trivial. The problem's specific features are simply substituted into the learned expression to calculate the recommended preconditioner parameter. This can be summarized with the following pseudo-code: (Python)\n\n```Plain\ndef calculate_preconditioner_param(problem_features, learned_expression):\"\"\"\n    Calculates the optimal preconditioner parameter using the learned symbolic expression.\n    \n    Args:\n        problem_features (dict): A dictionary of problem features, e.g., {'x1': 0.5, 'x2': 1.2}.\n        learned_expression (str): The formula discovered by SymMaP, e.g., \"1.8\" or \"0.5 + 0.9 * x1\".\n\n    Returns:\n        float: The computed preconditioner parameter.\n    \"\"\"\n    param_value = safe_evaluate(learned_expression, variables=problem_features)\n    return param_value\n```\n\n2. Seamless Integration with Standard Solver Libraries:\n\nThe computed parameter can be directly passed to mainstream solver libraries with no modification to the solvers themselves. We demonstrate this with two widely-used examples, as you suggested:\n\n- **Integration with PETSc (via command-line):** PETSc is a standard toolkit in HPC. A user can pass the computed relaxation factor `omega` directly to a PETSc-based application using the `-pc_sor_omega` flag: (Bash)\n\n```Plain\n# The value '1.8' is calculated by our method beforehand\n./my_petsc_solver -ksp_type gmres -pc_type sor -pc_sor_omega 1.8\n```\n\n- **Integration with PyAMG (via Python** **API****):** PyAMG is a popular Python library for Algebraic Multigrid methods. The computed strength threshold `theta` can be set directly in the solver configuration: (Python)\n\n```Plain\nimport pyamg\n# The value '0.25' is calculated by our method\noptimal_theta = 0.25 \n# Pass the parameter directly into the PyAMG solver configuration\nstrength_config = ('symmetric', {'theta': optimal_theta})\nml = pyamg.smoothed_aggregation_solver(A, strength=strength_config)\n```\n\nThe computational overhead of this integration is limited to a single evaluation of the symbolic expression, which is negligible compared to the cost of the iterative linear solve. The complete code for our experiments is already available in the supplementary material.\n\n**Commitment:** We agree that this is an important practical detail. In the final version of the paper, we will add a dedicated section to thoroughly describe this integration process. Furthermore, we will package our code into a user-friendly C/Python function to facilitate its adoption and potential integration into official preconditioning libraries.\n\n### **On Broadening the Evaluation Scope and Generalizability (Weakness 2)**\n\nOur initial focus was on SOR, SSOR, and AMG because they exhibited a non-monotonic relationship between their parameters and performance metrics (e.g., condition number or solution time). This presents a non-trivial optimization problem where a clear optimal parameter exists, making them ideal candidates for showcasing SymMaP's core capability.\n\nWe acknowledge that other important preconditioners, such as ICC and ILU, present a different challenge. For these, increasing a parameter like the fill-in level typically decreases the condition number but monotonically increases the computational cost of the factorization. A real-world application requires balancing this trade-off, which involves defining a custom, often heuristic, objective function. To avoid introducing ambiguity with such custom metrics in our initial submission, we prioritized the more clear-cut optimization cases.\n\nTo further address your concern and demonstrate the broader applicability of SymMaP, we have conducted a new experiment on **ICC** **preconditioning**. We defined a hybrid metric: `Objective = 0.03 * Condition_Number + Time(s)` to find the optimal fill-in level. The experiment was run on the symmetric second-order elliptic PDE (matrix size 40,000), with other settings consistent with paper.\n\n| Method          | None   | Default | Optimal Constant | SymMaP |\n| --------------- | ------ | ------- | ---------------- | ------ |\n| Objective Value | 227.66 | 24.01   | 22.74            | 22.27  |\n\nThese results confirm that SymMaP can effectively optimize parameters for different preconditioners, even under custom, multi-objective metrics.\n\n**Commitment:** In the final version, we will incorporate the analysis above and expand our experiments to cover more preconditioners (including ILU) with custom objective functions. We will also release the corresponding code to demonstrate SymMaP's flexibility. We believe this will substantially strengthen the paper's claims on generalizability.\n\n### **On Justifying Design Choices and Comparing with Baselines (Weakness 3 & 4)**\n\nOur primary focus was on demonstrating the efficacy of the symbolic learning framework itself, for which our chosen methods proved effective. The use of an RNN is a natural choice for processing sequential prefix notation data. However, the reviewer raises an excellent point that the data generation step could be performed by other established search methods.\n\nTo address this, we conducted a comparative experiment on the data generation step for the 2nd-order elliptic PDE with SOR (matrix size 40,000, dataset size 1,000). We compared our adaptive grid search against standard grid search, random search, and Bayesian optimization.\n\n| Data Generation Method       | Avg. Generation Time | Mean Absolute Error (vs. our method) | Final Avg. Solution Time (s) |\n| ---------------------------- | -------------------- | ------------------------------------ | ---------------------------- |\n| Standard Grid Search         | ~100h                | 7.80E-04                             | 16                           |\n| Random Search (2000 samples) | ~100h                | 1.74E-03                             | 16                           |\n| Bayesian Optimization        | ~11h                 | 9.80E-04                             | 15.9                         |\n| Our Adaptive Search          | ~25h                 | 0                                    | 15.9                         |\n\nThis analysis leads to three key conclusions:\n\n1. **Consistent Optimal Parameters:** All search methods, when run to sufficient precision, converge to nearly identical optimal parameters (MAE < 1.74e-3).\n2. **Consistent Symbolic Discovery:** Since the training datasets are virtually identical, SymMaP discovers the same symbolic expressions (differing only by minute constant values), resulting in nearly identical final performance (a 0.1s difference is within the margin of error).\n3. **Efficiency of Alternatives:** As you astutely pointed out, other methods like Bayesian optimization are significantly more efficient for generating the training data.\n\nThis is a valuable insight.\n\n**Commitment:** In the final version, we will update the paper to include this comparative analysis, explicitly justifying our design choices and contextualizing them with these standard baselines. We will also state that our framework is agnostic to the search method used for data generation and will provide multiple search strategies in our publicly released code to offer users more flexibility.\n\n----\n## **Thanks Again**\n\nWe sincerely thank you again for your constructive and detailed feedback, which has helped us identify clear pathways to improve our paper. Should you have any further questions or require additional discussion, please don't hesitate to reach out. If we have adequately addressed your concerns, we would be grateful for your consideration in adjusting your evaluation score accordingly."}, {"review_id": "SJ42CENJd3", "reviewer": "Reviewer_C25p", "summary": "This paper introduces SymMaP, a symbolic discovery framework for automatically learning interpretable, efficient symbolic expressions to predict matrix preconditioning parameters for linear solvers. Traditional fixed-parameter methods lack adaptability, while machine-learning-based approaches suffer from high inference costs and poor interpretability. SymMaP combines the benefits of both worlds by using a neural network to explore a symbolic search space and generate compact mathematical formulas describing optimal preconditioning parameters.", "questions": "1. Could it be extended for designing mechanisms for dynamic or adaptive parameter updating, which are often necessary in time-dependent or evolving problems? \n2. Please kindly refer to some questions in the weaknesses above.", "limitations": "yes", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 2, "originality": 4, "strengths_and_weaknesses": "Strengths\n- This paper aims to find a good tradeoff between interpretability and performance. Because the output is an interpretable symbolic model, users could understand and analyze the learned preconditioning strategies. Experimental results show consistent improvements over traditional fixed-parameter baselines on second-order elliptic PDE problems.\n- The learned expressions can be compiled into a shared library and integrated into existing solver frameworks with minimal engineering effort. Meanwhile, this method is computationally efficient at inference time.\n\nWeaknesses:\n- The paper claims that the framework is not strictly tied to SOR and could potentially generalize to other continuous-parameter preconditioners, but I am still worry about its generalization ability. I am wondering how well the learned symbolic expressions perform on linear systems from significantly different domains, like other type of PDE problems. \n- In real industrial scenarios, obtaining large-scale training data through extensive grid search can be computationally expensive or infeasible for new or large problems, which may limit the method’s practicality. \n- It remains unclear whether the method can be effectively applied to real-world problems such as linear programming, where the structure of the linear systems can be very different from PDE discretizations.\n- While the method produces symbolic expressions that are claimed to be interpretable, it is not always clear how these formulas correspond to the physical properties of the original problem. This may limit their practical usefulness in terms of interpretability. In some industrial settings, interpretability might be less prioritized compared to achieving faster convergence. Practitioners may run multiple preconditioners in parallel and select the fastest one, accepting the extra resource use to reduce overall end-to-end time. Therefore, there is a possibility that focusing on interpretability could trade off some potential performance benefits offered by less interpretable, purely empirical approaches.", "was_revised": true, "final_justification": "My final justification for the recommended score is still 4 (borderline accept). My questions have been answered, and I appreciate the authors' detailed responses.", "rebuttal": "**Dear Reviewer,**\n\nWe sincerely thank you for your thoughtful review and for recognizing the \"excellent\" originality of our work. Your insightful feedback is invaluable, and we are grateful for the opportunity to clarify our method's capabilities, address your concerns, and discuss its future potential. Your comments have helped us identify key areas for improvement and further exploration. \n\nWe offer the following points in response to your questions and weaknesses.\n\n---\n\n### **On Generalization and Applicability to Diverse Problem Domains (Weakness 1 & 3)**\n\nAs noted in the paper (Page 5, Line 157), the current implementation of SymMaP takes the parameters of the governing PDE as input features. This approach requires that different problems can be described by a consistent set of parameters. Therefore, if different types of PDEs cannot be modeled using a unified set of parameters, our current model cannot be trained on them simultaneously. This is a limitation of the current feature engineering, not a fundamental constraint of the SymMaP framework itself.\n\nIn principle, SymMaP is a symbolic discovery framework that is agnostic to the problem's origin (e.g., PDE discretization, linear programming). However, as you correctly imply, its practical application is currently bounded by the performance limitations of symbolic regression, which struggles with an excessively high number of input features (e.g., >1000). This makes SymMaP, in its current form, less suitable for problems characterized by a vast number of discrete and unstructured parameters, such as those found in large-scale integer programming.\n\nThis is a very important point, and your question highlights a crucial direction for our future research. We plan to address this by shifting our focus from high-level equation parameters to more fundamental matrix properties. Our next step is to use easily computable matrix features (e.g., matrix norms, trace) and perhaps a curated subset of key matrix elements as inputs to the symbolic learning process. While directly using all matrix elements is currently infeasible, we are actively exploring methods to automatically select the most informative matrix features. This will make SymMaP a truly universal tool for a much wider class of linear systems.\n\n### **On Data Efficiency and Adaptability to Dynamic Problems (Weakness 2 & Question 1)**\n\nYou are correct that SymMaP is a data-driven method. However, a key advantage of symbolic learning is its remarkable data efficiency. As shown in **Table 16 (Appendix, Page 29)**, SymMaP can achieve strong performance with a very small dataset. For instance, learning an effective AMG preconditioner for the Darcy Flow problem required **only 100 data points**. For the main experiments in the paper, we conservatively used a larger dataset (n=1000) to ensure the utmost rigor and stability of our results, but this is not the minimum requirement.\n\nYour suggestion to extend SymMaP to dynamic or adaptive scenarios is an excellent one and aligns perfectly with our vision for future work. We are planning to develop an online, adaptive learning extension of SymMaP. The idea is to first learn a foundational symbolic expression from a small, offline-generated dataset. Then, in a live operational setting, this expression would be continuously fine-tuned using reinforcement learning based on the performance on new, incoming problem instances. This approach would achieve two goals:\n\n1. It would further minimize the reliance on extensive initial training data.\n2. It would allow the preconditioning strategy to adapt to evolving problems, making it highly suitable for time-dependent and other dynamic industrial scenarios.\n\n### **On the Balance of Interpretability and Performance (Weakness 4)**\n\nWe appreciate the reviewer's nuanced perspective on the trade-off between interpretability and performance. We agree that in many industrial applications, computational efficiency is the highest priority. Our goal with SymMaP is to demonstrate that these two objectives are not necessarily mutually exclusive.\n\nThe primary advantage of SymMaP is not to propose a new standalone preconditioner, but rather to serve as a **versatile framework that can be** **embedded** **within any parameterized preconditioner to enhance its performance.** As we demonstrate in our **related work comparison (Appendix E.6, Pages 25-27)**, many learning-based methods focus on optimizing one specific algorithm (e.g., a particular type of AMG [1] or ICC [2]). If that underlying algorithm is a poor fit for a given problem, the learning-based enhancement has a limited ceiling. SymMaP, by contrast, offers broader applicability.\n\nMore importantly, SymMaP is not in opposition to other learning-based methods; it is complementary. Our results in **Table 12 and 13 (Appendix, Page 26)** show this clearly.\n\n- The **best overall performance** across our benchmarks was achieved by combining SymMaP with another learning-based method (\"NeurKItt [3] + SymMaP\"). This shows that SymMaP can be a critical component in a state-of-the-art pipeline.\n- When evaluating **standalone algorithms**, SymMaP consistently delivered the best performance, outperforming other individual methods.\n\nThese results strongly indicate that SymMaP provides a significant performance benefit. Interpretability is an additional, powerful advantage. It allows researchers to analyze the learned formulas (some of which are remarkably simple, as shown on Page 28) to gain deeper insights into the preconditioning process itself. In fact, we are actively pursuing this, with the goal of having SymMaP's discoveries inform the theoretical design of new and improved preconditioning algorithms.\n\n**References:** \n\n[1] Learning to optimize multigrid pde solvers, ICML 2019.\n\n[2] Neural incomplete factorization: learning preconditioners for the conjugate gradient method, TMLR 2024.\n\n[3] Neural krylov iteration for accelerating linear system solving, NIPS 2024.\n\n----\n## **Thanks Again**\n\nWe sincerely thank you again for your constructive and detailed feedback, which has helped us identify clear pathways to improve our paper. Should you have any further questions or require additional discussion, please don't hesitate to reach out. If we have adequately addressed your concerns, we would be grateful for your consideration in adjusting your evaluation score accordingly."}, {"review_id": "w6ldkxR8eC", "reviewer": "Reviewer_azLf", "summary": "This paper proposes SymMaP, a symbolic discovery framework for predicting efficient preconditioning parameters in linear system solvers. SymMaP uses a neural network to search the space of symbolic mathematical expressions, learning concise and interpretable formulas that map problem-specific features (such as PDE coefficients) to optimal preconditioning parameters. Unlike traditional fixed-parameter or black-box ML approaches, SymMaP produces symbolic expressions that are fast to evaluate, generalize well across problem instances, and are straightforward to deploy in CPU-based scientific computing environments. The framework is validated on a range of PDE-derived linear systems and preconditioners (SOR, SSOR, AMG), showing consistent improvements in computation time and condition number over both default and learning-based baselines.", "questions": "- The Font of Figure 1 has an inconsistency (non-Times New Roman) and misaligned \"78%\" labels (likely from unrounded decimals). Standardize fonts and round percentages uniformly.\n- Why omit benchmarks against high-performance solvers (e.g., Gurobi) or low-level libraries (e.g., Intel MKL)? If SymMaP targets generic linear systems, validate it on non-PDE matrices (e.g., SuiteSparse).\n- Does SymMaP require instance-specific training? If features are different (Darcy Flow Problem, Second-order Elliptic Partial Differential Equation, Biharmonic Equation, Poisson Equation, ...), what architectural changes ensure cross-instance generalization?", "limitations": "Yes", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 4, "originality": 3, "strengths_and_weaknesses": "Strengths:\n- The learned symbolic expressions are concise, interpretable, and incur almost no computational overhead, making them suitable for integration into existing solver libraries.\n- SymMaP demonstrates strong generalization across various PDE types, matrix sizes, and preconditioning methods, outperforming both fixed-parameter and neural network baselines in terms of both accuracy and speed.\n-  The symbolic expressions can be compiled and directly embedded into CPU-based solvers (e.g., PETSc), addressing a key limitation of ML-based approaches that often require GPU acceleration.\n- The approach enhances transparency in scientific computing by providing explicit formulas that can be analyzed for theoretical soundness and safety.\n\nWeaknesses:\n- While the paper compares with basic MLPs and some recent learning-based methods, it does not benchmark against the latest large-scale neural symbolic regression models or hybrid approaches that might offer stronger baselines.\n- Evaluated exclusively on PDE-derived systems; untested on broader linear algebra benchmarks (e.g., network problems, optimization tasks).", "was_revised": true, "final_justification": "Thank you for your detailed response. Your reply has addressed the questions I raised and clarified my understanding of your work.\n\nMy initial concerns have been resolved, and I do not have additional issues after the rebuttal and discussion.\n\nTherefore, I maintain my original positive score of 4.", "rebuttal": "**Dear Reviewer,**\n\nWe sincerely thank you for your thoughtful review. Your insightful feedback is invaluable, and we are grateful for the opportunity to clarify our method's capabilities, address your concerns, and discuss its future potential. Your comments have helped us identify key areas for improvement and further exploration.\n\nWe provide detailed responses to your points below.\n\n--- \n\n### **On Comparison with Advanced Symbolic** **Regression** **Baselines (Weakness 1)**\n\nWe thank the reviewer for this valuable suggestion. Our paper's primary contribution is the novel framework that applies symbolic discovery to the problem of matrix preconditioner parameterization, rather than proposing a new symbolic regression algorithm itself.\n\nOur main goal in the comparison with a standard MLP was to validate a core premise: that concise symbolic expressions possess sufficient expressive power to accurately model the relationship between problem features and optimal parameters in this domain. Our results confirm this viability.\n\nWe fully agree that the symbolic learning module within SymMaP is modular and can be replaced. The prefix notation + RNN approach we used is a well-established and effective paradigm. We believe our experiments demonstrate its strong performance for this task. We are excited by the prospect that future breakthroughs in neural symbolic regression can be readily integrated into the SymMaP framework, potentially unlocking even greater performance. We will clarify this aspect in the final version.\n\n### **On Generalization, Evaluation Scope, and Architectural Design (Weakness2, Questions 2 & 3)**\n\nWe thank the reviewer for these insightful questions regarding the evaluation scope and the mechanism for cross-instance generalization. We are happy to elaborate on our methodology and future vision.\n\n1. On the Choice of Problem Domain and Benchmarks (PDEs & PETSc):\n\nOur choice to focus on PDE-derived systems and use PETSc as the solver environment was a deliberate one, intended to ground our work in the context of recent, high-impact research. Our experimental design is closely aligned with that of prominent papers in the field, such as [1] (NeurKItt, NeurIPS 2024) and [2](ICLR 2024).\n\n- **PETSc** is a highly-optimized, professional library for solving linear systems, particularly those arising from differential equations. It is the foundation for industrial software like FEniCS and deal.II, making it a relevant and rigorous benchmark environment. Specialized solvers like Gurobi are designed for different problem classes (e.g., linear programming), while low-level libraries like MKL do not provide the same high-level preconditioning framework.\n- **Differential equations** represent a vast and critical class of problems in science and engineering. Many optimization problems, such as the **Markov chain** **(Boltzmann Motion)** example discussed on **Page 18 (Line 677)** and evaluated on **Page 26**, are also formulated or solved using differential equations.\n- In principle, SymMaP is a symbolic discovery framework that is agnostic to the problem's origin (e.g., PDE discretization, linear programming). However, as you correctly imply, its practical application is currently bounded by the performance limitations of symbolic regression, which struggles with an excessively high number of input features (e.g., >1000). This makes SymMaP, in its current form, less suitable for problems characterized by a vast number of discrete and unstructured parameters, such as large-scale integer programming problems in optimization tasks.\n\n2. On the Mechanism for Cross-Instance Generalization:\n\nAs noted in the paper (Page 5, Line 157), SymMaP is a data-driven framework that learns a mapping from a problem's input features to its optimal preconditioner parameters. For this learning process to succeed, the features must be semantically consistent across the dataset. Therefore, a single learned expression is specific to a class of problems that share a common parametric description (e.g., all second-order elliptic PDEs, characterized by their coefficients).\n\nThis requirement explains why we cannot directly train or test on a benchmark like SuiteSparse, which is a collection of individual, structurally diverse matrices without a shared, low-dimensional feature space. The architecture of SymMaP itself is general; it does not require changes between problem types. However, it does require a consistent feature representation for each training task.\n\n3. On the Path to Broader Generalization (Future Work):\n\nYou have raised an excellent point about extending SymMaP to handle a wider variety of linear systems and enable joint training across different problem types. Addressing this is a core part of our future research agenda. Our planned approach is to move from high-level equation parameters to more fundamental matrix properties as input features. We plan to learn the symbolic relationship between low-cost, obtainable matrix features (e.g., matrix norms, trace, sparsity patterns) and the optimal preconditioner parameters. This would create a universal feature space, breaking the dependency on problem-specific coefficients and enabling generalization across a much broader range of linear algebra problems. While directly using all matrix elements is currently infeasible, we are actively exploring methods to automatically select the most informative matrix features. This will make SymMaP a truly universal tool for a much wider class of linear systems.\n\n**References:** \n\n[1] Neural krylov iteration for accelerating linear system solving, NIPS 2024.\n\n[2] Accelerating data generation for neural operators via krylov subspace recycling, ICLR 2024\n\n### **On Revisions to Figure 1 (Question 1)**\n\nThank you for your sharp eye and for pointing this out. We have re-examined Figure 1 and confirmed that the misalignment of the \"78%\" labels was indeed due to a decimal rounding issue in the plotting script. We sincerely apologize for this oversight. In the final version of the paper, we will correct this and ensure all fonts are standardized to Times New Roman for consistency and professionalism.\n\n----\n## **Thanks Again**\n\nWe sincerely thank you again for your constructive and detailed feedback, which has helped us identify clear pathways to improve our paper. Should you have any further questions or require additional discussion, please don't hesitate to reach out. If we have adequately addressed your concerns, we would be grateful for your consideration in adjusting your evaluation score accordingly."}, {"review_id": "Fg3XE5lfGn", "reviewer": "Reviewer_A9dG", "summary": "This paper proposes a novel method for selecting parameters of preconditioning techniques for solving linear equations. The approach involves first generating a dataset by determining optimal parameters via grid search, then training a recurrent neural network within a reinforcement learning framework to produce symbolic expressions. These symbolic expressions are subsequently used to select preconditioning parameters, thereby improving the efficiency of solving linear systems. Since the symbolic representations are lightweight and can be easily computed on CPUs, they are well-suited for integration into existing numerical software. Moreover, their interpretability may assist in better understanding the problems and in facilitating mathematical analysis.", "questions": "- Why is the AMG case evaluated using condition number rather than computation time? The computation time would also be of interest.\n- Figure 1 compares with the optimal parameter, but what is the performance of the optimal parameter in the other experiments (Tables 1, 2, and 3)?\n- In line 118, AMG is cited as an example of low-cost matrix decomposition, but is this appropriate? ICC or ILU may be more suitable examples?\n- The input parameters in this method are low-dimensional, but can the method accommodate more complex problems? For instance, how would it handle variations in domain geometry or complex boundary conditions?\n- Many preconditioning parameters are discrete rather than continuous. Can this method be extended to handle such discrete parameter spaces?", "limitations": "yes", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "# Strengths\n- The approach of accelerating and improving traditional numerical algorithms and software through machine learning is practical and highly relevant to real-world problems.\n- The proposed method is original. The idea of using RL and RNN to produce lightweight symbolic models for parameter selection is interesting.\n- Performance is evaluated using CPU runtime within commonly used software, and the superiority of the proposed method over baselines and alternative approaches is demonstrated using this practical metric.\n\n# Weaknesses\n- The cost of generating training data is quite high. Creating a single data point requires solving the problem multiple times through grid search, which may become a bottleneck for larger problems due to inefficiency.\n- The paper does not sufficiently explain the necessity or advantage of interpretability in selecting preconditioning parameters. Since solution accuracy can be controlled via tolerance, it may be acceptable for preconditioning to function as a black box. The paper argues that opaque prediction may fail to satisfy theoretical constraints (e.g., 0 < \\omega < 2), but if the constraint is simply a matter of range, this could be enforced via a suitable activation function. Furthermore, symbolic expressions are not always interpretable. In fact, many expressions shown in Appendix E.1 appear difficult to interpret.\n- Even if the need for interpretability is established, the paper does not compare the proposed method against other interpretable and lightweight approaches. For example, how does the proposed RL+RNN approach compare in interpretability and performance to directly training simple models such as linear regression, polynomial regression, or decision trees? Also, if symbolic regression is the goal, how does this method compare to directly optimizing symbolic expressions without using an RNN? Such comparisons would further clarify the advantages of the proposed approach.", "was_revised": true, "final_justification": "Thank you for your detailed response. My questions have been addressed, and I updated the score.", "rebuttal": "**Dear Reviewer,**\n\nWe are sincerely grateful for your detailed and thoughtful review of our paper. Your constructive feedback is invaluable, and it has helped us identify several areas where we can strengthen our explanations and empirical validation. We hope that our detailed responses, including several new experiments conducted specifically to address your questions, will clarify the contributions of our work and merit a re-evaluation.\n\n### **On the Efficiency of Training Data Generation (W1)**\n\nYou are correct that SymMaP is a data-driven method. While an grid search can be computationally intensive, a key advantage of leveraging symbolic learning is its remarkable data efficiency. As we show in Table 16 (Appendix, P29), SymMaP can learn high-performing expressions from a very small dataset. For instance, learning an effective AMG for the Darcy Flow problem required only 100 data points to achieve strong results. For the experiments in our main paper, we conservatively used a larger dataset of 1000 samples to ensure the utmost rigor, but this should not be seen as the minimum requirement.\n\nFurthermore, your comment has inspired us to formalize our vision for future work. We plan to develop an online, adaptive learning framework. This would involve learning a foundational symbolic expression from a small initial dataset and then fine-tuning it in a live environment using RL based on the performance of new, incoming problem instances. This would not only further reduce the initial data requirement but also make SymMaP ideally suited for dynamic, real-world application scenarios.\n\n### **On the Value and Necessity of Interpretability (W2)**\n\nWe agree that for some applications, a black-box model is perfectly acceptable, and we can certainly enforce known constraints (like a parameter range) using methods such as activation functions. However, the value of interpretability in our scientific computing context is twofold:\n\n1. Discovering Unknown Constraints: While known constraints can be hard-coded, the primary advantage of an interpretable model is the ability to analyze its structure to identify potential unknown failure modes or numerically unstable regimes before deployment. A symbolic formula allows for direct mathematical analysis in a way a neural network does not.\n2. Fostering Scientific Insight: More importantly, interpretability serves as a bridge between machine learning and scientific discovery. By analyzing the simpler symbolic expressions discovered by SymMaP, researchers can gain a deeper understanding of the preconditioner itself. As we show on P21-22, some of the learned expressions are remarkably concise and reveal the underlying mathematical dependencies. We are actively pursuing this direction, with the goal of using the insights from these formulas to help design theoretically better preconditioning. \n\nWe will emphasize this motivation more clearly in the revised manuscript.\n\n### **Comparison with Other Interpretable Models (W3)**\n\nThis is an excellent suggestion. To better situate SymMaP and clarify its advantages over other lightweight, interpretable models, we have conducted a new experiment as you proposed.\n\nThe key advantage of symbolic regression over methods like linear or polynomial regression is its flexibility in both feature selection and functional form. Symbolic regression explores a much richer space of expressions. In fact, by constraining the library of operators (e.g., using only `+`, `*`, and polynomial terms as discussed in Section 4.2 (P5, Line 176) and Section 5.4 (P9,Line 331)), our framework can effectively subsume these simpler models; they are special cases of the broader search space SymMaP explores.\n\nTo provide a direct performance comparison, we ran an ablation on the second-order elliptic PDE problem (SOR, matrix size 40000), comparing SymMaP against Linear Regression and Quadratic Polynomial Regression.\n\n|  | Solution Time (s) |\n| -- | -- |\n| Linear | 20.9  |\n| Quadratic Polynomial | 17.0   |\n| SymMaP   | 15.8   |\n\nThe results clearly show that the flexible search of symbolic regression discovers a superior model, leading to faster solution times. The simpler models, while interpretable, lack the expressive power to capture the true optimal parameter landscape. We will add this experiment and analysis to the final version of the paper.\n\n### **On the Choice of Metric for AMG (Q1)**\n\nFor a fixed class of PDEs and a given discretization scheme, the computation-time-optimal threshold parameter often lies in a very narrow, near-zero range (e.g., 0-0.02 for the elliptic PDEs in our study). This makes it a less dynamic learning target.\n\nIn contrast, the threshold that minimizes the **condition number** varies more significantly with the PDE's specific coefficients, presenting a more meaningful optimization problem for our learning framework to solve. We acknowledge that a real-world application would likely balance solution time and numerical stability using a hybrid metric. However, to avoid introducing an arbitrary, custom metric in our main results, we chose to focus on the more dynamic and well-defined objective of minimizing the condition number for the AMG experiments.\n\nAs a proof of concept for custom metrics, we show in our response to **Question 5** that SymMaP performs excellently when optimizing a hybrid (time + condition number) objective for ICC.\n\n### **Performance of the Optimal Parameter Baseline (Q2)**\n\nThe \"Optimal Parameter\" represents the theoretical performance ceiling for any parameter-selection method and provides crucial context. We have compiled the data for the SOR experiments from Table 1 below.\n\n| Time (s)  | Theoretical Optimal | None | Default ω=1 | ω=0.2  | ω=1.8  | Optimal constant | SymMaP |\n| -- | -- | -- | -- | -- | -- | -- | -- |\n| Darcy Flow | 7.94  | 33.1   | 13.5  | 17.5  | 9.91 | 9.54 | 8.5 |\n| Biharmonic | 0.92  | 7.67   | 2.04  | 4.86  | 1.6 | 1.31 | 1.24   |\n| Poisson | 0.013  | 0.0412 | 0.0195  | 0.0215 | 0.0195 | 0.0138 | 0.0135 |\n| Thermal  | 0.0522  | 0.223  | 0.0598   | 2.07 | 1.18 | 0.594   | 0.0576 |\n\nAs the table shows, SymMaP's performance consistently approaches this theoretical ceiling, significantly closing the gap left by default or fixed-constant strategies. We will add the optimal parameter results for all relevant experiments to the final version of the paper to provide this valuable context.\n\n### **On the Example of a Low-Cost Decomposition (Q3)**\n\nYou are correct to point this out. Citing AMG as an example of a \"low-cost matrix decomposition\" in **Line 118** could be clearer. While AMG is computationally cheaper than a full factorization, ICC or ILU are indeed more direct and classic examples of low-cost decomposition-based preconditioners. This was meant only as an illustrative example. We will revise this sentence to include ICC and ILU for better clarity and accuracy.\n\n### **Applicability to Complex Geometries and Boundary Conditions (Q4)** \n\nThis is a very important question about the scalability of our feature representation. In principle, SymMaP is agnostic to the problem's geometry. For problems with varying boundaries, we can parameterize the boundary shape and include these parameters as inputs to the model.\n\nTo demonstrate this, we conducted a new experiment on a **second-order elliptic PDE with a variable triangular domain**. The triangle's vertices were randomly sampled, and these vertex coordinates were used as input features for SymMaP to predict the optimal SOR relaxation factor.\n\n|  | None | Default ω=1 | ω=0.2 | ω=1.8 | Optimal Constant | SymMaP |\n| -- | -- | -- | -- | -- | -- | -- |\n| Time (s) | 16.7 | 8.3  | 9.04  | 9.01  | 7.65 | 7.72  |\n\nThe results show that SymMaP performs robustly, achieving near-optimal performance even when the domain geometry varies.\n\nHowever, we acknowledge the limitation you insightfully identified: our feature engineering approach is currently limited to problems that can be described by a low-dimensional parameter vector ( <1000 parameters). For highly complex geometries or boundary functions, this would be a bottleneck. This is a key focus of our future work. We plan to move from high-level equation parameters to more fundamental matrix properties (e.g., matrix norms, trace) as inputs. This would create a universal feature space, making SymMaP applicable to a much broader and more complex class of problems. We will expand the limitations section to discuss this path forward.\n\n### **Extension to Discrete Parameter Spaces (Q5)**\n\nOur framework can indeed be extended to handle discrete parameters, such as the fill-in level for ICC. To demonstrate this, we performed a new experiment on the second-order elliptic PDE (matrix size 40000) with the **ICC**, where the goal is to find the optimal (integer) fill-in level. Since increasing the fill-in level improves the condition number but increases computation time, we defined a hybrid objective: `Metric = 0.03 * Condition_Number + Time(s)`.\n\n| Method | None  | Default | Optimal Constant | SymMaP |\n| --- | --- | -- | -- | -- |\n| Metric | 227.6 | 24.0 | 22.8  | 22.2 |\n\nThe results confirm that SymMaP can successfully optimize discrete parameters under a custom objective function. This demonstrates the flexibility of our framework beyond the continuous parameters presented in the main paper. We will add this experiment and a discussion on handling discrete parameters to the final version.\n\n----\n## **Thanks Again**\nWe sincerely thank you again for your critical and highly valuable feedback. Your review has pushed us to substantially improve the paper by conducting new experiments, clarifying our contributions, and strengthening our discussion of the method's scope and limitations. Should you have any further questions or require additional discussion, please don't hesitate to reach out. If we have adequately addressed your concerns, we would be grateful for your consideration in adjusting your evaluation score accordingly."}]

Author rebuttal:

Meta-review: This paper considers the important task of matrix preconditioning, a crucial step in linear system solving. While matrix preconditioners are frequently designed in an ad-hoc manner, the paper proposes to symbolically learn preconditioners to automate this process. Experimental results support the proposed method over prior strategies.

Reviewers and I find the problem relevant, and the proposed approach original and of interest to the community. I'm happy to recommend acceptance.
Decision: Accept (poster)