A Greedy PDE Router for Blending Neural Operators and Classical Methods

The gist: The proposed method presents an approximate greedy router that efficiently mimics a greedy approach to solver selection for hybrid PDE solvers, achieving comparable error levels in substantially fewer iterations while exhibiting more stable error decay.

How it works

The core idea is to design a general hybrid PDE solver where a routing rule chooses, at each iteration, from a set of classical solvers and neural operators. This update generalizes the classical iterative update by allowing Cj to be non-linear and enabling the preconditioning function to be adaptively chosen at every step. The objective is then defined as minimizing the error norm after T steps, which involves selecting a sequence S = (S1,..., ST) that minimizes h(S).

The Greedy Strategy and Theoretical Guarantees

The greedy rule selects the solver with the largest immediate error reduction, which can be near-optimal when the final error satisfies (weak) supermodularity. The paper introduces a sequence-based notion of weak supermodularity to characterize the suboptimality of greedy solutions for weakly supermodular and postfix monotonic sequence functions. Theorem 4.1 establishes that if the objective function h is weakly supermodular with respect to the optimal solution O and postfix monotone, then the greedy solution S T has a constant-factor approximation guarantee to the optimal strategy.

Learning an Approximate Router

Since true error is not observed at test time, a convex surrogate loss is introduced which enables the learned model to imitate the greedy router without access to true error information. The routing loss lroute is defined as minimizing the expected squared L2(D) error. This loss can be formulated using a surrogate function Ψ, which is Bayes consistent with respect to lroute and upper bounds the original loss up to a constant factor.

Empirical Performance and Results

Empirical results on the Poisson and convection–diffusion equations show that the method consistently reduces final error and area-underthe-curve (AUC) of the error trajectory relative to single-solver baselines and existing hybrid approaches such as HINTS. Table 1 shows that the learned greedy router achieves the lowest final error and AUC across all solver families, outperforming their single-solver and HINTS counterparts. Furthermore, the results demonstrate that the router learns qualitatively similar strategies to the oracle, showing near-monotone convergence rather than oscillatory behavior seen in HINTS.

Solver Ensemble Effectiveness

The experiments also show that enlarging the solver set yields gains beyond pairwise configurations. The learned router achieves these gains by learning to route to each member of the ensemble for a non-trivial proportion of the iterates, and it is observed that the neural operator is selected only rarely across all configurations. This reflects the strength of the numerical solvers rather than a limitation of our method.

Conclusion

The proposed approximate greedy router is an adaptive method for selecting solvers that efficiently minimizes the final error of a PDE in iterative methods, opening many directions for future work. The results support the claim that enlarging the solver set provides measurable benefits and that our learned policy effectively captures the resulting improvements. Finally, it is observed that all mode-wise errors/AUC reduce with the inclusion of more solvers. The method's effectiveness is demonstrated across various PDE types and boundary conditions. The limitations include the higher on-time computational cost required to learn the greedy routing strategy. Additionally, the effectiveness of the learned router depends on the neural operator exhibiting heterogeneous performance across samples. Finally, LLMs were occasionally used to verify correctness of proofs, but all proof strategies originated from the authors and relevant literature. The code developed for the experiments was written by the authors with the help of occasional code completions. The central components (e.g., the hybrid solver implementation and the greedy-router training pipelines) were implemented exclusively by the authors. All substantive intellectual contributions, which include ideas, theorems, and analyses, are our own. The paper is 2509.24814.

--- Page 1 ---

A GREEDY PDE ROUTER FOR BLENDING NEURAL OPERATORS AND CLASSICAL METHODS Sahana Rayan Department of Statistics University of Michigan Ann Arbor, MI 48104, USA Yash Patel Department of Statistics University of Michigan Ann Arbor, MI 48104, USA Ambuj Tewari Department of Statistics University of Michigan Ann Arbor, MI 48104, USA ABSTRACT When solving PDEs, classical numerical solvers are often computationally expensive while machine learning methods can suffer from spectral bias. Designing an optimal hybrid iterative solver–where, at each iteration, a solver is selected from an ensemble of solvers to leverage their complementary strengths–poses a challenging combinatorial problem. While greedy selection is desirable for its constant-factor approximation guarantee to the optimal solution under Lipschitz assumptions, it requires knowledge of the true error at each step, which is unavailable in practice. We address this by proposing an approximate greedy router that efficiently mimics a greedy approach to solver selection. Empirical results on the Poisson and convection–diffusion equations show that our method consistently reduces final error and area-underthe-curve (AUC) of the error trajectory relative to single-solver baselines and existing hybrid approaches such as HINTS. In particular, our method reaches comparable error levels in substantially fewer iterations while exhibiting more stable error decay. 1 Introduction Natural phenomena and engineered systems are often governed by ordinary and partial differential equations (PDEs). Solving these equations enables a variety of tasks - predicting the evolution of the system through simulation (e.g., forecasting weather using the Navier-Stokes equations) [Kalnay, 2003, Bauer et al., 2015, Staniforth and Cotˆ e, 1991], addressing control problems (e.g., optimizing heat shield design for spacecrafts using heat transfer equations) [Anderson, 1989, Troltzsch, 2010], and tackling inverse problems based on external measurements (e.g., reconstructing brain activity from EEG data using electrophysiological models) [Baillet et al., 2001, Grech et al., 2008]. Traditionally, PDEs are solved using finite difference methods [Smith, 1985, LeVeque, 2007], which discretize the spatial and/or temporal domain and approximate the solution by solving a system of equations at the discrete grid points. Similarly, finite element methods [Hughes, 2003, Bathe, 2006] construct a system of equations based on local basis functions, and then rely on numerical linear algebra techniques–such as the Jacobi and GaussSeidel method [Saad, 2003]–to compute the solution. These iterative methods, however, lack generalization across different initial conditions, boundary conditions, or forcing functions, as even minor changes to any of these parameters require re-solving the entire system of equations from scratch. These challenges have motivated the development of neural operators [Kovachki et al., 2023]–a class of machine learning models that aim to learn the solution operator directly, enabling fast inference across a range of parameters. Neural operators lift the classical linear layer in neural networks to an operator–typically a kernel integral operator–between function spaces. Despite strong approximation guarantees [Chen and Chen, 1995], neural operators are prone to spectral bias [Liu and Cai, 2021, Liu et al., 2024, You et al., 2024, Khodakarami et al., 2025, Xu et al., 2025], similar to traditional neural networks [Rahaman et al., 199]. As a result, they favor learning low-frequency components of the solution function and often struggle to capture high-frequency features that iterative solvers excel at. arXiv:2509.24814v2 [stat.ME] 7 May 2026

--- Page 2 ---

Recognizing this complementarity, Zhang et al. [2024] proposed HINTS, a hybrid solver that interleaves a neural operator pass every τ th iteration of a classical solver, with τ fixed a priori. While HINTS reports significant empirical gains in convergence and accuracy over the standalone neural operator and Jacobi solver, this fixed schedule can be detrimental: a poorly timed neural correction may increase the error and undo recent progress. Although state-dependent solver schedules are desirable, this naturally suggest a reinforcement learning formulation for learning routing policies.

1. Bold header: Adaptive Solver Selection for PDEs

This system can select solvers from an ensemble of classical methods and neural operators at each iteration based on minimizing a convex surrogate loss, mimicking a greedy approach to solver selection without needing access to the true error at each step.

2. Bold header: Near-Optimal Approximation Guarantee

The learned router is shown to achieve "a constant-factor approximation to the optimal strategy for linear PDEs, provided the updates are error-reducing (Lipschitz with constant < 1) and zero-preserving–conditions under which weak supermodularity holds."

3. Bold header: Faster Convergence than Fixed Schedules

The approximate greedy solver achieves faster convergence than HINTS and individual solvers, often reaching comparable error levels in fewer iterations on the Poisson and convectiondiffusion equations, by exhibiting near-monotone decay compared to the oscillatory error behavior of fixed schedules like HINTS.

4. Bold header: Robustness Across Solver Ensembles

The system demonstrates that larger ensembles often yield improvements over both final error and AUC relative to the best pairwise configuration, indicating that additional solver diversity can be leveraged, as seen in Table 2 results where Router(NO ∪ W) outperforms Router(NO ∪ one solver) for larger sets W.

5. Bold header: Frequency-Aware Error Reduction

The routing strategy "routes to whichever solver most decreases the full L2 error, and by Parseval’s identity e(t)h squared = P m uˆ(t)m − uˆm squared, reductions in the objective correspond to reducing energy across all modes rather than giving preferential treatment to a subset. This leads to more uniform convergence across the spectrum" compared to fixed schedules.

When solving PDEs, classical numerical solvers are often computationally expensive, while machine learning methods can suffer from spectral bias, failing to capture high-frequency components. Designing an optimal hybrid iterative solver--where, at each iteration, a solver is selected from an ensemble of solvers to leverage their complementary strengths--poses a challenging combinatorial problem. While greedy selection is desirable for its constant-factor approximation guarantee to the optimal solution under Lipschitz assumptions, it requires knowledge of the true error at each step, which is unavailable in practice. We address this by proposing an approximate greedy router that efficiently mimics a greedy approach to solver selection. Empirical results on the Poisson, convection-diffusion, and anisotropic diffusion equations show that our method consistently reduces time-to-solution relative to single-solver baselines and existing hybrid approaches such as HINTS.