Locally Subspace-Informed Neural Operators for Efficient Multiscale PDE Solving
Authors: Alexander Rudikov
Vladimir Fanaskov
Sergei Stepanov
Buzheng Shan
Ekaterina Muravleva
Yalchin Efendiev
Ivan Oseledets
Venue: ICLR 2026 Poster (2026)
Keywords: Neural Operators
heterogeneous PDEs
scientific machine learning
Generalized Multiscale Finite Element Method
localized spectral basis functions
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)

Abstract: We propose GMsFEM-NO, a novel hybrid framework that combines the robustness of the Generalized Multiscale Finite Element Method (GMsFEM) with the computational speed of neural operators (NOs) to  create an efficient method for solving heterogeneous partial differential equations (PDEs). GMsFEM builds localized spectral basis functions on coarse grids, allowing it to capture important multiscale features and solve PDEs accurately with less computational effort. However, computing these basis functions is costly. While NOs offer a fast alternative by learning the solution operator directly from data, they can lack robustness. Our approach trains a NO to instantly predict the GMsFEM basis by using a novel subspace-informed loss that learns the entire relevant subspace, not just individual functions. This strategy significantly accelerates the costly offline stage of GMsFEM while retaining its foundation in rigorous numerical analysis, resulting in a solution that is both fast and reliable. On standard multiscale benchmarks—including a linear elliptic diffusion problem and the nonlinear, steady-state Richards equation—our GMsFEM-NO method achieves a reduction in solution error compared to standalone NOs and other hybrid methods. The framework demonstrates effective performance for both 2D and 3D problems. A key advantage is its discretization flexibility: the NO can be trained on a small computational grid and evaluated on a larger one with minimal loss of accuracy, ensuring easy scalability. Furthermore, the resulting solver remains independent of forcing terms, preserving the generalization capabilities of the original GMsFEM approach. Our results prove that combining NO with  GMsFEM creates a powerful new type of solver that is both fast and accurate.

Published as a conference paper at ICLR 2026
LOCALLY SUBSPACE-INFORMER NEURAL OPERATORS FOR EFFICIENT MULTISCALE PDE SOLVING
Alexander Rudikov
Vladimir Fanaskov AXXX,
AXXX,
Institute of Numerical Mathematics
Laboratory of Computational Intelligence
ay- rudi kovê 1. com
Sergei Stepanov North-Eastern Federal University
Buzheng Shan
Department of Mathematics, Texas A&M University
Ekaterina Muravleva
Yalchin Efendiev
Institute of Numerical Mathematics
Department of Mathematics, Texas A&M University
Ivan Oseledets
AXXX,
Institute of Numerical Mathematics
ABSTRACT
We propose GMsFEM-NO, a novel hybrid framework that combines the robustness of the Generalized Multiscale Finite Element Method (GMsFEM) with the compu- tational speed of neural operators (NOs) to create an efficient method for solving heterogeneous partial differential equations (PDEs). GMsFEM builds localized spectral basis functions on coarse grids, allowing it to capture important multiscale features and solve PDEs accurately with less computational effort. However, com- puting these basis functions is costly. While NOs offer a fast alternative by learning the solution operator directly from data, they can lack robustness. Our approach trains a NO to instantly predict the GMsFEM basis by using novel subspace- informed loss that learns the entire relevant subspace, not just individual functions. This strategy significantly accelerates the costly offline stage of GMsFEM while retaining its foundation in rigorous numerical analysis, resulting in a solution that is both fast and reliable. On standard multiscale benchmarks--inhhun elliptic diffusion problem and the nonlinear, steady-state Richards equation- our GM method achieves a reduction in solution error compared stan- dalone NOs and other hybrid methods. The framework demonstrates effective performance for both 2D and 3D problems. A key advantage its discretization flexibility: the NO can be trained on small computational grid and evaluated on larger one with minimal loss accuracy, ensuring easy scalability. Furthermore, the resulting solver remains independent of forcing terms, preserving the gener- alization capabilities the original GMsFEM approach. Our results prove that combining NO with GMsFEM creates powerful new type of solver that is both fast and accurate.
a linear
INTRODUCTION
Many practical multiscale problems involve highly heterogeneous properties with high-contrast variations across multiple scales, posing significant challenges for the numerical solution of partial differential equations (PDEs). A well-established approach for such problems the Generalized
Multiscale Finite Element Method (GMsFEM) Efendiev et al. (201 2013); Chung et al. (2016),

---

Published as a conference paper at ICLR 2026
which constructs localized spectral basis functions coarse grids. By solving local eigenproblems, GMsFEM captures fine-scale information, enabling accurate coarse- scale solutions. However, this accuracy comes at a high computational cost due to the expense of solving these local cigenproblems.
Recently, data-driven solvers. particularly neural operators (NOs) like Fourier Neural Operators (FNOs) Li al. (020); Kovachki et al. (2023) ): |anaskov & Oseledets (023): Tran al. (2021) and DeepONets Lu et al. (2021): Wang et al. (2021), have emerged as powerful alternatiee for al. (2024); Karniadakis al. (2021). While
accelerating PDE simulations Azizzadenesheli
effective for problems with smooth coefficients, standard NOs often struggle to efficiently capture
the localized features of high-contrast heterogeneities, typically requiring extensive data and large network architectures.
In this work, we introduce GMsFEM-NO, a hybrid framework that combines the robustness of GMsFEM with the speed of neural operators. Our key innovation is subspace-informed NO that learns to map a heterogeneous coefficient field directly to the low-dimensional subspace spanned by the GMsFEM basis functions. Instead of learning individual basis canctions-white can be sensitive to small perturbations- design a novel subspace-aware loss function that enforces physical consistency at the subspace level. This approach offers several advantages: it is more data-efficient than learning the full PDE solution, as the basis functions are smoother and of lower dimension; and it is more robust than a pure NO, the final solution is obtained through GMsFEM, ensuring legitimacy even with imperfect basis predictions.
Our approach is distinct from existing hybrid methods Bhattacharya al. 2024 Vasilyeva et al. (2020): Wang et al. (2020); Liu et al. (2023); ): |röpf et al. (202) that combine machine learning with numerical
Those methods typically assume known macroscopic equation form and learn effective coefficients, which is infeasible for problems without scale separation and with high contrast. In contrast, GMsFEM-NO learns the macroscopic solution space itself, in the form of multiscale basis functions. making it suitable for these more challenging settings. A related approach (Spridonov al. (2025) used a fully connected neural network to predict an additional basis function for the steady-state Richards equation Richards (1931) Farthing & Ogden] (2017), supplementing an existing set of precomputed basis functions. While this approach enhanced prediction accuracy, it failed deliver computational efficiency gains because traditional methods still generated most basis functions. Furthermore, the simplicity of the fully connected architecture limited its ability to account for spatial variations, potentially compromising prediction accuracy for high-contrast data. Another category of related work aims to reduce the computational cost of PDE solving via reduced-order modeling (POD Volkwein (2013), DeepPOI) Franco et al.) and PCANet (Bhattacharya et al. (202T). DeepPOD and PCANet also leverage neural networks to learn compact solution representations, providing relevant baseline for comparing the efficiency of our method.
We validate GMsFEM-NO on two challenging benchmarks with high-contrast coefficients: linear elliptic diffusion problem and the nonlinear steady-state Richards equation. Results showed that our approach is better than NO in terms of solution accuracy and requires less training data achieve similar accuracy. Additionally, it reduces cosis-construction time by more than 60 times compared to traditional GMsFEM.
Our main contributions are:
1. We introduce novel hybrid approach (GMsFEM-NO) that combines the strengths of NOs with GMsFEM (see Fig. D.
2. A new subspace-informet loss function for learning stable and generalizable solution subspaces.
3. The approach is evaluated on high-contrast PDEs and shown deliver the same results as GMsFEM at a fraction of the computational cost.
4. Demonstration of resolution invariance of GMsFEM-NO: effective training on low- resolution data for application to high-resolution problems.
5. Superior in-distribution and diitr---ittiibttinn performance compared to standard NOs, without requiring domain adaptation.
2

---

Published as a conference paper at ICLR 2026
Nr
Input dataset,
K
(a)
Spectral problem
L
(aa }
Nor
NO
y=1
R
(b)
NO
Mus j=1
Solve the coarse-scale system: Aguo - So. Ao RAR So R
Figure 1: Illustration of training (a) and inference (b) stages of the proposed GMsFEM-NO method. NO is trained on heterogeneous fields n that defined on subdomain wi to predict subspace of basis functions (vlg)) 1ー=1, where Nar the number of basis functions. During training the subspace- NNx informed loss C is applied align predicted subspace {w5 Jj=1 with {v Jo=1 During inference stage (b), the predicted subspace forms the matrix R that projects matrix A and vectors the coarse space.
LOCALLY SUBSPACE-INFORMED NEURAL OPERATORS
2.1 DIFFUSION EQUATION
We consider the diffusion equation with heterogeneous coefficient
(r())Vu((z)) = f(x), I E Ω = (0, 1)P, u(x) = 0,
(I)
where A2 is a boundary of the unit hypercube Q, and x(xr) is heterogeneous field with high contrast. In partieular, assume that x(x) > 0, while x(x) have very large variations. For example, in this work we use k(x)) with values in range [1, 9600].
2.2 STEADY-STAR STATE RICHARDS' EQUATION
The steady-state version of Richards' equation, which describes water movement in unsaturated porous media, takes the following form:
(n(2,, (())) f(x), E = (0, 1)", u(x)) weon = 0,
(2)
where x <(z, u(x)) is unsaturated hydraulic conductivity, u(x) is the water pressure and f(x) is source or sink term.
We consider the Haverkamp model Haverkamp et al. (1977) to define sa(,, u(x)):
1
u(x)) K2(x) K,(u(x)) <(x)
+ |u]
where x(z) is heterogencous field with high contrast that denotes the permeability of soils, K+(u) represents the relative hydraulic conductivity, K2(x) stands for the saturated hydraulic conductivity,
1
All the multiscale heterogeneity is incorporated in s<x)) without regard to u, and the non-linearity.
includes all
1 + [u]

---

Published as a conference paper at ICLR 2026
2.3
GENERALIZED MULTISCALE FINITE ELEMENT METHOD
2.3.1 MULTISCALE SPACE APPROXIMATION
Multiscale methods Efendiev & Hou (009) form a broad class of numerical techniques. They are based on constructing multiscale basis functions local domains to capture fine-seale behavior.
Let TH be a coarse mesh of the domain Ω C R (with = 2 or 3), such that Th = U K,, where each K; is coarse cell and No is the number of coarse cells. Let Tn. be a fine grid obtained by refinement of Th, with h < H. We denote by {xx )= the nodes of the coarse mesh Th. where Ne the number of nodes of the coarse mesh. Let wi be the subdomain defined as the collection of coarse cells containing the coarse grid node 2; (see Fig. i Appendix A):
U U
K, TH 24 E {)
To ensure accurate approximations on the coarse mesh Th, we construct spectral multiscale basis functions following the GMsFEM. GMsFEM contains two stages:
Offline stage:
1. Coarse and Local Domain Definition: Define the coarse grid TH and generate the associ- ated local domains wi for i = 1....., N,,
2. Local Spectral Problem Solving: In each local domain W,, solve local spectral problem
to obtain a set of eigenvectors { j/v 1==1" where N is the number of coarse eigenvectors.
3. Multiscale Basis Function Construction: Select the first Nbr eigenvectors from each W; and multiply them by partition unity function Xi Babuska & Lipton (201T); [aabška N
et al. 2008); Strouboulis et al. 2000 to create the final multiscale basis functions (v)) J==1"
where Nbr N.
4. Global System Assembly: Map the local degrees of freedom to global and form a restriction matrix R.
Online stage:
1. Projection:
Use R to project the fine-scale system onto the coarse space.
2. Solution: Compute the solution within the coarse multiscale space.
3. Reconstruction: Obtain the fine-scale approximation by applying the prolongation operator
R to the coarse-scale solution.
2.3.2 SPECTRAL PROBLEM
We denote by (S2 the usual finite element discretization of piecewise linear continuous functions with respect to the fine grid Th. For each local domain wi, we define the Neumann matrix A by
T A wn f V(E) Vwn dr,
Von, Wy 4 (ws,
Uh
and the Mass matrix S by
Up S.h wh = (a(x)orw, da, Von, wn (wi).
We consider the finite dimensional symmetric eigenvalue problem A A.S
and denote its eigenvalues and eigenvectors by {At)) j=1 and (Φ) N J== respectively. Note that A O corresponds to the constant eigenvector 60 const. We order cigenvalues as
{
A
The cigenvectors (6))) N J==1 form S -orthonormal basis of 4 (w:).
4

---

Published as a conference paper at ICLR 2026
2.3.3 SOLVING OF THE COARSE-SCALE SYSTEM
For each local domain wi, we select eigenvectors corresponding the Nor N smallest eigenvalues and define a multiscale subspace
span 1/th Xedly
Nos,
No}
and define the restriction matrix R 1 ..... W
u N Nai
Coarse-grit solution
the finite element projection of the fine-scale solution into the space solution uo is given by
More precisely, multiscale
Aguio = fo.
where Ag RAR the projected system matrix, fo R is projected right-hand side. The reconstructed fine-scale solution is u = R to.
2.4 NEURAL OPERATOR
Here we consider one type of NO that employs Fourier modes, but there are no restrictions on using other types of NOs. Fourier neural operators (FNOs) class of NOs motivated by Fourier spectral methods. Originally, Li et al. 2020) formulate each operator layer as
(z = ( b(t) 1( (2(0))
(4)
where ((1 + b(0 is an affine point-wise map,
1CO) (= = IFFT (F (R(O.FFT(2))
a kernel integral operator. The Fourier domain weight matrices (R J==1 require O(LH²M
parameters, where H is the hidden size, M is the number of the top Fourier modes that are kept, and D is the dimension of the problem.
In Factorised FNO (F-FNO) Tran al. 2021). the operator layer in (4) is changed
C²
(0)
w (T) (2(0) + b b/
where (x(00) EAED [FFT ( (=)) In this case, the number of parameters MD. Therefore, the FFNO reduces model complexity and scales efficiently deeper networks.
2.5
PROPOSED METHOD
2.5.1
GMsFEM-NO ALGORITHM
We propose an efficient hybrid method for generating basis functions in the GMsFEM using NOs, significantly accelerating the offline stage.
Local domains vary in shape and orientation (see Appendix A where orientation refers to the relative placement of the coarse node x; shared by all cells the local domain. We address this variability by categorizing the local domains based on their geometry: into full, half, and corner types in 2D, and into full, half, quarter, and corner types in 3D (see Appendix AA). Before training, we normalize the orientation of each local domain by rotating both the input data and the target basis functions, ensuring a standardized coarse node z; position within each group. This preprocessing step guarantees consistency in the input structure for the NO.
We train separate NOs, each specialized for one domain group (see Appendix B). Each NO predicts Nor basis functions for local domains within its assigned category. This group-specific approach improves prediction accuracy by accounting for geometric variations across local domain types.
For test data, we first decompose the computational domain into local domains. The corresponding NO then generates the required basis functions. The predicted basis functions are extended to the domain Ω (with zeros padded outside their respective local domains) and vectorized to construct the
5

---

Published as a conference paper at ICLR 2026
restriction matrix R. Finally, the online stage of GMsFEM is executed to compute the multiscale solution.
This approach substantially reduces offline computational costs while maintaining the accuracy and flexibility of GMsFEM, making it particularly suitable for problems with heterogeneous or highly varying coefficients.
2.5.2 SUBSPACE-INFORMNN LOSS FUNCTIONS
The selection of an appropriate loss function critical when training NOs. propose a Subspace Alignment Loss (SAL) that directly optimizes the geometric consistency of the learned subspaces. Let R .i...... Vin represent the target subspace basis and It' denote the predicted subspace. The SAL measures alignment between subspaces using their orthonormalized bases Qz and Q jr:
SSAL = E2 Nox - 1 Q 11]
(5)
where the Frobenius norm term 11 Q Rc 11 quantifies the subspace overlap, achieving its maximum value No when subspaces are perfectly aligned. Appendix[[provee SAL's theoretical foundation and provides error bounds connecting subspace alignment to solution accuracy.
While SAL ensures subspace coherence, it may overlook finer discrepancies in how functions are projected onto the subspaces. To enforce consistency in projection behavior, we introduce Projection Regularization term. This term evaluates the discrepancy between projections of randomized test vector o² onto the target and predicted subspaces, governed by their projection matrices Pn: and Pp::
SSAL LSAL + x R,, (PR -Pi.) 1l122 u N~,(0,
where v = Ek cuy PR. = Qx: ni: Pp. = -Q Q and À is hyperparameter.
We compare proposed loss functions (5) 66) with conventional one which is L2 loss. Since basis functions are defined only up their sign (see Appendix [E), the conventional Lg loss is adapted to account for this invariance, resulting in the Relative Basis Function L2 Loss (RBFL2):
[
1 + +011 ]
LRBFL = Esg min
where vy and و denote the j-th target and predicted basis functions for the i-th local domain wi The minimization over - ensures invariance to sign permutations.
3 RESULTS
We use datasets of 2D coefficients at resolutions of 1002 and 250 and 3D coefficients at 50 and 1003 (see example in Appendix F). The domain is partitioned into N, subdomains corresponding
to coarse grids (e.g., N, = 36 for 5 x 5, 121 for 10 x 10, 216 for x 5 x 5, 729 for 8 x 8 x 8). For grid sizes except 100 x 100 x 100, we used 1000 samples (800 train, 200 test). For the 100 x 100 x 100 grid, we used 150 train and 50 test samples due to computational constraints. The different local domain types occur with varying frequencies within a single sample (see Appendix For training NOs, we utilize the first 8 basis functions (Nar per subdomain as training targets.
To evaluate method robustness, we consider two right- side configurations:
Uniform unit forcing term
( 1.
(8)
Spatially variable forcing (see Appendix F defined by
f(x) ~ NN (I- -))
(9)

---

Published as a conference paper at ICLR 2026
To measure quality of the obtained solutions fine grid, the following metrics:
So 2 u//l 17]2dd
dx
Jo
Vuy
Vil
V
dx
L2 = En
H1 E.,
11 dee
Sn
Vu
dx
All experiments performed single Nvidia Tesla H100 80Gb HBM3. The comparison of our approach with baseline methods is presented in Appendix G
3.1 LRAFLY VS. LSAL. LRAL.
To determine the optimal training configuration for the NOs predicting basis function subspaces, we performed grid search over architectural parameters and training hyperparameters. Full specifica- tions (except loss function type) are in Appendix H Loss function results for Nor = 8 are shown in For Nor 4, results are in Appendix
Table I
As shown in Table CRBFL, underperforms compared to our proposed subspace alignment losses (CSaL. LSAL-PR). For the Richards equation with simple right- side ( and Nor 8. our proposed loss improves the relative L2 metric by a factor of 1.8. Notably, the projection regularization term in CSAL.PR yielded nearly identical results to LSAL. While projection regularization had a minimal impact on smaller grids-likely because the subspace alignment term alone suffices-it effect became significant for larger problems. For the 2502 grid using Richards' equation with right-hand side it reduced the L2 error from 1.82% to 1.72% (see Table 3)
Table I: Performance comparison of loss functions for NO training (100 x 100 grid, No 36).
LRBFL
LSAL
LSAL PR
Nbr Dataset
L2
H1
L2
H1
L2
H1
Diffusion,
1.75%
14.83%
1.06%
11.57%
1.06%
11.65%
Diffusion,
3.53%
21.77%
2.82%
19.07%
2.81%
19.03%
Richards,
3.46%
15.04%
1.88%
11.10%
1.87%
11.25%
Richards,
3.77%
22.38%
2.99%
19.61%
2.99%
19.60%
3.2 GMsFEM vs. GMSFEM-NO
In this section, we compare the performance of the original GMsFEM and our proposed GMsFEM NO methods in terms of solution accuracy (quantified by La and Hi metrics) and computational efficiency for basis functions generation. For large-scale 3D simulations (100 x 100 x 100 grid), we employed the GMsFEM-NO method with the SAL loss. Although the SAL-PR loss offers benefits for smaller problems (see Appendix D) its computational cost becomes prohibitive at this scale, making the standard SAL loss the practical choice.
As shown in Tables GMsFEM across all datasets and grid sizes. While GMsFEM-NO shows slightly better results for some configurations, this is likely due to statistical variation. The same behaviour is observed the time-dependent equation and the diffusion equation with mixed boundary conditions (see Appendix and M); these experiments were conducted solely on the 250 x 250 grid.
and 5 GMsFEM-NO achieves nearly identical L2 and Hy errors to
Table 6 compares the time required to generate basis functions using the GMsFEM offline stage and GMsFEM-NO for different grid sizes and N, values. GMsFEM-NO employs several NOs, one for each local domain type. The proposed method achieves more than 60 speedup, demonstrating its computational superiority. Basis calculation speedup grows with grid size and dimensionality.
While GMsFEM-NO reduces the cost of the traditional GMsFEM offline stage over many subsequent simulations, it does not eliminate cost of train data generation and training. The breakeven point, defined as the number of inference samples required to offset these initial costs, is derived in the Appendix N
7

---

Published as a conference paper at ICLR 2026
Table Performance comparison of GMsFEM and GMsFEM-NO for 2D (100 100, NN =6))
GMsFEM
GMsFEM-NO
Not
Dataset
L2
H1
L2
H1
Diffusion,
1.15%
11.68%
1.06%
11.57%
Diffusion,
2.82%
19.07%
2.81%
19.03%
Richards,
2.03%
11.68%
1.87%
11.25%
Richards,
3.09%
20.20%
2.99%
19.60%
Table 3: Performance comparison GMsFEM and GMsFEM-NO for 2D (250 250, N, 121).
GMsFEM
GMsFEM-NO. LSAL
GMsFEM-NO. LSAL
Nbt
Dataset
L2
H1
L2
H1
L2
H1
Diffusion,
1.12%
14.01%
1.13%
13.92%
1.05%
14.07%
8 Diffusion,
1.60%
21.49%
1.62%
22.50%
1.60%
22.03%
Richards, 8
1.79%
14.57%
1.82%
14.27%
1.72%
14.53%
Richards,
1.62%
21.76%
1.62%
22.62%
1.60%
22.13%
Table 4: Performance comparison of GMsFEM and GMsFEM-NO for 3D (50 x 50 50, N, 216).
GMsFEM
GMsFEM-NO. LSAL
GMsFEM-NO, LSAL
Not
Dataset
L2
H1
L2
H1
L2
H1
Diffusion,
3.07%
20.72%
3.10%
20.39%
3.10%
20.39%
8 Diffusion,
5.08%
25.26%
5.01%
24.92%
5.00%
24.86%
Richards,
4.04%
15.43%
4.12%
15.37%
4.14%
15.39%
Richards,
5.02%
24.89%
5.02%
24.92%
5.00%
24.89%
Table 5: Performance comparison of GMsFEM and GMsFEM-NO for 3D (100 100 X 100,
N 729).
GMsFEM
GMsFEM-NO, LSAL
N Dataset
L2
H1
L
H1
Diffusion,
1.62%
14.76%
1.68%
14.98%
8 Diffusion, 9
3.0%
12.55%
3.04%
12.84%
Richards,
2.54%
17.83%
2.57%
17.91%
Richards,
2.55%
17.85%
2.57%
17.91%
Table 6: Basis generation time: GMsFEM-NO vs. standard GMsFEM offline stage.
Grid
N₂
GMsFEM, sec.
GMsFEM-NO, sec.
100 x 100
36
16.87
0.28
250 250
121
210.5
0.31
50 x 50 x 50
216
935.4
0.84
100 100 100
729
10547.2
1.33
8

---

Published as a conference paper at ICLR 2026
3.3 STANDALONE NOS vs. GMsFEM-NO
In this section, we compare GMsFEM-NO against standalone SOTA NOs, including F-FNO, GNOT HHoo et al. (2023) and Transolver++ Luo al. (2255. While these standalone models offer fast inference, their accuracy deteriorates significantly on high-contrast datasets. As shown in Table 7 (250x250 grid), GM>FEM-NG achieves superior accuracy; for Richards' equation with complex source term ( it reduces the relative L2 error by 2.8x compared to the best standalone NO. Full training details are in Appendix !
A critical advantage of GMsFEM-NO over the standalone NO lies in its independence from the right-hand side terms of the PDE. The standalone NO exhibits catastrophic failure when tested on out-of-distribution forcing terms, as evidenced in Table 13 in Appendix K
Since each coefficient contains multiple local domains of each type, GMsFEM-NO requires fewer samples than F-FNO for training. As shown in Table 8 when Nuain is reduced below 800, the error for F-FNO begins to increase significantly. In contrast, GMsFEM-NO'S accuracy remains stable across the range of 800 to 400 samples. Even with only 200 samples, the performance degradation for GMsFEM-NO remains small; for example, on Richards' equation with the simple right-hand side (8) the error increases only modestly from 1.85% to 2.07%.
Table 7: Performance comparison of NOs and GM>FEM-NO (250 250 grid)
No1 Dataset
F-FNO
GNOT
Transolver++
GMsFEM-NO
Diffusion, 1.02%
1.26%
1.15%
1.05%
8 Diffusion,
4.51%
14.29%
6.63%
1.60%
Richards,
2.44%
2.34%
2.17%
1.72%
Richards,
4.45%
14.69%
8.82%
1.60%
Table 8: Comparison of F-FNO and GMsFEM-NO performance across different training dataset sizes for 250 x 250.
D train
Diffusion,
Diffusion, g
Richards,
Richards, 9
GMsFEM-NO
1.33%
1.77%
2.07%
1.77%
200
F-FNO
2.85%
11.56%
6.52%
11.49%
GMsFEM-NO
1.15%
1.63%
1.85%
1.62%
400
F-FNO 1.60%
8.16%
4.17%
8.41%
GMsFEM-NO
1.12%
1.61%
1.78%
1.61%
600
F-FNO
1.21%
4.92%
3.27%
5.43%
GMsFEM-NO
1.13%
1.62%
1.82%
1.62%
800
F-FNO
1.02%
4.51%
2.44%
4.45%
3.4 GMsFEM-NO FOR DIFFERENT GRIDS
Table demonstrates the resolution invariance of GMsFEM-NO by training the model on coarse grid and testing it finer grid (500 500), with results compared against the standard GMsFEM solution computed directly on the fine grid. We used a 10 x 10 coarse grid (121 subdomains) for all experiments. The results demonstrate the stability of the proposed method. GMsFEM-NO performs effectively when evaluated on a grid resolution higher than its training resolution, a key advantage enabled by the neural operator's ability to generalize to different discretizations.

---

Published as a conference paper at ICLR 2026
Table 9: Evaluation of GMsFEM-NO trained on coarse grid and tested on finer grid, with comparison to standard GMsFEM.
Diffusion,
8
Diffusion,
Richards,
Richards,
Train grid
Test grid
GMsFEM-NO
100
500
2.42%
2.97%
4.70%
3.49%
250
500
1.45%
1.79%
2.25%
1.97%
GMsFEM
500
1.17%
1.46%
1.93%
1.66%
4 CONCLUSION
In this work, we propose GMsPEM-NO, novel method for solving multiscale PDEs that employs NOs to predict the multiscale basis function subspaces in the GMsFEM offline stage, replacing the conventional solution of local eigenvalue problems. validated the method on standard 2D and 3D benchmarks: linear elliptic diffusion problem and the nonlinear steady- Richards equation. Additionally, we demonstrated its efficacy for time-dependent equations and problems with mixed boundary conditions in 2D. GMsFEM-NO achieves more than 60x speedup in basis generation compared to standard GMsFEM.
A key contribution is novel subspace alignment loss function, which enables direct learning of the basis function subspace and improves the L2 accuracy over conventional CREFL, loss. The GMsFEM- NO framework remains independent of the PDE's right-hand side, allowing maintain consistent performance across varying forcing terms. This contrasts with standalone NOs, which exhibit errors exceeding 100% on out-of-distribution data. Furthermore, GMsFEM-NO demonstrates greater data efficiency, requiring half the training samples of a comparable NO. A significant ad advantage is the method's discretization invariance: GMsFEM-NO performs effectively when evaluated on grid resolutions higher than those used for training. demonstrating strong generalization across different computational meshes. By preserving the mathematical structure multiscale methods while leveraging NO speed, this work establishes practical paradigm for heterogeneous PDE simulation.
The primary limitation of our method is current restriction to structured grids due the chosen NO architecture. Additionally, our experiments focused relatively small grid sizes, which may not fully represent large-scale applications. While we successfully tested our approach on dimendependent equations and problems with mixed boundary conditions in 2D, the study primarily focused on steady-state problems with Dirichlet boundary conditions.
Future work will expand this framework in several key directions. First, we will target complex PDEs. Second, extending the framework to irregular domains is a critical next step. This is well- supported theoretically, as the GMsFEM methodology established for unstructured meshes, and can be integrated with geometry-aware NOs like Transolver++. Alongside these goals, we will also investigate performance on finer grid resolutions and the impact of coarse-grid sizing to fully realize the method's 's potential for large-scale, real-worid simulations.
REFERENCES
Jørg E Aames and Yalchin Efendiev. Mixed multiscale finite element methods for stochastic porous media flows. SIAM Journal on Scientific Computing, 30(5):3319-239, 2008
Kamyar Azizzadenesheli, Nikola Kovachki, Zongyi i Li, Miguel Liu-Schiaffini, Jean Kossaifi, and Anima Anandkumar. Neural operators for accelerating scientific simulations and design. Nature Reviews Physics, 6(5):320-328, 2024.
Ivo Babuska and Robert Lipton. Optimal local approximation spaces for generalized finite element methods with application to multiscale problems. Multiscale Modeling & Simulation, 9(1):333-44 2011.
10

---

Published as a conference paper at ICLR 2026
Ivo Babuška, Victor Nistor, and Nicolae Tarfulea. Generalized finite element method for second-order elliptic operators with dirichlet boundary conditions. Journal of Computational and Applied Mathematics, 218(1):175-18) 2008
Pau Batlle, Matthieu Darey, Bamdad Hosseini, and Houman Owhadi. Kernel methods are competitive for operator learning. Journal of Computational Physics, 496: 12549, 2024.
Thomas Bendokat, Ralf Zimmermann, and P-A Absil. A grassmann manifold handbook: Basic geometry and computational aspects. Advances in Computational Mathematics, 50(1);6, 2024.
Kaushik Bhattacharya, Bamdad Hosseini, Nikola B Kovachki, and Andrew M Stuart. Model reduction and neural networks for parametric pdes. The SMAI journal of computational mathematics, 121-157, 2021.
Kaushik Bhattacharya, Nikola B Kovachki, Aakila Rajan, Andrew M Stuart, and Margaret Trautner. Learning homogenization for elliptic operators. SIAM Journal on Numerical Analysis, 62(4): 1844 1873, 2024.
Xiangning Chen, Chen Liang, Da Huang, Esteban Real, Kaiyuan Wang, Hicu Pham, Xuanyi Dong, Thang Luong, Cho-Jui Hsieh, Yifeng Lu, et al. Symbolic discovery of optimization algorithms. Advances in neural information processing systems, 36:49205 49233, 2023.
Eric Chung, Yalchin Efendiev, Guanglian Li, and Maria Vasilyeva. Generalized multiscale finite element methods for problems in perforated heterogeneous domains. Applicable Analysis, 95(10): 2259-2279, 2016.
DeepMind, Igor Babuschkin, Kate Baumli, Alison Bell, Surya Bhupatiraju. Jake Bruce, Peter Buchlovsky, David Budden, Trevor Cai, Aidan Clark, Ivo Danihelka, Antoine Dedieu, Clau- dio Fantacci, Jonathan Godwin, Chris Jones, Ross Hemsley, Tom Hennigan, Matteo Hessel, Shaobo Hou, Steven Kapturowski, Thomas Keck, Iurii Kemaev, Michael King, Markus Kunesch, Lena Martens, Hamza Merzic, Vladimir Mikulik, Tamara Norman, George Papamakarios, John Quan, Roman Ring, Francisco Ruiz, Alvaro Sanchez, Laurent Sartran, Rosalia Schneider, Eren Sezener, Stephen Spencer, Srivatsan Srinivasan, Miloš Stanojević, Wojeiech Stokowiec, Luyu Wang, Guangyao Zhou, and Fabio Viola. The DeepMind JAX Ecosystem, 2020. URL http://github. com wooggoog//oooooooooooooooooooooo goog e-deepmi nd
Yalchin Efendiev and Thomas Y Hou. Multiscale finite element methods: theory and applications, volume 4. Springer Science & Business Media, 2009.
Yalchin Efendiev, Juan Galvis, and Kiao-Hui Wu. Multiscale finite element methods for high high-contrast problems using local spectral basis functions. Journal of Computational Physics, 230(4))937 2011.
Yalchin Efendiev, Juan Galvis, and Thomas Y Hou. Generalized multiscale finite element methods (gmsfem). Journal of computational physics, 251:116-135, 2013.
Vladimir Sergeevich Fanaskov and Ivan V Oseledets. Spectral neural operators. In Doklady Mathematics, volume 108, pp. $226-S232. Springer, 2023.
Matthew W Farthing and Fred L Ogden. Numerical solution of richards' equation: A review of advances and challenges. Soil Science Society of America Journal, (166)1127-1126, 2017.
Nicola Rares Franco, Andrea Manzoni, Paolo Zunino, and Jan S Hesthaven. Deep orthogonal decomposition: a continuously adaptive neural network approach to model order reduction of parametrized partial differential equations.
Zhongkai Hao, Zhengyi Wang, Hang Su, Chengyang Ying, Yinpeng Dong, Songming Liu, Ze Cheng, Jian Song, and Jun Zhu. Gnot: general neural operator transformer for operator learning. In International Conference Machine Learning, pp. 2556611266 PMLR, 2023.
Roland Haverkamp, Michel Vauclin, Jaoudat Touma, PJ Wierenga, and Georges Vachaud. A comparison of numerical simulation models for one-dimensional infiltration. Soil Science Society of America Journal, 41 (1(2)285-299) 1977.
11

---

Published as a conference paper at ICLR 2026
Jan S Hesthaven and Stefano Ubbiali. Non-intrusive reduced order modeling of nonlinear problems using neural networks. Journal of Computational Physics, 363:55-78, 2018
George Em Karniadakis, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. Physics-informed machine learning. Nature Reviews Physics, 3(6):422-440, 2021.
Patrick Kidger and Cristian Garcia. Equinox: neural networks in JAX via callable PyTrees and filtered transformations. Differentiable Programming workshop Neural Information Processing Systems 2021, 2021.
Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Learning maps between function spaces with applications pdes. Journal of Machine Learning Research, 24(89):1-97, 2023.
Fabian Kröpfl, Roland Maier, and Daniel Peterseim. Operator compression with deep neural networks. Advances in Continuous and Discrete Models, 2022(1):29, 2022.
Fabian Kröpfl, Daniel Peterseim, and Elisabeth Ullmann. Neural network localized orthogonal decomposition for numerical homogenization of diffusion operators with random coefficients. arXiv preprint arXiv:2509. 12896, 2025.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametrie partial differential equations. arXiv preprint (0t::2000008899 2020.
Burigede Liu, Eric Ocegueda, Margaret Trautner, Andrew M Stuart, and Kaushik Bhattacharya. Learning macroscopic internal variables and history dependence from microscopic models. Journal of the Mechanics and Physies of Solids, 178: 178:105329, 2023.
Ilya Loshchilov and Frank Hutter. Decoupled weight decay regularization. arXiv preprint arXiv:111.05101, 2017.
Lu Lu, Pengzhan Jin, and George Em Karniadakis. Deeponet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators. arXiv preprint arXiv: 0900093/1990000 2019.
Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, and George Em Karniadakis. Learning nonlinear operators via deeponet based on the universal approximation theorem of operators. Nature machine intelligence, 3(3):218-229) 2021.
Huakun Luo, Haixu Wu, Hang Zhou, Lanxiang Xing. Yichen Di, Jianmin Wang, and Mingsheng Long. Transolver++: An accurate neural solver for pdes on million-scale geometries. arXiv preprint wwXiv:2502.02414, 2025.
André LG Mandolesi. Asymmetric geometry of total grassmannians. arXiv preprint arXiv:2310. 17865, 2023.
Brek Meuris, Saad Qadeer, and Panos Stinis. Machine-learning custom-made basis functions for partial differential equations. arXiv preprint arXiv:2111.05307, 2021.
Brek Meuris, Saad Qadeer, and Panos Stinis. Bechbine-learing-bbee based spectral methods for partial differential equations. Scientific Reports, 13(1):1739, 2023.
Lorenzo Adolph Richards. Capillary conduction of liquids through porous mediums. physics, 1(5): 318-333, 1931.
Denis Spiridonov, Sergei Stepanov, and Tina Mai. Prediction of discretization of online gmsfem using deep learning for richards equation. Journal of Computational and Applied Mathematics, 454:116167, 2025.
Theofanis Strouboulis, Ivo Babuška, and Kevin Copps. The design and analysis of the generalized finite element method. Computer methods in applied mechanics and engineering, 181(1-3):43-66 69, 2000.
12

---

Published as a conference paper at ICLR 2026
Alasdair Tran, Alexander Mathews, Lexing Xie, and Cheng Soon Ong. Factorized fourier neural operators. arXiv preprint arXiv:2111.13800 2021.
Maria Vasilyeva, Wing T Leung, Eric T Chung, Yalchin Efendiev, and Mary Wheeler. Learning macroscopic parameters in nonlinear multiscale simulations using nonlocal multicontinua upscaling techniques. Journal of Compulational Physics, 412:10932, 2020.
Maria Vasilyeva, Aleksei Tyrylgin, Donald L Brown, and Anirban Mondal. Preconditioning markov chain monte carlo method for geomechanical subsidence using multiscale method and machine learning technique. Journal of Compulattonnal and Applied Mathematicx, 392:1 13420, 2021.
Stefan Volkwein. Proper orthogonal decomposition: Theory and reduced-order modelling. Lecture
Notes, University of Konstanz, 4(4):1 2013.
Min Wang. Siu Wun Cheung, Wing Tat Leung, Eric T Chung, Yalchin Efendiev, and Mary Wheeler. Reduced-order deep learning for flow dynamics. the interplay between deep learning and model reduction. Journal of Compulational Physics, 401:108939, 2020.
Sifan Wang, Hanwen Wang, and Paris Perdikaris. Learning the solution operator of parametric partial differential equations with physics-informet deeponets. Science advances, 7(40):eabi8605, 2021.
Eugene Wong. Stochastic processes in information and dynamical systems. (No Tirle), 1971.
A COARSE GRID
The notation refers to the i-th local domain, where the index corresponds to the numbering of points on the coarse grid. Fig. shows examples local domains wo, W20. and W31, representing the full. half, and corner types in 2D. In 3D, there are four types: full (8 cells), half (4 cells), quarter (2 cells), and corner (1 cell), where the cell cube. Each local domain is discretized with a fine grid.
H (Coarse Grid)
W34
W20
w2
Figure 2: Illustration of 5 x coarse grid Th showing local domains of different types: the corner type wp (1 cell), half type Wa4 (2 cells), and full type who (4 cells), where the cell is square.
B TRAINING GMSFEM-NO
We train separate specialized NOs for each geometric domain type: three for 2D problems (full, half, corner) and four for 3D problems (full, half, quarter, corner), as illustrated for the 2D case in Fig. 3 (a-c). Each NO predicts the Nor basis functions for all local domains of its assigned type.
The number of local domains for each geometric type can be calculated based on the coarse grid dimensions. For 2D grid with 36 domains (5 5 cells), the counts are: 16 full, 16 half, and 4 corner
13

---

Published as a conference paper at ICLR 2026
domains. For finer 2D grid with 121 domains (10 x 10 cells), the counts are: 81 full, 36 half, and corner domains. 3D, for a grid with 216 domains (5 x 5 x cells). the distribution is: 64 full, 96 half, 48 quarter, and corner domains. In 3D, for finer grid with 729 domains (8 8 x cells), the distribution is: 343 full, 294 half, 84 quarter, and 8 corner domains.
Input dataset, k
Multiscale basis functions, 3
NO
(a)
NO
(b)
NO3
(c)
Figure 3: Multiscale basis generation algorithm for three subdomain types wyt (a) full, (b) half, (c)
corner using dedicated NOs per type with further extension Q.
C SUBSPACE ALIGNMENT Loss (SAL)
To understand the relationship between the proposed LSAL (5) and classical Grassmannian geometry
Bendokat et al. (2024); : Mandolesi] (2023), we begin with the orthogonal projection matrices. For subspace R spanned by a set of basis vectors, we compute an orthonormal basis QR via the thin QR decomposition. The orthogonal projection matrix onto R is then given by PR = QRQ
The Grassmannian distance between two d-dimensional subspaces R and R is defined using these projection matrices. The distance derivation proceeds as follows:
Pn PRll2= (r(Pp) 2tr (PP ir(PR) dim R) + dim( R?) 2tr (e(P
The matrix QhQR contains the cosines of the principal angles between the subspaces. Therefore,
PR Pa 1 2
Consequently, the Grassmannian distance simplifies to:
d(R, R
1 22
PR
Pall = k 12
-
This derivation confirms that minimizing LSAL is equivalent to minimizing the expected Grassmannian distance between the true and predicted subspaces.
In considering the theoretical soundness of our approach, we note that formal upper bound linking principal angles the final error would further strengthen the framework. When the SAL loss is small, it ensures that the learned solution remains close to the exact solution. To formalize this, let Q nc represent the multiscale solution computed using GMsFEM and Q denote the learned
14

---

Published as a conference paper at ICLR 2026
solution. The error then be bounded as follows:
||u Qadl² 2 I Iu Qnell 2 + 1Q Qaol==
||u - Qucll+ + (Qn - Qn) (Qu QnC-=
(10)
|| Qxcl 2c (QR) [Q
The first term is GMsFEM error (which is small), and the second term is small if I (Qn) is small. The norm above can play an important role in the proof. Since the GMsFEM error is done in energy norm, in general, one needs to take the energy norm, which non-local and can slow down computations.
This error estimate gives bound between the learned solution and the angle between the spaces. To estimate the angle via the solution error is more difficult. Indeed, these questions need to be addressed in the future and will help to choose the appropriate loss functions.
D CSAL-PR
The primary role of LSAL 5 to enforce geometric alignment between the learned subspace It' and the target GMsFEM subspace R'. However, this loss has specific limitation: it is invariant to rotations within the subspace. If Qu. = QR J for any unitary matrix U (s.t. U. = I), the subspaces are considered identical under this metric, as:
No 11Q Q NM - |11| =0
This geometric alignment may overlook finer discrepancies in how specific functions are projected. The LSAl-PR (6) term was introduced to address this theoretical gap by directly enforcing projection consistency. It tests the subspaces with random vectors vf drawn from the target subspace R'. Since o' lies in R' its projection via the true basis is itself: PR. v² v². The discrepancy (I P v thus represents the projection error. Minimizing this forces the learned subspace to correctly capture arbitrary vectors from the target subspace, going beyond mere geometric overlap.
E THE SIGN INVARIANCE OF BASIS FUNCTIONS
Multiscale basis functions w are derived from local eigenvectors of' via 1 =x where Xi is partition of unity function and df solve the symmetric generalized eigenvalue problem:
A Φ
A. Φ
where A and S are symmetric matrices. These eigenvectors satisfy the orthogonality relations:
0 A du
Azdjer
0f S du 6jk-
Since eigenvectors are defined only up to a scalar multiple, if Φ; is an eigenvector corresponding to eigenvalue Aj. then -Ø; also valid eigenvector for the same eigenvalue. Both choices satisfy the orthogonality and normalization conditions above, meaning the sign of each basis function is arbitrary and does not affect its mathematical properties.
To address this sign ambiguity, define the RBFL2 loss as:
J
LRBFL = Es, min
where bf are the final multiscale basis functions (typically obtained by multiplying eigenvectors by partition of unity functions). This loss compares the prediction وif against both vi and - -b. selecting the smaller emrr account for the sign invariance.
15

---

Published as a conference paper at ICLR 2026
F INPUT DATA
We the Karhunen-Loewe expansion (KLE) Wong (1971); Aarnes & Efendiev (2008): Vasilyeva
et al. (222)) to generate stochastic permeability fields. This method decomposes a random field into deterministic spatial functions and random coefficients.
1. Covariance Function. We assume the covariance function has an exponential form:
R(x, v) oR exp(-vA2)
with
A
[1
D2
12
11
1/2 1
12
12
for 2D case and
A
21
22
12
1//1
1/2
12
|2]
22
1
12
for 3D case with correlation lengths I, Ly, 1 and variance ce
For the 2D case: I, = 0.02, ly = 0.6, o 22
For the 3D case: bg = 0.02, by = 0.6, 1 = 0.2, o =
2. Eigenvalue Problem. The eigenfunctions Φu and eigenvalues As are obtained by solving the homogencous Fredholm integral equation:
J R(x, (ddn))dd (x) k k= 1,2,.....
Random Field Construction. The random field is represented as:
(x,,w) Σ (w)x(
where (w) are scalar random variables, and L is chosen capture most of the field's energy by retaining the largest eigenvalues.
Permeability Field Generation. Each stochastic permeability field is defined as:
<(z, w) exp(ax ((x,w)):
where ((x, w) represents the heterogeneous porosity field derived from Yz (z( (w, and ax > O scaling parameter that controls the contrast.
This KLE framework provides systematic approach for generating realistic permeability fields with prescribed spatial correlation structures. An example of 2D input coefficient field x(z)) is shown in Fig.
The spatially variable forcing term is defined by
f(x) ~
where N denotes Gaussian random field. The parameters are set follows:
For the 2D case: 7 2000, a 1, and |3 0.5;
For the 3D case: 7 2000, a 2, and 8 1.
An example of 2D right-hand side fxx) shown in Fig.
16

---

Published as a conference paper at ICLR 2026
9600
8400
2,7
7200
1.0
4800
3600
-0.9
2400
-188
1200
-3.6
(a) Input coefficient, <(x) with values in range [1, 9600]
(b) Right-hand side, f(x)
Figure 4: Example of input coefficient and right-hand side.
G BASELINES VS. GMSFEM-NO
We compare GMsFEM-NO with several baselines:
1. POD Volkwein 2013). Classical global intrusive POD.
2. Intrusive POD with DeePONet/FFNO or POD basis Meuris et al. (2021). Meuris et al. (2023. First selected neural network is trained on standard regression problem. After that extract basis from trained network and uses similar to intrusive POD to form reduced model. For FFNO Tran et al. (2021) basis is extracted from the last hidden layer, for DeepONet Lu al. (2019) basis is extracted from trunk net.
3. PCA-Net Hesthaven & Ubbiali 2018), Bhattacharya et al. (2021). POD is used compress features and targets, MLP is used as processor.
4. Kernel Batlle al. (0044. Vector RKHS method is used map sampled input functions to sampled output functions.
5. DeepPOD Franco et al. projector-based loss.
A DL-based techniques used to directly learn optimal basis with
We use dataset with spatially variable forcing term (9) covered in more details in Appendix Neural networks was trained and evaluated on grid 100 x I00.
For each selected baseline we perform sweep over hyperparameters:
Intrusive POD with FENO basis. Architecture defined by the number of features in the hidden layer, number of modes used by spectral convolution, and number of layers. Number of features in the hidden layer was fixed to 64, number of modes was selected from the set (10, 14, 16). number of layers from the set [3, 4, 5). Optimisation was performed with Lion optimiser Chen et al. (2023) with exponential decay 0.5 with number of transition steps selected from [00,, 200), and learning rate selected from [5 10 5 10_4]. We optimise for 1000 epoch with batch size 10. In all architectures we used GELU activation function.
2. Intrusive POD with DeepONet basis. Architecture is defined by trunk and branch nets. As trunk net used convolution architecture with spatial downsampling by factor of 2 along each dimension after each layer, simultaneously, the number of channels was multiplied by 2 after each layer, branch net used standard MLP. We apply optimisation similar the one of FFNO, but select learning rate from [10-3 10 10-4J. Number trunk network layers was fixed
features, kernel size of convolution in trunk was selected among [3, 7]. In the branch net we vary number of layers [3, and the number of basis vectors [100, 200) in the last layer.
trunk encoder transformed 2 input features
either
3. DeepPOD Grid search for DeepPOD exactly the same for Intrusive POD with FFNO.
17

---

Published as a conference paper at ICLR 2026
4. PCANet. For PCANet the optimisation similar to Intrusive POD with DeepONet, but with 3000 epochs. We vary the sizes of POD encoder and decoder among [100, 500, and [100, 300, 500]. For MLP processor vary the number of layers [3. 4, 5] and the number of hidden neurons [100, 300, 500]-
5. Kernel. closely followed code provided by authors. As kernels we used Matern, RBF We combined the method with POD and performed grid search over the number of modes: [50, 100, 150, 200] for both features and targets.
Table 10: Regression-based methods.
method
train emror
test erTor
GMsFEM-NO
2.6%
2.8%
PCANet 6%
24%
kernel
7%
100%
Comparison of regression-bedd baaed approaches with GMsFEM-NO appears in Table TO] We observe significant overfitting for kernel-based method and PCA-Net.
Intrusive techniques are compared in Figure We see that bases extracted from DeepONet and FFNO generally not appealing. FFNO slightly improves over global POD (weak baseline) for = 64 basis functions. DeepONet fails to reach accuracy of global POD. The most competative approach is DeepPOD. Note however, that DeepPOD becomes comparable to GMsFEM-NO with 8 sparse (localised) basis functions only when it uses 20 dense basis functions. With 30 basis functions DeepPOD outperforms GMsFEM-NO with 8 basis functions.
81
POD
DeepONet
FFNO
DeepPOD
relative error,
23
GMsFEM(8)
6.5
2.8
2.2
10
50
100
150
200
Naasis
Figure 5: Comparison of accuracy for intrusive techniques. Number of basis functions for GMsFEM fexed to 8. DeepONet and FFNO mean Intrusive POD with DeepONet/FFNO basis.
H GMsFEM-NO TRAINING DETAILS
For F-FNO training to predict basis functions subspace, we used AdamW optimizer Loshchilov Hutter (2017) with cosine decay learning rate scheduler. The initial learning rate was 10 We trained NO 600 epochs. The best number of random test vectors n² 10 (see the Table For the rest of the hyperparameters,
performed grid search:
1. Batch size [8, 32].
2. Number of operator layers [4, 5].
3. Number of modes used in F-FNO kemnel:
For 2D:
18

---

Published as a conference paper at ICLR 2026
for full domains [16, [18,18]] 18]
for half domains [[8,8], [10, 10), [14, 8), [14, 10]]
for corner domains [[6, 6), [8, [10, 10]].
For 3D:
for full domains [[6,6, 6,. [8, 8,8]
for half domains [[8,8,4] [6, 6,3]
for quarter domains [[8, 4, 4), [6,333]
for corner domains [[4, 4, 4], [3,
4. Number of channels in the FENO kernel [64, 128].
Table 11: Performance GMsFEM-NO for 2D (100 x 100, N. 36).
GMsFEM-NO, LSAL
v²
5
v¹
10
v²
15
20
Dataset
L2
H1
L2
H1
L2
H1
L
H1
Diffusion,
1.09%
11.90%
1.06%
11.65%
1.31%
12.34%
1.10%
12.13%
Diffusion,
2.84%
19.14%
2.81%
19.03%
2.88%
19.38%
2.87%
19.37%
Richards, 1.91%
11.21%
1.87%
11.25%
1.98%
11.75%
1.94%
11.64%
Richards, 9
2.91%
20.06%
2.99%
19.60%
2.96%
20.33%
2.94%
20.26%
The source code containing the optimal parameters will made publicly available upon acceptance. We use JAX, Optax DeepMind et al. (2020) and Equinox Kidger & Garcia (021) in all experiments.
I F-FNO, GNOT, TRANSOLVER++ TRAINING DETAILS
I.1 F-FNO TRAINING DETAILS
For F-FNO Tran et al. (2021) we used the following training protocol. We employed the AdamW optimizer Loshchilov Hutter (2017) with a cosine decay learning rate scheduler and trained for 600 epochs, using base learning rate of 10-3, We performed grid search over the following hyperparameters:
1. Batch size: [8, 16,32] 32];
2. Number of modes in F-FNO kernel: [14, 16];
3. Number of operator layers: (4,5];
4. Number of channels in F-FNO kernel: [64, 128).
The optimal hyperparameters were: batch size 8, 5 operator layers, 16 modes, and 128 channels.
1.2 GNOT TRAINING DETAILS
For GNOT Hao et al. (2023), we employed the following training protocol: AdamW optimizer with onecycle learning rate scheduler for 600 epochs, using base learning rate 10-3 We conducted grid search over hyperparameters that were previously found optimal for different 2D problems in
Hao et al. (023)::
1. Batch size: [4, 8]
2. Number of attention layers: [3, 4]
3. Hidden size of attention and input embeddings: [96, 128, 192]
4. Number of MLP layers: [3, 4]
19

---

Published as a conference paper at ICLR 2026
5. Hidden size of MLP: [128, 192]
6. Number of heads: [4, 8]
7. Number of experts: [3, 4
The optimal configuration was: batch size 4, 4 attention layers, hidden size of 128 for attention, MLP, and input embeddings, 8 heads, 4 experts, and 4 MLP layers.
1.3 TRANSOLVER++ TRAINING DETAILS
For Transolver++ Luo et al. 2025 we employed the following training protocol: Adam optimizer with onecycle learning rate scheduler for 600 epochs, using base learning rate of 10-3. We conducted grid search over hyperparameters previously identified as optimal for 2D problems in Luo al. 2025):
1. Batch size: [4,8] 8]
2. Number of layers: [4, 8]
3. Hidden size: [128, 256]
4. Number of MLP layers: [1,2]
5. Number of heads: [4, 8]
6. Slices: [32, 64]
The optimal configuration was: batch size 4, 8 layers, hidden size 256, 2 MLP layers, heads, and 64 slices.
J RESULTS OF GMsFEM-NO FOR NBF 4
Table 12: Performance comparison loss functions for NO training (100 x 100 grid, N, 36).
LRBFL
CSAL
ESAL-PR
Nbr
Dataset
L2
H1
L2
H1
L2
H
Diffusion,
2.72%
19.10%
2.39%
18.01%
2.40%
18.01%
4 Diffusion,
6.03%
29.10%
5.76%
28.40%
5.78%
28.43%
Richards, 3.87%
16.26%
3.14%
15.01%
3.17%
15.02%
Richards,
9.78%
34.39%
6.11%
29.37%
6.13%
29.42%
For the Richards equation with complex right-hand side 99 and Nor = 4 basis functions, our proposed loss improves the relative L2 metric by factor of 1.6.
We did not conduct further experiments with GMsFEM-NO using Nor = 4 because its performance was insufficient and underperformed compared the standalone neural operator.
K OUT-OF-DISTRIBUTION RESULTS
Unlike standalone NOs, which suffer from catastrophic failure when applied to out-of out-oo-diistriution forcing terms (Table 13,, GMsFEM-NO is fundamentally independent of the right-hand side, ensuring robust performance. Retraining the NO for new right-hand side terms requires computationally expensive recalculation of solutions, highlighting key limitation of standalone NO learning.
20

---

Published as a conference paper at ICLR 2026
Table 13: Out-of-distribution results for the NO: training and testing on PDEs with different right-hand sides.
Train, Drain
Test, Prest
100 x 100
250 250
Diffusion, Diffusion, Diffusion.
Diffusion,
218%
174%
1392%
1632%
Richards, Richards, Richards, Richards,
196%
113%
6503%
6554%
L HEAT EQUATION FOR 2D
We consider the heat equation with heterogeneous coefficient
du
Öt V (k())Vu(x)) = f(x), IE x (0, Tmax];
= 0, E N x (0, Tmma
h= 0, E Ωt
where time parameters: At = 2.5 10 Tmax = 10 10 We consider two right-hand side configurations:
Uniform unit forcing term
f(z)
(11)
Spatially variable forcing defined by
[(x) sin (TTr)) cos(TTy)
(12)
Table 14: Performance comparison of GMsFEM and GMsFEM-NO for Heat equation for 2D (250 250, N₂ = 121).
x
GMsFEM
GMsFEM-NO, LSAL
Nbr
Dataset
L2
H1
L2
H1
Heat,
0.72%
10.81%
0.78%
13.06%
8
Heat, 12
1.13%
19.92%
1.17%
21.57%
The results, presented in the Table [4] show that our proposed GMsFEM-NO method with LSAl maintains high accuracy for this time-dependent problem, with relative L2 errors below 1.2%. This demonstrates the successful application of our framework to time-dependent problems.
M DIFFUSION EQUATION WITH MIXED BOUNDARY CONDITIONS
We consider the diffusion equation with heterogeneous coefficient
V (x(x)Vu(x)) f(x), .EΩ= (0, 1)
0, I E TED.
Du a (u - Mconst) x EETR.
k(x) On
where the boundary is partitioned into Dirichlet part Tp fI 02 y = 0} and Robin part I'R = {i 0Q y = 1} with parameters: Robin coefficient a 1, Mconst = 1.
Our framework demonstrates robust performance when extended to problems with mixed boundary conditions, maintaining high accuracy even in these more complex scenarios. As shown in Table [5] the GMsFEM-NO method with LSAL achieves relative L2 errors of approximately 1.2% for diffusion problems with Robin boundary conditions. This represents a significant extension beyond the homogeneous Dirichlet conditions typically considered in multiscale methods.
21

---

Published as a conference paper at ICLR 2026
Table 15: Performance comparison of GMsFEM and GMsFEM-NO for equation with mixed boundary conditions for 2D (250 x 250, N. 121).
GMsFEM
GMsFEM-NO, LSAL
Nbr
Dataset
L
H1
L2
H1
Diffusion,
1.11%
7.98%
1.22%
12.10%
Diffusion,
1.09%
7.72%
1.23%
12.31%
N TIME
The core objective of GMsFEM-NO is reduce the cost of the traditional GMsFEM offline stage over many subsequent simulations. To provide full transparency, we have performed detailed analysis that compares two practical scenarios: training models sequentially (one after another on one GPU) and in parallel. The breakeven point x is the number of inference samples where the total time of GMsFEM-NO (including training data generation and training) equals that of using standard GMsFEM for all samples:
datta + Tusin Tiat T TGMsFEM
where Tama is the wall-clock time for generating training data via traditional GMsFEM, Train the total wall-clock time for training NO (which varies by scenario), TeM.FEM the wall-clock time for single offline GMsFEM basis generation, and Ti the wall-clock time for single GMsFEM-NO inference. The results are in the Table [6
Even with sequential training, the method pays off for multi-query scenarios. The benefit is larger for large-scale 3D problems, where the breakeven point remains low (159 samples) due to the high cost of the traditional GMsFEM solver.
Table 16: Performance comparison of sequential and parallel implementations
Problem
Performance Metries
Configuration
TOMPFE
Tiat
Samples
Implementation
T Train
100 x 100
16.87
0.28
800
Sequential
37800
-3270
16.87
0.28
800
Parallel
12600
-1570
250 x 250
210.5
0.31
800
Sequential
43200
1040
210.5
0.31
800
Parallel
14400
870
50 x 50 x 50
935.4
0.84
800
Sequential
86400
1090
935.4
0.84
800
Parallel
21600
820
100 X 100 X 100
10547.2 2
1.33
150
Sequential
100800
-159
10547.2
1.33
150
Parallel 25200
-152
22


Reviews:
[{"review_id": "3SVBrL48Tk", "reviewer": "Reviewer_zYEo", "summary": "The paper tackles multiscale, high-contrast PDEs by replacing the expensive local eigenvalue problems of GMsFEM with a neural operator that predicts the local coarse spaces directly. Instead of regressing individual basis functions, the method learns the subspace they span and trains with a Grassmann-geometry–aware subspace alignment loss (with an optional projection regularizer). The learned spaces assemble into a restriction operator that preserves the classical GMsFEM solve, yielding GMsFEM-level accuracy while reducing offline basis construction by more than an order of magnitude. Experiments on diffusion and steady Richards equations in 2D and 3D show strong accuracy, data efficiency, right-hand-side robustness, and resolution transfer.", "questions": "Could you provide a bound that maps principal angles between true and predicted coarse spaces to a bound on the $𝐻^1$ relative error of the final solution, at least under simplifying assumptions?\nFor 3D large-scale runs, what are the peak memory, wall-clock breakdown, and parallelization strategy during inference and coarse solves?", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The conceptual shift from matching basis functions to matching their span is clear and well motivated, aligning the learning objective with what GMsFEM truly needs. The subspace alignment loss directly optimizes principal angles on the Grassmann manifold, removing sign and ordering ambiguities and leading to stable training. The framework preserves the classical coarse-solve pipeline, so accuracy remains comparable to GMsFEM while offline cost is drastically reduced; the empirical speedups are practically meaningful. Results demonstrate robustness to right-hand-side changes and cross-resolution transfer, where direct solution-predicting NNs typically struggle. The paper is clearly structured, with a clean train→assemble→solve story that is easy to reuse in other PDEs.", "weaknesses": "Theoretical guarantees connect subspace error to final solution error only implicitly; a formal upper bound from principal angles to \n$𝐻^1$error would strengthen soundness. Evaluation focuses on structured grids with Dirichlet boundaries and steady problems; non-structured meshes, mixed or Robin boundaries, and time-dependent or strongly nonlinear systems remain open. Comparisons omit graph- or mesh-based neural operators and learning-augmented multigrid with learned prolongation/restriction, which are natural baselines here. The model zoo created by per-region-type networks increases operational burden at scale; an ablation on parameter sharing or conditional modulation would be helpful. Wall-clock reporting includes basis generation speedup but a comprehensive end-to-end budget (training time, energy, peak memory, coarse-solve cost) would improve practical transparency.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "tc5TDohzEn", "reviewer": "Reviewer_9cmw", "summary": "The paper proposes **GMsFEM-NO**, a hybrid framework that uses a neural operator to accelerate the **offline stage** of the Generalized Multiscale Finite Element Method (GMsFEM). Instead of solving many local eigenvalue problems to construct multiscale basis functions, the method trains an NO to map a heterogeneous coefficient field (\\kappa(x)) directly to the **subspace** spanned by the GMsFEM bases. The key idea is a **Subspace Alignment Loss** ((\\mathcal{L}_{\\text{SAL}})) that aligns predicted and true multiscale subspaces on a Grassmannian, avoiding the sign/permutation ambiguity of individual eigenvectors. The approach preserves two important GMsFEM properties: (i) independence from the forcing term (f(x)), so it generalizes to OOD right-hand sides, and (ii) **resolution invariance** (train on coarse, use on fine). Experiments on 2D/3D high-contrast diffusion and steady-state Richards equations show **>60× speedup** in basis generation while matching GMsFEM-level accuracy.", "questions": "1. Several tables show GMsFEM-NO outperforming the classical GMsFEM it is meant to approximate. Can you provide a more concrete explanation than “statistical variation” (e.g. averaging over multiple local domains, extra smoothing from the NO, or differences in quadrature/projection)?\n\n2. For $(\\mathcal{L}_{\\text{SAL-PR}})$, how sensitive is the reported improvement to the choice and number of random test vectors \n$v^{i}$? \n\nWould the simpler $\\mathcal{L}_{\\text{SAL}}$ suffice in most practical cases?\n\n3. Can you add comparisons to more recent neural operators (beyond F-FNO / POD-style baselines) to better justify that GMsFEM-NO is competitive not only with classical GMsFEM but also with current NO-based multiscale solvers?", "limitations": "", "rating": 4, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "1. **Targets the real bottleneck.** It identifies the most expensive part of GMsFEM (local spectral problems in the offline stage) and replaces it with an NO, which is a practical and impactful acceleration.\n2. **Subspace-level supervision.** The proposed $(\\mathcal{L}_{\\text{SAL}})$ is technically meaningful: learning the *space* of local bases is more robust than learning each basis function, and it naturally resolves sign ambiguity; this is the main methodological contribution.\n3. **Strong empirical claim: fast but still GMsFEM.** The method demonstrates >60× speedup in basis construction while retaining essentially the same $(L_2/H^1)$ errors as classical GMsFEM on multiscale benchmarks, which makes it attractive for large runs or parameter sweeps. \n4. **OOD and resolution generalization.** Because the solver structure of GMsFEM is kept, the approach remains independent of $(f(x))$ and supports train–coarse / test–fine usage, which standard NOs typically fail to do.", "weaknesses": "1. **Marginal gain from $(\\mathcal{L}_{\\text{SAL-PR}})$.** The paper introduces a more complex variant with projection regularization, but the improvement over plain $(\\mathcal{L}_{\\text{SAL}})$ is small (e.g. 1.82% → 1.72% on 250×250), and only on a single setting; the extra term is not clearly justified. \n2. **Systematically “better than” the target.** In several tables, GMsFEM-NO slightly outperforms the original GMsFEM it is approximating; calling this “statistical variation” is not entirely convincing because it appears repeatedly, suggesting either measurement noise, data leakage, or a smoother inductive bias from the NO. This needs a clearer explanation. \n3. **Baselines not fully up to date.** Most comparisons are to F-FNO and classical reduced-order/ POD-style methods; stronger modern neural operators for multiscale/high-contrast settings are missing, so it’s hard to tell where GMsFEM-NO sits against the current SOTA.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "tRgAi8ieBb", "reviewer": "Reviewer_EehQ", "summary": "The paper extends the GMSFEM by incorporating a neural operator to accelerate the process. The work proposes a subspace alignment loss function to learn coherent subspaces, and introduces projection regularization terms to enforce consistent projections. Experiments compare the proposed loss function with the conventional loss function, GMSFEM-NO with GMSFEM and F-FFO, for evaluations.", "questions": "Are there failed cases that don't converge to the solution?\n\nIs it able to extend the method to time-dependent equations?\n\nIs it possible to further reduce the error if more basis functions are used? What's the limits of the error?", "limitations": "", "rating": 4, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The combination of GMSFEM with NO is an interesting idea.\n\nThe proposed SAL and PR losses can train the NO effectively.", "weaknesses": "The method is only tested with m moderate-scale data in rectangular domains. Is it able to solve larger-scale problems in non-rectangle domains?\n\nThe generalization ability of the learned NO is not evaluated sufficiently. The usability of the method is unclear.\n\nOnly zero Dirichlet boundary condition equations are tested? Can it solve Dirichlet boundary condition equations with other values?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "A7lhSbP3fa", "reviewer": "Reviewer_E9D5", "summary": "The manuscript introduces the GMsFEM-NO framework, which integrates the Generalized Multiscale Finite Element Method with NOs to improve the efficiency of existing NO frameworks when solving multiscale and high-contrast PDEs. Legacy GMsFEM requires eigen solvers, which are computationally expensive. GMsFEM-NO legacy eigen solvers with a NO trained to efficiently predict the local multiscale basis subspaces. A Subspace Alignment Loss, with a regularized variant, is also proposed to enhance global reception and ensure geometric and physical consistency. Experiments on 2D and 3D benchmarks show significant speedup over legacy GMsFEM with a comparable accuracy. Other technical merits include better generalization to OOD enforcing terms (v.s. NOs) and resolution invariance.", "questions": "- Could the same SAL losses be extended to time-dependent GMsFEM formulations or to non-elliptic PDEs?\n- Can GMsFEM-NO handle non-rectangular or unstructured meshes if implemented with GNOs or other NO variants, e.g., GNOT, Transolver, etc.?", "limitations": "", "rating": 6, "confidence": 5, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- **Mathematical grounds**: Links subspace alignment to Grassmann manifold distance, enhancing theoretical credibility. The new SAL and SAL-PR losses are well-motivated and mathematically tied to Grassmannian geometry. \n- **Comprehensive experiments**: Includes linear/nonlinear PDEs, 2D/3D cases, various forcing terms, and extensive comparisons vs. multiple baselines.\n- **Performance boost**: 60x speedup v.s. legacy GMsFEM with comparable accuracy; better OOD generalization vs. standalone NOs; resolution invariance.", "weaknesses": "- **Limited scope**: The experiments are limited to elliptic problems. The study does not include time-dependent PDEs, limiting the demonstrated generality. FEMs and NOs are known to be robust to hard problems, e.g., Navier-Stokes. This limits the scope of the proposed framework. \n- **Assumption of regular grids**: Current implementation works only for regular grids due to the use of F-FNO; performance on irregular geometries or adaptive meshes remains unexplored.\n- **Dependence on GMsFEM setup**: The framework still relies on a coarse grid and precomputed spectral problems for training data. It does not remove the expensive stage, but only accelerates it.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)