SymMaP: Improving Computational Efficiency in Linear Solvers through Symbolic Preconditioning

We appreciate the reviewer's insightful and valuable comments. We respond to each comment and sincerely hope that our rebuttal could properly address your concerns. If so, we would deeply appreciate it if you could raise your score. If not, please let us know your further concerns, and we will continue actively responding to your comments and improving our submission.

Weaknesses 1 & 2

SymMaP enhances traditional preconditioning algorithms by discovering symbolic expressions for optimal parameters. In our experimental design, we intentionally avoid comparisons with other learning-based approaches for the following reasons:

1. Focus on Preconditioning Enhancement: SymMaP is designed to improve traditional preconditioners without altering their core algorithms. A fair comparison would require other learning-based methods that also preserve the original preconditioning process—yet no such methods currently exist.
2. CPU Compatibility: As noted the Introduction (line 26, right column), linear solver environments are predominantly CPU-based. SymMaP’s symbolic expressions integrate seamlessly into CPU workflows, whereas neural network-based approaches lack efficient CPU deployment.
3. Generality: Section 3.2 (line 152, left column) highlights SymMaP’s adaptability to diverse preconditioners and linear solvers. In contrast, existing learning-based methods either: ‘combine with specific preconditioners’ (e.g., [4], [6]), or ‘target narrow use cases’ (e.g., [8], [2] for ICC-preconditioned CG; [1] exclusively for CG).
4. Interpretability & Safety: Numerical algorithms demand rigorous analysis (Section 3.2, line 129, right column). Opaque predictors (e.g., neural networks) risk violating theoretical constraints (e.g., avoiding ω ≈ 0 or 2 in SOR). SymMaP’s symbolic expressions enable proactive avoidance of such issues through analytical guarantees.

To further address your concerns, we can compare SymMaP with other learning-based preconditioning methods (disregarding the aforementioned factors). Let us first categorize existing works:

1. Category 1: Methods that enhance Krylov solvers by adding optimization modules or modifying their iteration process (e.g., [4], [6]). These works target different components than SymMaP and could potentially be combined with our approach for further improvements.
2. Category 2: Methods that use neural networks to predict a preconditioning matrix (e.g., [8], [1], [7]).

• Experiment: SSOR preconditioning for elliptic PDEs (same settings as the paper):

  • method	none	PETSc default SSOR	[4]	[4]+SSOR	[4]+SymMaP+SSOR	[1]	[8]	SymMaP+SSOR
    time (s)	23.9	10.5	11.2	7.5	5.45	8.6	8.4	7.7

  • Training times: [1] (0.5h), [8] (0.5h), [4] (2h), SymMaP (0.25h).

  • Key takeaway: SymMaP outperforms standalone learning-based methods and can combine with them for further gains.

• We reiterate that direct comparisons are inherently unfair—akin to contrasting symbolic (Mathematica) and numerical (MATLAB) PDE solvers. If reviewers identify similar works to ours, we welcome suggestions and will conduct comparative experiments accordingly.

[8] Learning from Linear Algebra: A Graph Neural Network Approach to Preconditioner Design for Conjugate Gradient Solvers

Weaknesses 3

• Dataset size impact: SymMaP is fundamentally distinct from neural network training, as it seeks a symbolic expression to approximate the mapping relationship in the data rather than optimizing parameters. Our preliminary tests show that even small datasets (~100 samples) suffice to learn high-performance symbolic expressions. For consistency, all experiments in the paper use a dataset size of 1,000.
• Additional experiments: To further address this concern, we conducted supplementary tests on SOR preconditioning for elliptic PDEs (settings identical to the main experiments), varying the training set size (100, 300, 500, 1,000). Results show that regardless of the size of the data set, SymMaP can find the best symbolic expression in the experiment in our paper within 1000 seconds.

Weaknesses 4

• Current design: Experiments intentionally avoid presupposing symmetry among input parameters, reflecting real-world scenarios where such relationships are unknown.
• Future direction: We agree this is valuable and plan to develop a module that analyzes mathematical properties (e.g., symmetries) to constrain symbolic search spaces, improving efficiency and performance.

We appreciate the reviewer's insightful and valuable comments. We respond to each comment and sincerely hope that our rebuttal could properly address your concerns. If so, we would deeply appreciate it if you could raise your score. If not, please let us know your further concerns, and we will continue actively responding to your comments and improving our submission.

Methods And Evaluation Criteria 1

• In preliminary tests, we observed that for certain preconditioners (e.g., SOR), the optimal preconditioning parameters are largely independent of the matrix size(resolution) for a given problem. Thus, we fixed the matrix size in our experiments for consistency. To address your concern, we conducted additional experiments:
  • We evaluated SOR preconditioning for a second-order elliptic PDE, testing matrix sizes of ( $1 \times 10^4, 2 \times 10^4, \dots, 6 \times 10^4$ ) while keeping other settings identical to the main experiments.
  • The results show that the optimal parameters for the same problem but different sizes vary negligibly (difference < 0.01).
• This confirms that matrix size has minimal impact on SOR’s optimal parameters. For preconditioners where size does affect parameters, SymMaP can also easily incorporate the size as an additional input.

Methods And Evaluation Criteria 2

• To further validate SymMaP’s versatility, we conducted experiments on AMG preconditioning with an SOR smoother, optimizing two parameters simultaneously: AMG’s threshold parameter θ (default in PETSc: 0) and SOR relaxation factor ω.

  • All other settings match the AMG preconditioning and second-order elliptic problem in the paper. We jointly optimize both parameters, with the condition number as the evaluation metric.

    • method	None	θ= 0 ω=1	θ= 0 ω=0.1	θ= 0 ω=1.9	Optimal constant	SymMaP
      Condition number	6792	163	168	163	159	156

  • This demonstrates SymMaP’s ability to handle multiple preconditioner parameters. We will include detailed results and analysis in the final version.

Essential References Not Discussed

Please see the response to reviewer Wyqu's Weaknesses 1 & 2

Q1

• As noted in Table 1’s caption (line 275), SymMaP 1 and 2 are the two highest-scoring expressions identified during symbolic learning.
• The disparity arises from the stochastic nature of the learning process: Symbolic learning iteratively refines expressions by introducing new operators, which may non-monotonically improve performance (some additions help, others degrade it). Thus, the two expressions represent local optima with similar reward scores but differences in structure and performance due to the randomness of operator selection.

Q2

To explain the significant variation in condition numbers in Table 3 (row 343), we analyze from the following perspectives:

1. AMG’s Impact on Condition Number:AMG reduces matrix condition numbers through a multigrid strategy, leveraging coarse-grid correction and high-frequency error smoothing. On coarse grids, restriction and interpolation operators transform low-frequency errors (small eigenvalues) into high-frequency ones, which are then rapidly damped by fine-grid smoothing (e.g., Gauss-Seidel). This multiscale decomposition concentrates the eigenvalue spectrum, significantly lowering the condition number (ratio of largest to smallest eigenvalues).

2. Role of Threshold Parameter θ:θ controls strong connectivity during coarse-grid generation. If the coupling strength (off-diagonal entries) exceeds θ, the connection is retained; otherwise, it is discarded, affecting grid sparsity and approximation accuracy.

   1. Too small θ: Retains excessive weak connections, leading to dense coarse grids that inadequately address low-frequency errors (small eigenvalues).
   2. Too large θ: Over-sparsifies the grid, losing critical connections and increasing approximation error, hindering high-frequency error (large eigenvalue) reduction.

3. Matrix Properties:The matrices derive from differential operators, where the smallest eigenvalue (in matrix norm) reflects the operator’s minimal eigenvalue. Such matrices typically require stronger large-eigenvalue suppression, favoring smaller θ (0–0.5). The wide θ range in Table 3, tailored to operator specifics, explains the condition number disparity.

4. Mitigation Strategy:While condition numbers vary sharply, they change continuously with θ. Moreover, optimal θ values vary smoothly with operator parameters. Here, SymMaP effectively approximates these optimal parameters, resolving instability. This scenario motivated our studies.

We appreciate the reviewer's insightful and valuable comments. We respond to each comment and sincerely hope that our rebuttal could properly address your concerns. If so, we would deeply appreciate it if you could raise your score. If not, please let us know your further concerns, and we will continue actively responding to your comments and improving our submission.

Methods 1

And, as mentioned above, there is a lack of experiments with more generic datasets. (e.g. [1] [2] [3])

• In our main experiments (Section 5.1, line 374), we evaluate 5 classes of physical equations. For context: The cited works study fewer equation classes: [1] (2 classes), [2] (3 classes), [3] (5 classes, two of which are special cases of second-order elliptic equations).
• Other related works are even more limited: [4] (3 classes), [6] (2 classes), [7] and [8] (1 class each).
• Among our five classes, Poisson’s equation is discussed in [1–4, 6, 7], thermal problems in [2, 4], and second-order elliptic equations in [3].
• These datasets represent distinct, scientifically significant problem classes with varied mathematical properties (see Sections 5 and D.1). Crucially, SymMaP outperforms baselines across all dataset-preconditioner combinations (Tables 1–3), demonstrating its generalization capability. We will add citations to these works in the final version. If the reviewer suggests additional equations to test, we are happy to include them.

Methods 2

Comparisons with other methods in the literature are also lacking (e.g. [1][2][3]).

• Please see the response to reviewer Wyqu's Weaknesses 1 & 2
• Additional Note: You mentioned [3], a neural network-based PDE solver. Similar works include PINN [5], which employs neural networks as optimizers to solve PDEs via gradient descent. These methods do not involve traditional linear solvers or preconditioning optimization.

Theoretical Claims

• As noted in Section 3.1 (line 100) and Figure 1 (line 67), we observe that SOR, SSOR, and AMG preconditioner parameters exhibit continuous relationships with performance metrics, typically with one or two local minima. Thus, grid search effectively approximates the optimal parameters.
• Theoretical derivation of optimal preconditioner parameters is often intractable, justifying our grid-search approach.

Q 1

• We agree. As noted in Section 3.2 (C1, line 154), preconditioner objectives vary by scenario:
  • For stable systems, minimizing runtime is prioritized.
  • For ill-conditioned systems (e.g., certain elliptic PDEs), reducing condition numbers is critical to ensure solvability.
• Our AMG experiments (Section B.2.2) focus on threshold parameters, which filter small values during graph aggregation. For some problems (e.g., second-order elliptic PDEs), setting this parameter to 0 minimizes runtime, so we prioritize stability here.

Q 2

• SOR: The relaxation factor in SOR ranges between (0, 2). For well-conditioned matrices, larger values (e.g., 1.5–1.9) are typically stable. For matrices with significant off-diagonal elements (e.g., nonlinear PDEs from turbulence), smaller values (e.g., 0.1–0.5) are preferred. Note that SOR reduces Gauss-Seidel when the factor is 1, which is also the default in libraries like PETSc. Thus, we test three fixed values: 0.2, 1, and 1.8. The same logic applies to SSOR.

• AMG: Our preliminary tests show that the optimal threshold for our physical equations lies below 1.5. A threshold of 0 (or lower) disables coarsening, which is the default in PETSc and similar libraries. We therefore test three fixed values: 0, 0.2, and 0.8.

Q 3

• Tolerance is the relative convergence criterion: ∥b − Ax∥/r₀, where r₀ is the initial residual. We apologize for the ambiguity and will define this explicitly in the final version.

Q 4

• SymMaP uses RNN-based symbolic learning. Training time: 800s (non-polynomial symbols) or 2600s (with polynomials) (see Section "Computational Time," line 926).
• In Section 5.2 (line 376), we compare SymMaP to an MLP (trained for 1.5 hours to convergence).

[1] Learning from Linear Algebra: A Graph Neural Network Approach to Preconditioner Design for Conjugate Gradient Solvers

[2] Learning Preconditioners for Conjugate Gradient PDE Solvers

[3] A Finite Element based Deep Learning solver for parametric PDEs

[4] Neural Krylov iteration for accelerating linear system solving

[5] Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations

[6] Neural incomplete factorization: learning preconditioners for the conjugate gradient method

[7] Learning to Optimize Multigrid PDE Solvers

[8] Learning Neural PDE Solvers with Convergence Guarantees

We appreciate the reviewer's insightful and valuable comments. We respond to each comment and sincerely hope that our rebuttal could properly address your concerns. If so, we would deeply appreciate it if you could raise your score. If not, please let us know your further concerns, and we will continue actively responding to your comments and improving our submission.

• Claims And Evidence 1

  • As noted, data generation can be time-consuming. However, symbolic learning in our framework requires significantly less data compared to other AI-based approaches. For details, please refer to our response to Reviewer Wyqu’s Weaknesses 3.

Claims And Evidence 2

To address this concern, we provide additional experimental results at https://anonymous.4open.science/r/rebuttal2-534E/rebuttal2.2.png. The data demonstrate the stability of our algorithm’s performance. We will include a comprehensive distributional analysis (e.g., median, quartiles) in the final version.

Claims And Evidence 3

• To clarify, we use the SOR-preconditioned second-order elliptic PDE experiment as an example:
• As stated in Line 795 (Page 15), the input parameter α for SymMAP is sampled uniformly at random from its defined interval.
• In this experiment, the average discrepancy between the predicted symbolic ω and the ground truth is 0.03, indicating no memorization or overfitting.

Claims And Evidence 4

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, we conducted additional experiments:

• Settings: Equation: Second-order elliptic, grid size 40,000. Test range: α ∈ (−2, 2), coupling term ∈ (−0.5, 0.5) (vs. training range: α ∈ (−1, 1), coupling term ∈ (−0.01, 0.01)). Tolerance: 1e-3.

  • None	PETSc default 1	Fixed constant 0.1	Fixed constant 1.9	Optimal constant	SymMaP
    16.98	4.04	1.26	15.2	0.94	0.86

• The results confirm strong generalization. We will include a detailed analysis in the final version.

Generalization Across PDE Types and Discretization Schemes

As correctly pointed out, SymMAP searches for mappings from input parameters to optimal preconditioner parameters. Changing the PDE type or discretization scheme alters the optimal parameters, limiting direct generalization.

However, our interpretable expressions reveal that optimal parameters depend on specific matrix features (e.g., sparsity, diagonal dominance). In future work, we aim to generalize SymMAP by directly ingesting matrix properties, enabling adaptation to diverse PDEs and discretizations.

Experimental Designs Or Analyses 2

Moreover, trajectories of iterations could be provided for samples of PDEs.

We include iteration trajectories for the SOR-preconditioned biharmonic equation https://anonymous.4open.science/r/rebuttal2-534E/rebuttal2.1.png, demonstrating SymMAP’s superior convergence. A full analysis will be added to the final version.

Essential References Not Discussed

The paper fails to consider adaptive methods that tune ω as the residual Ax(k)−b shrinks over iterations. For instance [1].

We appreciate this insightful suggestion. The work in [1] presents an elegant approach for adaptively adjusting the SOR relaxation parameter ω based on residual reduction. We will explore further in future work to potentially integrate with our framework, and we will ensure proper citation.

However, we note a fundamental distinction between the application scenarios:

1. SOR as a Standalone Solver (as in [1]): • This is a fixed-point iteration method where optimal ω can be derived theoretically (albeit computationally expensive).
2. SOR as a Preconditioner (our focus): • Here, SOR assists Krylov subspace methods (e.g., GMRES, MINRES, CG). The optimal ω lacks theoretical guidance and depends on the solver choice (e.g., GMRES-preconditioned ω ≠ CG-preconditioned ω). • Our grid-search-based approach addresses this gap, as no existing adaptive methods target preconditioner tuning in this context.

Thus, while [1] is highly relevant to SOR solvers, its direct applicability to preconditioned Krylov methods remains limited.

Other Comments Or Suggestions

There are typos in Figure 3.

• We apologize for this oversight and will correct the typographical errors in the final manuscript.

[1] Adaptive SOR methods based on the Wolfe conditions, 2019, https://dl.acm.org/doi/10.1007/s11075-019-00748-0