Neural Green’s Functions
Authors: Seungwoo Yoo
Kyeongmin Yeo
Jisung Hwang
Minhyuk Sung
Venue: NeurIPS 2025 poster ()
Keywords: Neural Operators
Primary area: machine_learning_for_sciences

Abstract: We introduce Neural Green’s Function, a neural solution operator for linear partial differential equations (PDEs) whose differential operators admit eigendecompositions. Inspired by Green’s functions, the solution operators of linear PDEs that depend exclusively on the domain geometry, we design Neural Green’s Function to imitate their behavior, achieving superior generalization across diverse irregular geometries and source and boundary functions. Specifically, Neural Green’s Function extracts per-point features from a volumetric point cloud representing the problem domain and uses them to predict a decomposition of the solution operator, which is subsequently applied to evaluate solutions via numerical integration. Unlike recent learning-based solution operators, which often struggle to generalize to unseen source or boundary functions, our framework is, by design, agnostic to the specific functions used during training, enabling robust and efficient generalization. In the steady-state thermal analysis of mechanical part geometries from the MCB dataset, Neural Green’s Function outperforms state-of-the-art neural operators, achieving an average error reduction of 13.9% across five shape categories, while
being up to 350 times faster than a numerical solver that requires computationally expensive meshing.

Neural Green's Functions
Seungwoo You
Kyeongmin Yeo
Jisung Hwang
Minhyuk Sung
KAIST
{dreamy1534, aaaaa, 4011hjs, wwsung ac.kr
GT (Sduiioo)
Trunshee [Wa
Wangl
Figure 1: Neural Green's Function. We propose neural solution operator designed for linear PDEs whose differential operators admit eigendecompositions. Our framework effectively handles irregular geometries and diverse source and boundary functions. In the steady-state thermal analysis used model problem, method demonstrates superior generalization performance compared to state-of-the-an neural operators in predicting ground-truth solutions (column 2), as evidenced by the error maps (columns 3-6).
Abstract
We introduce Neural Green's Function, neural solution operator for linear partial differential equations (PDEs) whose differential operators admit eigendecompo- sitions. Inspired by Green's functions, the solution operators of linear PDEs that depend exclusively imitate their behavior, achieving superior generalization across diverse irregular ge- ometries and source and boundary functions. Specifically, Neural Green's 's Function extracts per-point features from volumetric point cloud representing the problem domain and uses them to predict decomposition of the solution operator, which is subsequently applied to evaluate solutions via numerical integration. Unlike recent learning- baaeing-basedd solution operators, which often struggle generalize unseen source or boundary functions, our framework is, by design, agnostic to the specific functions used during training, enabling robust and efficient generalization. In the steady-state thermal analysis mechanical part geometries from the MCB dataset, Neural Green's Function outperforms state-of-the-arr neural operators, achieving an average error reduction of 13.9% across five shape categories, while being up 350 times faster than numerical solver that requires computationally expensive meshing.
the domain geometry, we design Neural Green's 's Function to
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

Introduction
Linear partial differential equations (PDEs) play a central role in many scientific and engineering disciplines, such as thermal analysis, electrostatics, fluid dynamics, and elasticity. Similar to other PDEs, their ubiquity has driven the development of numerical methods, such as finite difference methods (FDMs) and finite element methods (FEMs), to solve the equation defined over complex geometric domains in real-world applications where analytical solutions are often unavailable. These techniques discretize problem domains into linear systems that can be constructed and solved efficiently. However, the reliance on meshes poses significant limitation, as generating them requires running computationally expensive meshing algorithms [13 [14), making it challenging to rapidly evaluate solutions across multiple problem domains with different boundaries. This bottleneck is especially pronounced during the early design phase of engineering workflows, where rapid iteration is essential for achieving optimal results.
meshes of lattices polygons, resulting in
Recently, learning-basse -bsedd solvers have emerged as promising surrogates to conventional numerical solvers for solving PDEs. They predict solutions without mesh construction, either by optimizing neural network parameterize the solution function for given problem instance [30, 19 or by learning solution operators that map input functions to solutions, enabling inference with a single
While these approaches significantly improve efficiency compared to numerical solvers and provide less expensive approximations, little attention has been paid to their generalization capabilities under simultaneous variations problem domains, source functions, and boundary functions. Some previous methods focus on learning the solutions to PDE problems within a single domain [30 25 27,, requiring re-training whenever the problem domain or even a specific problem instance changes. Approaches such as Boullé et al. (6], Teng al. (31], and Negi et al. [28] address this by learning the solution operator for a given domain, enabling the handling of different problem instances. However, they still require re-training when applied to different domains. Others 20 26 23
forward pass 20 22 26 23
36
34
3)
36 35 attempt to learn solution operators that generalize across collection of shapes but often struggle to handle unseen source functions and boundary functions. This limitation arises because these inputs are directly provided to the neural networks that regress the solution operators, requiring sufficient training examples to achieve generalization.
In this work, we propose a learnable solution operator capable of handling various shapes and func- tions, with a specific focus on linear PDEs whose differential operators admit eigendecompositions. such as the Poisson and Biharmonic equations. The foundation of our approach is rooted in the mathematical definition of the solution operator for a linear PDE, known as the Green's function, which only depends on the geometry of the problem domain. Once the Green's 's function of domain known, the linearity allows solutions to be computed for arbitrary source and boundary functions. This unique property motivated the design of our framework, which predicts neural features solely from the domain geometry. These features are used to approximate the eigendecomposition of the corresponding Green's 's function and to predict related differential quantities, such as per-vertex masses, that are essential for computing solutions via numerical integration. This design imposes a strong prior on the solution operator learned by the neural network, making the model independent of the specific source and boundary functions used during training.
In our experiments, we first empirically validate the generalizability of the proposed framework using simple examples of the Poisson and Biharmonic equations. We then extend the evaluation to a practical setting of steady-state thermal analysis on complex 3D geometries, using a benchmark we built from the MCB dataset [16]. This benchmark comprises a diverse collection of mechanical part geometries, each paired with a steady-state thermal distribution obtained by solving Poisson's 's equation. The dataset features a wide variety of shapes with significant intra-class variations, along with triplets of source, boundary, and solution functions generated by numerical solver, providing challenging benchmark setup for both baselines and our method. In comparisons, Neural Green's Function demonstrates superior generalizability compared to sfate-of-the-ar neural operators [35] 34 3).
Notably, our framework reduces the error metric by 13.9 % compared to Transolver [35]. while sharing the same backbone. These results highlight the significance of explicitly formulating and predicting solution operators in achieving superior generalization capabilities.
2 Related Work
Physics-Ifformed Neural Networks (PINNs). Physics-Informed Neural Networks (PINNs) 30
parameterize the solution of each PDE instance directly using a neural network, and optimize its
2

---

parameters by minimizing objective functions derived from the governing equations and empirical observations. They can be easily implemented using automatic differentiation, which supported by deep learning frameworks 22 2 71. However, PINNs are trained separately for each problem instance, with each neural network approximating single solution at time. This approach requires retraining whenever the problem domain, source function, or boundary condition changes.
Neural Operators. Unlike PINNs, Neural Operators learn function-to-function mappings that act as solution operators for PDEs. They achieve this by parameterizing (nonlinear) kernel functions and implicitly perform kernel integration across subsequent network layers. Specifically, GNO [20] employs graph neural networks (GNNs) to approximate integral kernels through local aggregation via cassage-passing FNO 21 improves efficiency and performance by utilizing fast Fourier transform (FFT) [00 in regular domains and learning global integral kernels in the spectral domain. Follow-up work [22 26] combine GNO [20] and FNO [2T] to handle irregular problem domains by efficiently
modeling both local and global interactions within systems. Several work [24] 23 222 33 utilizing
Transformers 33 architecture have also been proposed, leveraging self-attention mechanism to
capture interactions among mesh points. These work bypass the quadratic complexity of self- self-attention layers using linear Transformers [8] 8 Recent work further reduce the computational complexity by incorporating compact latent representations with dimensionalities significantly smaller than mesh sizes [35] 34 3 Neural operators are a versatile framework for operator learning over irregular geometries with diverse source and boundary functions. However, they often couple input meshes with sampled function values, which makes generalization to new functions challenging.
Learning Green's Functions. Several recent work focus linear PDEs and propose data-driven approaches to discovering the Green's function of a problem domain, enabling solutions for varying source and boundary functions. Boullé et al. (6] propose a method to learn Green's 's functions for linear PDEs using pairs of forcing and solution functions, leveraging rational neural networks [5]. Teng et al. K extend this approach by enabling the evaluation of solutions with arbitrary source and boundary functions, regressing Green's 's functions through the approximation of Dirac delta functions with Gaussian distributions. Similarly, Negi et al. [28] approximate free-space Green's functions using radial basis function (RBF) kernel-based neural network. By directly learning Green's 's functions, the aforementioned approaches can predict solutions for PDEs with varying source and boundary functions. However, these methods are restricted to individual problem domains and require retraining for unseen geometries, as the networks learn domain-specific Green's functions and their gradients. Additionally, they focus on simple domains (e.g.. circular domains in 2D), where numerical quadrature for computing solutions with predicted Green's 's functions relatively straightforward. Our work addresses these limitations through a framework design inspired by the eigendecomposition of Green's 's functions for linear PDEs. In particular, the framework learns to extract geometric features from problem domains and composes them reconstruct the corresponding Green's 's functions. Moreover, our approach predicts differential quantities needed to approximate numerical integration over complex. irregular geometries, extending its applicability to real-world scenarios, as demonstrated in Sec.
3 Background
Consider a linear PDE subject to Dirichlet boundary conditions defined over a continuous domain D C I with boundary aD:
[Cu(x)=f(a) fx) (u(x) = h(x))
x x
(1)
x xEAD,
where L is linear differential operator, and u, f. and are functions that reside in Banach spaces U4, F. and H. of functions, respectively. The solution operator & for Eqn. I that computes the solution
u = 9 (D, f. h) is the mapping:
G: x F x H U,
where D is set of problem domains. Since Eqn. iis linear PDE, & is explicitly given as: u (x) GD (f,h)
/ GD
(3)
+
h(y) v GD (z. v) nly)dy
BD
3

---

where GD == .D)...) is an instantiation of G in the domain D and n2y) is the outward normal vector at y E OD, respectively. The integral kernel Gn D x D - R, referred as the Green's function of the operator L in D, represents its impulse response and is a solution of:
(G)( v) = 6(z v), I DED D Gz(zzy)=
(4)
E aD
where d is the Dirac delta function. One important property of 9D and Green's function GD is that they only depend on the geometry of the domain D, regardless of the specific source function or boundary function h. This implies that, once the Green's 's function Gp for a given domain D known, the PDE in Eqn. T can be solved for arbitrary source and boundary functions ʃ and h. However, solutions based on Eqn. havee only been applied to simple domains where closed-form expressions for Green's 's functions are available.
Numerical solvers, on the other hand, utilize volumetric mesh representing D to discretize the operator C in Eqn. forming a discrete counterpart represented as linear systems. Assume that D is represented as mesh D = (V, T) of a single connected component, consisting of vertices V and faces T. We define selection matrix S E [0,1]
to indicate N, boundary vertices among the No vertices, and its complement K E {0. (N (N) to represent the remaining interior vertices. The scalar-valued source function and the boundary function h are discretized by sampling their values at the vertices, resulting in two vectors: E N. and h E NN, respectively. Then, Eqn. can be written linear system in its discrete form:
Lu Mf
s.t Su=h. h,
(5)
where L E IRN.XN, X. is the discretization of L and M E NoXN. is the mass matrix of the mesh D, respectively. The solution to this system can be written as:
u == K (KLK))) (KMF - KLST h) + S
(6)
= K" (G(MM KLSTh)) + S
where G (KLK?) R.N.-NN) (No. represents the inverse of the matrix L restricted to
the interior of the domain. The detailed derivation can be found the Appendix. The matrix G
serves the same role as the Green's function discussed in Sec. 3 Notably, the terms G, M, and K
in Eqn. 6 are independent of the given functions and h, analogous to the Green's ' function GD
in Eqn. This independence implies that the solution operator for Eqn. [5 can be approximated by
predicting these terms using only the geometric information of the mesh D.
Assuming that L is symmetric, its submatrix KLKT,, obtained by eliminating the rows and columns corresponding to the boundary vertices, also symmetric. Furthermore, the imposed boundary conditions are sufficient, KLK becomes full-rank and admits the following eigendecomposition:
G --
where & E (MN.- is the matrix whose columns the eigenvectors of KLKT, and A E (N.. is the diagonal matrix of positive eigenvalues corresponding to the eigenvectors in d. These properties form the basis of our novel framework, designed to predict solutions of linear PDEs across diverse irregular geometries, detailed in the following section.
4 Neural Green's Function
Numerical solvers are well-established tools for solving linear PDEs; however, they rely the volumetric mesh D generated from surface representations (e.g., boundary representations in CAD), which can be computationally expensive. This poses challenge rapid iteration, which is essential in the early stages of design exploration.
may help mitigate this issue by enabling coordinate-based solution queries, eliminating the need for the mesh D. These models effectively capture complex behaviors of nonlinear systems by encoding both query points, along with source and boundary functions, into learned latent representations. However, for the class of linear PDEs discussed in Sec. 3 we claim that stronger prior can be incorporated into the
Directly applying neural surrogate models 22 26 24 23
36 351

---

Go
D
M
f
Lg
Figure 2: Method overview. Given query points representing the geometry of the problem domain where the solution function to be predicted, Neural Green's 's Function extracts neural features that are subsequently used to approximate the Green's 's function and to predict the differential quantities required to evaluate Eqn. The example rendered using a tetrahedral mesh, which is used solely for visualization purposes.
operator design to improve its generalization capability. In the following, we discuss how the design of Neural Green's s Function, inspired by Green's 's functions, enables generalization to unseen source and boundary functions. Although we mainly focus on problems 3D space, our approach is not limited to particular dimension.
Our framework, as outlined in Fig. 2 takes set of query points Q in 3D space that represents the geometry of the domain D. These points can be obtained through interior point queries 4 thereby avoiding the costly domain meshing. From these query points, we extract features that are used to approximate Green's 's functions at the corresponding points, rather than directly predicting the solutions. While the query points are associated with source function f evaluated at Q and boundary function h sampled at the boundary points, our network remains agnostic to these functions. This independence the key for achieving generalization unseen functions.
Specifically, we employ neural network with parameters @ predict features b R Q| for all query points, whose coordinates are provided as inputs to the network. These predicted features are then utilized construct the ground-truth Green's 's function G in Eqn. 6 Inspired by Eqn. [ we compute each element of Ga the dot product of feature vectors predicted for interior points. The computation of Go for all pairs of points is efficiently parallelized through a single matrix multiplication:
neural Green's 's function Ge, which serves drop-in replacement for
G = (K (K
(8)
Since the mass matrix M and the submatrix KLST of the operator L are unavailable without a mesh D, we also predict these quantities to evaluate Eqn. 6 In particular, we employ MLP to decode per-vertex features p into per-vertex mass values Mg, using the Softplus activation to ensure that the predicted masses are positive. Similarly to Gg. the submatrix L = KLS" predicted by computing dot products between features V E IR Q xd which are decoded from & by another MLP:
Ly = (o) (SMa))
(9)
Lastly, these components are combined to evaluate Eqn. 6
u K (Go [KMff Loh) S
(10)
yielding the predicted solution ug.
Given dataset of N examples {u', f' (('',f,h'))) h²) generated by solving Eqn. 5 using an FEM solver meshes (D²) Jim along with their corresponding
matrices {M+) ] i=11 the network parameters
0 are optimized end-to-end by minimizing the empirical risk:
arg min E,
[]] [| u 112 + k >| diag (M') [1]]
(11)
where À > 0 is regularization weight applied to the predicted mass values ensure consistency with the ground truth. We observe that regularizing the predicted masses is essential for the convergence of network training, discussed further in Sec. 5.4 In all our experiments, set À 1.
5

---

5 Experiment
5.1 Experiment Setup
In this section, we outline the baselines, evaluation metric, and implementation details of our experiments, while deferring detailed descriptions of the problem setups their respective sections.
Baselines. For 2D Poisson and Biharmonic examples in Sec. 5.2 we compare Neural Green's 's Function against Transolver [35 the state-of-the-art neural operator. In Sec. 5.3 we further expand the set of baselines to conduct a more comprehensive evaluation in steady-state thermal analysis. This includes Transolver (35], Latent Neural Operator (L.O)) (44, and Universal Physics Transformer (UPT) (3) For all baselines, we use their official implementations and train the networks under the same configuration as ours. For neural operators that require conditioning on source and boundary functions, we follow the approach of Transolver [35] and LNO [3],, concatenating the source and boundary function values at the query points as inputs to the baseline models.
the error e(u, uz) between the predicted solution ug and the corresponding ground truth u is measured as the relative L2 distance:
Evaluation Metrics. Following the baselines 35
34 3
| u//2
e(u, un) =
(12)
|uu|a
which is averaged over the entire test set for quantitative comparisons.
Implementation Details. We adopt the network architecture of Transolver 335 a state-of-the-art neural operator and one of our baseline methods, to extract the features $g. encouraging feature sharing across query points. Specifically, our network consists of an MLP that maps coordinates of query points to latent representations. The layer is followed by eight Transolver blocks (35] and MLP heads that decode them into the components required to evaluate Eqn. [O Unless otherwise specified, all models in our experiments are trained for 40 epochs using the ADAM optimizer [17 with a OneCycleL.R learning rate scheduler, setting the maximum learning rate to 1 x 10-1, For the steady-state thermal analysis experiments in Sec. [533 we use batch size of and accumulate gradients over 8 training steps to stabilize training, to handle meshes with large number of vertices.
5.2 Two-Dimensional Poisson and Biharmonic Equations
We first assess the generalization capability Neural Green's 's Function to unseen source and boundary functions, by fixing the problem domain to unit square [0, 1] x [0, 1] in 2D. The domain is discretized with the resolution of 100 x 100 by uniformly placing grid points along the two axes. As model problems, we consider Poisson and Biharmonic equations defined over this domain.
Problem instances of Poisson's 's equation are generated using boundary functions instantiated from the template that satisfies Au 0:
u(x, v) A(x² 3xy²) + Bly 33x + x²,
(13)
with the source function set to zero within the domain interior. We use 100 training and 100 test examples, each generated by sampling tuples of coefficients (A, B). To ensure that the boundary functions at test time do not overlap with those used for training, we sample coefficients from the uniform distribution U[[ 1] for the training set, and from 21[1, 2] for the test set.
For Biharmonic equation, we consider the template
u(x.y)= A(z 62 + Br Cy
(14)
which satisfies 2 = O. Analogous to the Poisson's 's equation setup, we use 100 training and 100 test examples, generated by sampling coefficients over different ranges. For the training set, the coefficients (A, B, C)) drawn from 2/|11, 1], while for the test set, they sampled from 2/(1, 2].
For both setups, we train our Neural Green's 's Function and Transolver 35) for 40 and 200 epochs, respectively. We increase the number of training epochs for Transolver to account for its slower convergence in training loss. After training, we evaluate both models on the test set by computing the relative La error across all examples. As summarized in Tab. [. Neural Green's 's Function, which

---

Poisson's is Equation
Biharmonic Equation
Train Set
Text Set
Triin Set
Test Set
Transoiver (Wu Ours
0.053
0.372
0.025
0.337
0.014
0.012
0.010
0,009
Table I: Relative L2 errors measured on the training and test sets 2D Poisson and Biharmonic equation examples. In both setups, our method achieves comparable errors the training and test sets, whereas Transolver [35 exhibits a notable gap due to its network being jointly conditioned on both domain geometry and boundary functions.
SCREWS & BOLTS
NUT
MOTOR
FITTING
GBAR
Transalver at. 35
0,221
0.320
0.407
0,180
0.281
(Wang Wang
0.239
0.372
0.528
0.259
UPT (Alkin al. g
0.358
0.516
0.765
0.392
0.507
Ours
0.189
0.275
0.338
0,160
0.243
Error Reduccion (55
14.7
14.1
16.9
10.8
13.3
Table 2: Relative L2 errors measured on the test sets of the steady-state thermal analysis benchmark. Our method achieves the lowest relative L2 error across five shape categories. Notably, it improves the metric of Transolver [35] by an average of 13.9%, despite utilizing the same network architecture. The error reduction is the difference between the two errors, divided by the second-best best error.
relies solely on features extracted from the domain geometry while remaining agnostic to the specific boundary functions used during training, generalizes well compared to Transolver 335 whose network is jointly conditioned both domain geometry and boundary functions.
5.3 Steady-State Thermal Analysis
Building the observations from Sec. 5.2 extend our experiments to a more challenging and practical task, steady- state thermal analysis complex 3D geometries. This task is appropriate for evaluating the generalizability framework to unseen problem domains, beyond the source and boundary functions considered in the previous section. To this end, we construct new PDE benchmark using the MCB dataset [16), which contains a variety of 3D mechanical part shapes. We generate tetrahedral meshes for shapes in five categories (SCREWS & BOLTS, NUT, Figure 3: Example shapes from our
dataset. Our dataset comprises diverse of unit-cube- coreecbbmooomaeeddd shapes using fTetWild [14). The mechanical part shapes from the MCB shape collection is divided into 200 shapes for training and 20 shapes for testing to evaluate generalization to unseen problem domains. We illustrate several example shapes
MOTOR, FITTING, and GeAR) by meshing the interiors
dataset (16). designed evaluate the generalizability of learning-based PDE solvers across shape variations.
in Fig. As shown, the shape collection used in the experiments includes shapes with thin structures and intricate details. Additionally, shapes within the same category can have significantly different geometries, presenting a challenge for learned solution operators to generalize effectively across diverse geometries.
To generate PDE examples, we employ an off-the-shelf FEM solver [1] to solve Poisson's 's equation defined on the shapes our dataset. These problems are constructed using different source and boundary functions derived from predefined templates, details of which are provided in the Appendix. Specifically, we use 8 source functions for training and reserve an additional for testing. For bound- ary functions, we use 2 for training and 2 for testing. This setup results in 16 unique combinations of source and boundary functions for training, while entirely unseen combinations are used for testing. It offers greater problem diversity than other benchmarks for PDEs on irregular geometries, such as ShapeNet-CFD [2,, which assumes a fixed driving velocity. Overall, the dataset comprises 3200 chape-problen pairs for training and 320 pairs for testing across all categories.
7

---

GT
Transberr
LNO
UFF
LND
LPT
Onx
(Hang Wong
(Recial
(Wary Wing
[Akancial
O Fo
Figure 4: Qualitative comparison. For each ground- truth solution in columns and six, the per per-point L2 errors from various methods are visualized in the error maps across the remaining columns. Red denotes higher values in the columns labeled "GT", while in the error maps, brighter colors indicate greater errors. Neural Green's Function demonstrates its ability to accurately predict solutions irregular geometries from diverse domains. Compared state-of-the-an neural operators trained to directly map source and boundary functions to solutions, our method achieves lower errors, indicated by the darker regions in the error maps. The examples are rendered using tetrahedral meshes, that are used solely for visualization purposes. Best viewed when zoomed-in.
we summarize the errors measured on the test sets of five shape categories using the baselines and our method. Neural Green's 's Function consistently achieves lower errors compared the baselines, demonstrating the effectiveness of the proposed method that explicitly predicts the solution operator. Compared to Transolver [35] which achieves the second-best errors across all categories, our method achieves an average error reduction of 13.9%, calculated as the difference between the best and second-bess best errors, divided by the latter. Since our framework utilizes the network architecture of Transolver (35) as backbone, this improvement highlights the significance of incorporating prior on the solution operator.
Quantitative Analysis. Tab. 2

---

GT
Ours (No Mass Reg.)
Ours
GT
Ours Mass Reg.)
Ours
000000
Figure 5: Qualitative results from the ablation study, For each ground-truti solution in columns one and four, we visualize the per-point 12 errors from variants of our method in the remaining columns. Red denotes higher values in the columns labeled "GT", while in the error maps, brighter colors indicate greater errors. The examples are rendered using tetrahedral meshes, that are used solely for visualization purposes. Best viewed when zoomed-in.
Qualitative Analysis. Qualitative results from five categories-Frriiee (rows 1-2), SCREWS
BOLTS (rows 3-5), MOTOR (rows 6-7), NUT (rows 8-10), and GEAR (rows 11 11-13) are presented in Fig- 4 Additional results can be found in the Appendix. For each example, we display color maps of the ground-truth solutions from the dataset, labeled as "GT"", alongside error maps of the predictions generated by different methods. In the solution visualizations, red indicates higher values, while in the error maps, brighter colors represent greater errors. We visualize the slices taken along one the canonical axes (z, y/, 2) of each shape for clarity. The darker error maps, compared to those of the baselines, demonstrate that method produces more accurate predictions, validating our observations from the quantitative analysis. Our method better handles irregular shapes with intricate details, such as thin structures of MOTOR shapes (rows 6-7), or the cogs of GEAR shapes (rows 113)).
Runtime Analysis. We evaluate the efficiency of our framework by comparing it with an FEM measuring the total runtimes required to process the
solver and neural operator baselines [35 34 test sets. The runtimes are summarized in Tab. 6 All reported runtimes are measured in seconds, with each sample processed sequentially to minimize system load and prevent interference that could impact the timing of other samples. Compared the FEM solver, Neural Green's Function benefits greatly from not requiring mesh generation, which accounts for the majority of the processing time in the FEM pipeline. Notably, achieves substantial reduction in runtime, being up to 350 times faster than FEM. Furthermore, Neural Green's s Function introduces only minimal computational overhead compared to neural operator baselines [35] 34 3 while delivering superior generalizability and performance, discussed previously. Details of this analysis, including the breakdown of the runtimes, are provided the Appendix.
3
Transolver (Wu et 35
LNO (Wing Wing
(s)
FEM Tonal (4)
31
(Alkin al. 3}
Ours (s)
SCrEWS & BOLTS
12,956
0.033
0.022
0.056
0.039
NUT
12,238
0.022
0,020
0,076
0.051
MOTOR
$0.095
0.090
0.088
0.166
0.140
0.012
0.059
GEAR
46.384
0.188
0.187
0.246
0.225
Table Runtime comparison against an FEM solver and neural operators. Our approach achieves up to a 350x speedup over the FEM solver, which requires mesh generation, while main- taining comparable runtime to neural operator baselines and delivering improved generalizability and performance. All runtimes are reported in seconds.
5.4 Ablation Study
We analyze the impact of mass regularization during training and the effect of feature dimensionality on model performance. In Tab. 4 we report the test errors of our model and its variant trained without mass regularization, denoted as "Qurs (No Mass Reg.)", on the SCREWS & BOLTS and GEAR categories. As reflected in the errors, regularizing the mass prediction outputs from model

---

ScReWS & BOLTS
GEAR
64
128
256
Ours (No Mass Reg.) Ours
0.285
0.411
Relative L.2
0.180
0.189
0.206
0.243
Table 5: Analysis of feature dimension. The model's performance remains consistent across different feature dimension choices.
Table 4: Quantitative results from the ablation study. Regularization of the mass predictions crucial for performance.
is crucial for performance. This is because the per-vertex masses involved in Eqn. TO are typically very small (approximately on the order of x 10-4), leading instability during the early stages of network training when the initial mass predictions deviate significantly in scale from the ground ground-truth values. Qualitative results showcased in Fig. 5 further support this observation: the model trained with mass regularization produces more accurate predictions, as evidenced by the darker regions in the error maps (columns 2-3 and 5-6).
On the other hand, we examine whether the model's performance sensitive to the dimensionality of the feature representation 4g. Intuitively, & plays role of the eigenvectors in Eqn. 7 serving as low-rank bases for approximating the operator matrix. To assess this, we train variants of our framework, which by default uses 128ddimensional feature vectors, with alternative configurations employing 64 and 256 dimensions. All models are trained on the SCREWS & BOLTS class using the same setup as in our main experiments, and the relative L2 errors on the test set are reported in Tab. [5] As shown, the model's 's performance remains insensitive to the choice of feature dimension, highlighting that even 64 learned bases are sufficient to accurately approximate the solution operator.
6 Conclusion
We present Neural Green's Function, solution operator specialized in linear PDEs whose differential operators admit eigendecomposstions. capable of generalizing across diverse domains, source and boundary functions. Our framework incorporates the principles of Green's 's functions into its design by extracting neural features solely from the geometries of problem domains. These features are then used approximate the corresponding Green's ' functions and related differential quantities, which are subsequently utilized to perform numerical integration for solution prediction. We empirically validate our design choice by applying Neural Green's 's Function solve the Poisson and Biharmonic equations on 2D domains, and further demonstrate that it outperforms -tate-of-the--e the-art neural operators on a challenging 3D benchmark encompassing a wide variety of shapes and Poisson equation instances defined over them. In particular, Neural Green's 's Function achieves a 13.9% improvement over the best performing baseline, while accelerating solution evaluation by 350x in steady -tate thermal analysis involving complex 3D geometries. These results highlight the superior generalizability of our framework across problem domains, as well as diverse source and boundary functions.
Limitations and Future Work. While we believe our approach represents a step toward developing generalizable, data-driven solution operators for PDEs, it is not without limitations. In particular, extending our framework to broader class of linear PDEs and incorporating additional boundary conditions, such as Neumann and Robin, are promising future directions. Moreover, although framework achieves runtimes comparable to existing neural operator baselines, further acceleration of the numerical integration step during the forward pass would facilitate its practical application.
Acknowledgments
This work was supported by the NRF of Korea (RS-2023-00209722); ITT grants (RS-2022- 11220594, RS- RS-2023-00227599 2025-026531 13); and the Technology Innovation Program (RS-2025-02317326 all funded by the Korean government (MSIT and MOTTE), as well as by the DRB-KAIST Fuutur Research Center.
RS-2024-003998817 RS-2025-254413111 RS-2025-25443318, RS-
10

---

References
[1] cholespy: An easily integrable cholesky solver on epu and gpu. URL httpp:://gitiiie www.e//////... -plf-
choleapy Software available GitHub.
[2] Abadi, M.. Agarwal, A., Barham. P.. Brevdo, E., Chen, Z.. Citro, C.. Corrado, G. S., Davis, A.. Dean, J., Devin, M., Ghemawat, S., Goodfellow, I., Harp, A., Irving. G., Isard, M., Jia, Y., Jozefowizz, R., Kaiser, I.,, Kudlur, M., Levenberg, Mané, D., Monga, R., Moore, Murray, D., Olah, C., Schuster, M., Shlens, J. Steiner, B., Sutskever, I,, Talwar, K., Tucker, P., Vanhoucke, V., Vasudevan, V., Viégas, F., Vinyals,
O., Warden, P.. Wattenberg, M., Wicke, M.. Yu, Y., and Zheng, X. TensorFlow: Large-scale machine learning on heterogeneous systems, 2015. URLhttps: https: /twp://wwww.co/ httpss//www.tensorfloor org/ Software available
from tensorffow.orr.
[3] Alkin, B., Fürst, A., Schmid, S., Gruber, I.,, Holzleimerr, M., and Brandstetter, J. Universal physics transformers. In NeurlPS, 2024.
[4] Barill, G., Dickson, N., Schmidt, R., Levin, I., and Jacobson, Fast winding numbers soups and clouds. ACM TOG, 2018.
[5] Boullé, N., Nakatsukasa, Y., and Townsend, Rational neural networks. NeurIPS, 2020.
[6] Boullé, N., Earls, C. J., and Townsend, A. Data-driven discovery of green's functions with human- understandable deep learning. Scientific Reports, 2022.
[7] Bradbury, J., Frostig, R., Hawkins, P., Johnson, M. J., Leary, C.,, Maclaurin, D., Necula, G., Paszke, A.,
VanderPlas, J., Wanderman-Milne, S., and Zhang, Q. JAX: composable transformations of Python+NumPy programs, 2018. URL/tttp: httttp//gittube com/ ax-ml/jax
[8] Cao, S. Choose transformer: fourier or galerkin. In NeurlPS, 2021.
Choromanski, K., Likhosherstov, Dohan, D., Song, X., Gane, A., Sarlos, Hawkins, Davis, Mohiuddin, A., Kaiser, L., Belanger, D., Colwell, L., and Weller, A. Rethinking attention with performers. In ICLR, 2021.
[10] Cooley, J. and Tukey, J. An algorithm for the machine calculation complex fourier series, Math. Comput., 1965.
[11] Deep-MI. Lapy: Toolhox for differential geometry triangle and tetrahedra meshes, 2023. URL https://gitae com//eep-
Software available on GitHub.
[12] Hao, Z., Ying, C., Wang, Z., Su, H., Dong. Y., Liu, S., Cheng, Z., Zhu, J., and Song, J. Gnot: A general neural operator transformer for operator learning. In ICML, 2023.
[13] Hu, Y., Zhou, Q., Gao, X., Jacobson, A., Zorin, D., and Panozzo, D. Tetrahedral meshing in the wild, ACM TOG, 2018.
[14] Hu, Y.. Schneider, T., Wang, B., Zorin, D.,. and Panozzo, D. Fast tetrahedral meshing in the wild. ACM TOG, 2020.
[15] Jacobson, A., Panozzo, D., et al. https://ibigll.ittuu githubb.o..
libigl:
simple C++ geometry processing library, 2018.
[16] Kim, S., Chi, H.-g., Hu, X., Huang, Q., and Ramani, K. large-scale annotated mechanical components benchmark for classification retrieval tasks with deep neural networks. ECCV, 2020.
[17] Kingma, D. and Ba, J. Adam: method for stochastic optimization. In ICLR, 2015.
[18] Kitaev, N., Kaiser. and Levskaya, A. Reformer: The efficient transformer. In ICLR, 2020.
[19] Krishnapriyan, A. Gholami, A., Zhe, S., Kirby, R., and Mahoney, M. W. Characterizing possible failure modes in physics-informed neural networks. NeurlPS, 2021.
[20] Li, Z.,, Kovachki, N. B., Azizzadenesheli, K., Liu, B., Bhattacharya, K., Stuart, A., and Anandkumar, A. Neural Operator: Graph Kernel Network for Partial Differential Equations. arXiv, 2020.
[21] Li, Z., Kovachki, N. B., Arizzadenesheli, K., Liu, B., Bhattacharya, K., Stuart, A., and Anandkumar, A. Fourier Neural Operator for Parametric Partial Differential Equations. In ICLR, 2021.
[22] Li, Z., Kovachki, N. B., Choy, C., Li, B., Kossaiti, J., Otta, S. P., Nabian, M. A.. Stadler, M., Hundt, C., Azizzadenesheli, K., and Anandkumar, A. Geometry Informed Neural Operator for Scale-Scale- 3D PDEs. In NewrIPS, 2023.

---

[23] Li. Z., Meidani, K., and Farimani, A. B. Transformer for partial differential equations' operator learning. TMLR, 2023.
[24] Li, Z., Shu, D., and Farimani, A. B. Scalable transformer for pde surrogate modeling. In NeurlPS, 2023.
[25] Li, Z.,, Yang. G., Deng, X., De Sa, C., Hariharan, B., and Marschner, S. Neural Caches for Monte Carlo Partial Differential Equation Solvers. In SIGGRAPHI Asia, 2023.
[26] Li, 7., Huang, D. Z., Liu, B., and Anandkumar, A. Fourier neural operator with learned deformations for pdes general geometries. JML.R, 2024.
[27] Nam, H. C., Berner, J., and Anandkumar, A. Solving Poisson Equations Using Neural Walk-on-Spheres. In ICML, 2024.
[28] Negi, P., Cheng, M., Krishnamurthy, M., Ying, W., and Li, S. Learning domain indeeenden green's 's func- tion for elliptic partial differential equations. Computer Methods in Applied Mechanics and Engineering, 2024.
[29] Paszke, A., Gross, S., Chintala, S., Chanan, G., Yang, E., Vitio, Z., Lin, Z., Desmaison, A., Antiga, I.,, and Lerer, A. Automatic differentiation pytorch. In NeurlPS Workshop, 2017.
[30] Raissi, M., Perdikaris, P., and Kariiadakis, G. E. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686-707, 2019.
[31] Teng, Zhang, X., Wang, Z., and Ju, Leamingg green's 's functions of linear reaction-diffusio equations with application to fast numerical solver. In Proceedings of Mathematical and Scientific Machine Learning. 2022.
[32] Umetani, N. and Bickel, Learning three-dimensional flow for interactive aerodynamic design. ACM TOG, 2018.
[33] Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L.,, Gomez, A. N., Kaiser, L. u., Polosukhin, I. Attention is all you need. In NeurlPS, 2017.
[34] Wang, T. and Wang, C. Latent neural operator for solving forward and inverse pde problems. NeurlPS, 2024.
[35] Wu, H., Luo, H., Wang. H., Wang. J., and Long, M. Transolver: transformer solver for pdes on general geometries. In ICML, 2024.
[36] Xiao, Z., Hao, Z., Lin, B., Deng, Z., and Su, H. Improved operator lcarning by orthogonal attention. In NewrIPS, 2024.
12

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: Yes, the abstract and introduction conveys the core idea of the main text.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
It is fine to include aspirational goals as motivation long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [No]
Justification: While we do not discuss the limitation in the main paper, we will include it in the revision.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g.. independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't It acknowledged judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
the paper. The authors should use their best
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [NA]
13

---

Justification:
The paper does not make theoretical claims that require proofs.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: The paper provides detailed experiment setup required to reproduce the results presented in the paper.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, No answer this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how reproduce that algorithm. (b) If the contribution is primarily
model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome deseribe the particular way they provide for reproducibility, In the case of closed-source models, it may be that access to the model is limited in some way (e.g., registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
14

---

Answer:
[No]
We plan to release the code and data required to reproduce the results upon acceptance.
Justification:
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips.cm ec/
for more details.
While encourage the release of code and data, understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https ioweeeomiisiiiioooio
//nips. //nipscc//pie wwcppbblc/gguiicuceccccc
for more details.
The authors should provide instructions data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: The paper discuss the settings of the experiments reported in the main text.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [No]
Justification: No, the reported results do not include error bars.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the enror bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
15

---

It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out range (e.g. negative enror rates).
If error bars reported in tables plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification:
We provide an runtime analysis.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary failed experiments that didn't make into the paper).
9. Code ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics http::
Answer: [Yes]
Justification:
We find no significant ethical issue while working on this work.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special cireumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., if there is special consid- eration due to laws or regulations their jurisdiction..
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [No]
Justification: While the paper does not include statements on broader impacts, we plan to include in the revision.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
16

---

The conference expects that many papers will be foundational research and not tied particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, scraped datasets)?
Question:
Answer: [NA]
Justification: We believe that the models and data used in this work have low potential of misuse.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g.. code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification:
We have included references to the original datasets used in the experiments.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (c.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the Datasets.papersstt paperswi thcode. cww/ datasets package should be provided. For popular datasets,
has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
17

---

If this information is not available online, the authors are encouraged to reach out to the asset's creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: We will release the dataset and documentation upon acceptance.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/cooe/mmoel as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsoureing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification: The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included the main paper.
According the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question:
Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification:
The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and expect authors to adhere to the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such the institution conducting the review.
16. Declaration of LLM usage
18

---

Question: Does the paper describe the usage of LLMs if it is important, original, or non-standard component of the core methods this research? Note that if the LLM used only for writing, editing, or formatting purposes and does not impact the core methodology. scientific rigorousness, or originality the research, declaration is not required.
Answer: [NA]
Justification: The core method development in this research does not involve LLMs any important, original, or non-standard components.
Guidelines:
The answer NA
that the method development in this research does involve LLMs as any important, original, or non-standard components.
Please refer to our LLM policy for what should or should be described.
19

---

Appendix
In the following, present the derivation of Eqn. Sec. A details of the problem templates
used to define the source and boundary functions for the steady-state thermal analysis experiment
(Sec. 5.3 in Sec. B details of the runtime analysis presented Sec. 5.3 Sec. C and additional
qualitative results in Sec. D
A Derivation of the Solution for the Linear System
While the derivation of Eqn. from Eqn. 5 is straightforward, we include the step-by-step derivation for completeness.
Let us decompose the solution into interior and boundary components using the section matrices K and S:
K" Uint + S"h.
(15)
where Uint (No denotes the solution Eqn. 5 within the domain.
As u satisfies Eqn. 5 we have:
L K Unt + S"h) = Mf.
(16)
Expanding the left-hand side and rearranging terms yields:
LK Un
Mf LS h.
(17)
By left-multiplying K, we obtain
KLK Uimi = KMf - KLS"
(18)
Assuming KLK is full-rank, which is ensured by imposing appropriate boundary conditions, Uint can be written as:
Uni
G (KMf
hLS"h)
(19)
where G (KLKT)-L Substituting this result into Eqn. 15 gives Eqn.
B Details of Problem Templates
When constructing the benchmark used for the steady-state thermal analysis experiment in Sec. 5.3 source functions f are populated from the following template:
[x, y, A (sin Aux cos ACTy
+ (1 - cos ATæ) (1 - sin-inn-ee ABTy)
(20)
+ sin? ADT2),
with the coefficients A = {1.25, 2.5), ;, B {1.5, 3.5). C = (1.5,335,, and D = {1.5, 3.5}. resulting in 16 different combinations. Similarly, the boundary functions are generated using the template:
h(x, y. 2) = E(x 3ry²) + F(y² - -
(21)
+ (x2 - =2),
using 4 distinct combinations of the coefficients: E = 1-1.0, and F = (0.0.1.0)
C Details of Runtime Analysis
use the binary compiled from the official implementation of fTetWild [14] for generating tetrahedral meshes. The implementation supports multi-core processing, and we used default parameters with an edge length of 0.02 for mesh generation. To populate the query points used for measuring runtimes,
In this section, we provide additional details of runtime analysis Sec. 5.3
20

---

Meshing (6)
Solve (5)
FEM Tocal (4)
(9)
Churs Forward (9)
Ours Total (S)
SCREWS & BOLTS
12.579
0.577
12.956
0,016
0.025
0.039
Nurr
11.601
0.637
12.238
0.014
0.038
0.051
MOTOR
1.434
50.095
0.082
0.140
GEAR
45,803
0.581
46.384
0.180
0.044
0.225
Table 6: Breakdown runtime comparison with the FEM solver. The total runtime of the FEM solver includes both the meshing and linear solve times, whereas our runtime is computed as the sum of the time required for the interior point query (denoted IPQ) and the network forward pass. All runtimes reported in seconds.
Transober eWo all (35)
1ND
1PP
(Wing & Wang M
(Alkin al. 3
Forward (s)
(s)
(s)
(6)
Forward (s)
(5)
(s)
0.033
0.022
NUT
0.014
0.008
0.022
0.006
0.020
0.062
0.076
0.038
0.051
MOTOR
0,082
0.088
0.090
0.006
0.088
0,084
0.156
0.058
0,140
FITTING
0,024
0.008
0.032
0.006
0.030
0:052
0.076
0.035
0.059
0.180
0.225
Table 7: Breakdown of runtime comparison with neural operators. The total runtime of each neural operator computed as the sum of the time required for the interior point query (denoted IPQ) and the network forward pass. All runtimes are reported in seconds.
uniformly sample the 3D space and determine whether each point lies inside the domain using the Fast Winding Number [4] implemented in libigl [5]. To ensure fair comparison with our method, re-implemented the solver algorithm from [1 to support GPU acceleration. In particular, we utilize Cholespy[I [. GPU-accelerated
which leverages GPU acceleration during network inference,
Cholesky solver, for matrix prefactorization to solve Eqn. 5 The runtime breakdowns corresponding to the results in Tab. 3 are presented in Tab. and 7 respectively, each of which summarizing the time spent in each stage of the FEM solver- -inlurninlluding meshing and solving linear systems-at well as, for the neural operators, the time required for interior point queries (IPQ) and network forward passes, All runtimes are measured on system equipped with an Intel Xcon Gold 6442Y processor with 24 cores and an NVIDIA RTX 3090 GPU with 24 GB of VRAM.
21

---

D
Additional Qualitative Results
GT
Transolver
LNO
UPT
Ours
GT
Transolver
LNO
UPP
Ours
HOHT
Figure 6: Additional qualitative comparison. For each ground-truth solution in columns one and six, the per-point Lg errors from various methods visualized in the error maps across the remaining columns.
22


Reviews:
[{"review_id": "VFt96MQ1Vq", "reviewer": "Reviewer_SXcQ", "summary": "This paper presents a method that learns to predict neural approximate Green's functions for Poisson's equation. The approximate Green's function is a low-rank approximation formed from a learned domain-dependent pointwise feature map. The predicted Green's function yields an approximate solution operator that generalizes under changes in domain shape, boundary conditions, and source function.", "questions": "- What do you mean by \"leveraging the eigendecomposition of Green's functions [by] reconstructing them from features extracted from problem domains\" (111-112)? Are the features the eigenvalues/eigenfunctions? Or some approximation thereof?\n- Does your method generalize beyond Dirichlet boundary conditions?\n- How does your method compare to boundary element / boundary integral methods? Or Monte Carlo/Walk-on-spheres methods, which also support pointwise queries without requiring mesh generation or a global linear solve? When is it advantageous to learn a proxy solution operator vs., for example, using (possibly denoised) Monte Carlo estimates?\n- The method represents the geometry of the domain through the sampled query points. Thus the query points are doing double duty as points at which the Green's function and solution are evaluated, and as an input that conditions the predicted solution. How much does the learned solution depend on the details of this sampling? Ideally, the predicted solution at a fixed point $p$ in the domain should be invariant to the choice of other query points $Q$. Does the choice of query points matter at test time? Only during training? Or are the results nearly invariant to these choices, as one would hope?", "limitations": "N/A", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 2, "strengths_and_weaknesses": "The results seem pretty impressive. However, the evaluation is on only a limited number of shapes (test set of 20). Since the main point of the paper is generalization to arbitrary shapes, it might be good to stress-test the method on more shapes (e.g., from a completely different collection or category). I am also not sure what you mean in composing combinations of domain, boundary function, and source function. E.g., a given boundary function is defined on the boundary of a particular domain. What does it mean to test the same boundary function across two different domains?\n\nIt might be good to evaluate the results quantitatively using the H^1 metric, rather than just L^2 as in equation (12). Given a source function $f$ in $L^2$, the solution to Poisson's equation should lie at least in H^1 in the interior, and it would be good to ensure that not only do the pointwise values of your approximate solution match the ground truth well, but also that their derivatives match well.\n\nTo give contex to your runtime comparison in Table 2, it would be fair to note also the total time used for training. (Though I understand the point that *at test time* for an individual domain, the neural proxy is much faster.)", "was_revised": true, "final_justification": "The authors' rebuttal addressed most of my points. I also see that they did an experiment with a biharmonic equation. I would still like to see evaluation of solution error in $H^1$ rather than just $L^2$.", "rebuttal": "Dear Reviewer SXcQ,\n\nThank you for your positive feedback and for finding our results impressive. We appreciate your constructive comments, which have helped us improve our work. Below, we address your questions in detail:\n\n**Generalization Across Shape Categories and Boundary Conditions**\n\nWe agree that the practicality of Neural Green’s Function would be improved if a model trained on one shape category could generalize across different categories. In this work, similar to our baselines [1,2] that learn solution operators over geometries within a specific category (e.g., cars, airfoils), we focused on a relaxed setting where all shapes belong to the same category but exhibit intra-category variations, as shown in Fig. 5 of the main paper. We acknowledge that enabling inter-category generalization would substantially broaden the method’s applicability, and we plan to explore this direction in future work. Similarly, we agree that extending the method to support diverse boundary conditions is an important and valuable direction for future work. We plan to explore this in more depth, and in the upcoming revision, we will clarify the problem statement to more accurately reflect the current scope and limitations of our approach. \n\n**Training Time**\n\nUsing the *Screws and Bolts* class from the MCB-B dataset, which contains point clouds with an average of 8.6K points, our model takes approximately 15 minutes to process 3,200 examples per training epoch. On the *Gear* class, with a higher average of 19.8K points, a single epoch takes about 30 minutes. We also evaluated Transolver with a comparable number of parameters and observed similar runtimes across both classes. Since all models are trained for 40 epochs in our experiments, the total training time amounts to approximately 10 hours for the *Screws and Bolts* category and 20 hours for the *Gear* category. We will include a discussion of the runtime in the revised manuscript.\n\n**Clarification of L111-112**\n\nWe apologize for the ambiguity in our previous explanation. Our model takes a point cloud as input and is trained to approximate the eigenvalues and eigenfunctions used to construct the operator defined in Eqn. 8 of the main paper. At inference time, our model does not require performing any eigendecomposition of operators defined over the domain. We will clarify this point in the upcoming revision.\n\n**Comparisons Against Boundary Element/Integral Methods and Monte-Carlo Methods**\n\nThank you for bringing this to our attention and giving us the opportunity to discuss mesh-less PDE solvers that are closely related to our work. We believe each approach has its own strengths, and we would like to note that Neural Green’s Function is designed to complement existing solvers rather than replace them.\n\nStarting with the Boundary Element Method (BEM), to the best of our knowledge, it requires factorizing a dense system upfront. While this enables efficient solves for multiple boundary conditions on the same domain, the factorization must be recomputed when the domain changes. Monte Carlo and Walk-on-Spheres methods, which have recently gained attention in the graphics community, avoid both meshing and matrix factorization by estimating solutions through pointwise random walks. However, these methods typically require a large number of samples to achieve low-variance estimates, resulting in longer runtimes compared to learned proxy operators like ours, which produce solution estimates in a single forward pass.\n\nAlthough our method requires an initial training phase, it enables fast inference through a single forward pass—a key distinction from the aforementioned methods. In real-world applications, we envision a hybrid workflow where Neural Green’s Function is used for rapid prototyping or early-stage exploration, providing candidate solutions that can be further refined using more precise solvers such as BEM- or WoS-based approaches.\n\n**Impact of Choosing Query Points**\n\nOur analysis conducted on the test set, summarized in Tab. 1, implies that the model can handle point clouds of varying sizes—and consequently, densities—since all shapes are normalized to the unit cube. However, as our approach is data-driven, excessively dense or sparse query point clouds may present challenges. A practical workaround is to match the density of query points to that of the training data, either by selecting a representative subset that captures the overall shape (if input is dense) or by sampling additional interior points using uniform sampling combined with winding-number-based rejection (if input is sparse).\n\n**Table 1. Statistics of Point Cloud Sizes**\n\n| Category          | $V_{\\text{avg}}$    | $V_{\\text{std}}$    | $V_{\\text{min}}$ | $V_{\\text{max}}$  |\n|-------------------|----------|----------|--------|--------|\n| *Nut*               | 18086.10 | 7492.09  | 6849   | 34189  |\n| *Screws and Bolts*  | 8666.00  | 5836.00  | 1733   | 21679  |\n| *Gear*              | 19790.30 | 11614.31 | 5051   | 49634  |\n| *Fitting*           | 13962.35 | 3854.70  | 4978   | 20406  |\n| *Motor*             | 27106.25 | 10352.99 | 9942   | 49454  |\n\n**Composing Combinations of Domain, Boundary Function, and Source Function**\n\nBy using the same boundary function across different domains, we meant that the coefficients of the template function—evaluated at the boundary points of each domain—are kept identical. We apologize for the confusion and will revise the text to clarify this aspect of our experimental setup.\n\n**Reference**\n\n[1] Transolver: A Fast Transformer Solver for PDEs on General Geometries: https://arxiv.org/abs/2402.02366  \n[2] Universal Physics Transformers: A Framework For Efficiently Scaling Neural Operators: https://arxiv.org/abs/2402.12365"}, {"review_id": "q3bAWvb9fQ", "reviewer": "Reviewer_fcRy", "summary": "The paper introduces Neural Green's Function (NGF), a mesh-free neural operator for linear Poisson PDEs with Dirichlet boundary conditions. NGF predicts a low-rank approximation of the Green's kernel through spectral feature learning, uses geometry-only conditioning for inference independence, and employs auxiliary MLP heads for solution evaluation. Tested on 3,400 heat-diffusion problems, NGF achieves 14% lower error than Transformer baselines and 350× speedup over FEM.", "questions": "1. What value of feature dimension d is used? How sensitive are results to it?\n\n2. Does the model ever predict negative or near-zero eigenvalues causing numerical issues?\n\n3. How are interior query points chosen during training and inference? Uniform sampling? Poisson-disk?\n\n4. Could NGF be combined with reduced-order models (e.g., POD) to further cut costs?\n\n5. Are there cases where NGF underperforms (e.g., extremely thin shells)? A qualitative failure example would be instructive.\n\n6. How difficult would extension be to Neumann/mixed boundary conditions or other linear elliptic operators (bi-Laplacian, anisotropic Laplacian)? \n\n7. For CAD models containing multiple disjoint solids, would spectral factorization remain stable?\n\n8. Since $G_{\\theta}$ is built from learned features, how do authors guarantee $KLK^{\\top}$ invertibility and positive mass prediction (beyond Soft-plus non-negativity)?\n\n9. Would performance degrade when evaluating on denser grids than training grids?\n\n---\n\nOverall, This paper presents a reasonable operator-learning framework that applies classical Green’s function analysis to solve Poisson’s equations. While the method is clearly explained and the results are promising on the presented benchmarks, the overall contribution is incremental rather than groundbreaking. Some aspects, such as the limitations to specific PDE types, lack of uncertainty quantification, and limited discussion of generalizability, prevent a higher recommendation. With further empirical validation and broader evaluation, the impact and significance of the work could be strengthened. At its current stage, I consider the contribution modest but relevant for the NeurIPS community.", "limitations": "The paper does not explicitly articulate its limitations, but the proposed method is confined to specific types of partial differential equations (PDEs) and lacks uncertainty quantification, offering only limited discussion on generalizability.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "**Strengths:** Unlike prior neural-operator works that learn a black-box mapping, their method explicitly reconstructs the Green’s kernel from geometry alone and exploits linear superposition. Demonstrates significant speed improvements compared to TetWild + Cholesky method. The paper is well organized as well. \n\n**Weaknesses:** Currently limited to Poisson + Dirichlet cases. Framework assumes single connected components. No discussion of stability/positive-definiteness. Interior-point sampling density is fixed (unspecified).", "was_revised": true, "final_justification": "I am quite positive about this paper and believe it would make a valuable contribution to the field upon acceptance.", "rebuttal": "Dear Reviewer fcRy,\n\nWe appreciate your recognition of our approach—explicitly reconstructing the Green’s kernel from shapes—as a strength compared to existing learning-based solvers that treat the solution operator as a black-box mapping. We believe this direction can contribute to the development of more interpretable learning-based solvers. Below, we have reviewed your questions and compiled our responses accordingly.\n\n**Extension to Other PDE Problems**\n\nWe extended our framework to the biharmonic equation, based on the observation that our formulation is applicable to PDEs whose differential operators admit an eigendecomposition. We empirically verified that our method can approximate the corresponding solution operator using a toy example. The detailed experimental setup and results are provided below.\n\nWe consider a triangular mesh obtained by discretizing the unit square [0, 1] x [0, 1]. Similarly to our experiment setup involving shapes from the MCB-B dataset, we define a template function of form\n$$\nf(x, y, z) = A(x^4 - 6x^2 y^2 +y^4) + B(x^4 - 6x^2z^2 + z^4) + C(y^4 - 6y^2z^2 +z^4).\n$$\nWe use this template function—whose Bilaplacian is zero—to generate our training and test samples by varying the coefficients $A$, $B$, and $C$ of boundary functions. For the training set, all coefficients are sampled from $\\mathcal{U}[-1, 1]$, while for the test set, they are sampled from $\\mathcal{U}[1, 2]$. We generated 100 samples each for training and testing. Along with the corresponding triangular mesh, these samples are used to train both Neural Green’s Function and Transolver. During training, we observed that Transolver required a larger number of epochs to reach comparable accuracy on the training set. Accordingly, we trained our model for 40 epochs, while the Transolver was trained for 200 epochs. \n\nAfter training, we evaluated both models on the test set by measuring the relative $\\mathcal{L}_2$ errors, as summarized in the table below. We note that both Transolver and Neural Green’s Function converge on the training set, as evidenced by their low training errors. However, Transolver shows a decline in performance on unseen examples compared to the training set, which may be attributed to the increased difficulty of generalizing to unseen boundary functions. This is due to the joint conditioning over both geometry and boundary function values. In contrast, our method remains agnostic to specific boundary functions and relies solely on geometric features, achieving similar accuracy on the test set. We will incorporate these results in future updates and revise the Methods section accordingly. While the presented experiment demonstrates the potential of our framework to generalize to a broader class of PDEs, handling diverse boundary conditions remains an open challenge. We plan to explore extensions of our framework to address this in future work.\n\n| Method            | Transolver | Ours   |\n|-------------------|------------|--------|\n| Relative $\\mathcal{L}_2$ (Train) | 0.0252     | 0.0102 |\n| Relative $\\mathcal{L}_2$ (Test)  | 0.3373     | 0.0093 |\n\n**Extension to Shapes with Multiple Components and Stability**\n\nWe empirically observed that our model remains effective when applied to shapes with multiple disconnected components. This robustness is attributed to the model’s ability to directly approximate the solution operator while capturing point interactions through dot products between feature vectors. Specifically, we analyzed the *Motor* class in the MCB-B dataset and identified 15 shapes with more than two disconnected components. Using our pre-trained model, we predicted solutions for PDE problems generated with the same coefficient configurations as the test set. Across 240 examples (15 shapes × 16 problems each), the model achieved a relative $\\mathcal{L}_2$ error of 0.278. \n\n**Discussion of Numerical Stability**\n\nSince our framework does not require explicit prediction of eigenvalues to compute the terms in Eqn. 10, we have not encountered issues related to negative or near-zero eigenvalue predictions. Similarly, since the model directly approximates the inverse operator of $K L K^T$—which is assumed to exist in our problem setup and during dataset construction—no explicit matrix inversion is required during inference. Accordingly, we have not encountered any issues related to invertibility during training or inference.\n\n**Dimensionality of Feature Dimension**\n\nWe agree that examining the sensitivity of our model to the feature dimension is important, especially given that, as you correctly noted, our approach can be interpreted as learning a low-rank approximation of the Green’s function. By default, we used 128-dimensional feature vectors. To assess the impact of feature dimensionality, we trained models with 64 and 256 dimensions on the *Screws and Bolts* class, using the same training setup for 40 epochs. As summarized in the table below, the model’s performance remains insensitive to the choice of feature dimension. We will include this analysis in the upcoming revision.\n\n| Feature Dimension | 64   | 128   | 256   |\n|-------------------|------|-------|-------|\n| Relative L2       | 0.18 | 0.189 | 0.206 |\n\n**Sampling Method for Interior Query Points**\n\nFor computing interior points, we first perform uniform sampling within the bounding box—centered at the surface mesh’s centroid—with a slight margin of 10%, and use fast winding number computation to identify interior points. Since our approach is data-driven, excessively dense or sparse query point clouds may present challenges. A practical workaround is to match the density of query points to that of the training data, either by selecting a representative subset that captures the overall shape or by sampling additional interior points using uniform sampling combined with winding-number-based rejection.\n\n**Failure Cases of Neural Green’s Function**\n\nWe apologize for not including a discussion of failure cases in the paper. As you rightly pointed out, thin regions in the input geometry can pose challenges during inference. While the rebuttal policy prohibits uploading visuals to OpenReview, we will include a detailed discussion of such failure cases, along with qualitative examples, in the revised version.\n\n**Limitations**\n\nWe appreciate you for pointing this out. While we have demonstrated that our framework extends beyond the originally considered Poisson’s equation (e.g., Biharmonic equation), we acknowledge that its applicability could be further improved by supporting a broader range of boundary conditions. We will discuss the limitations of our work, including failure cases, in the upcoming revision."}, {"review_id": "1KZA6Eqq7d", "reviewer": "Reviewer_KHba", "summary": "The authors introduce an operator learning framework - Neural Green's function. The idea is to embed Green's functions in the learning procedure for operators to enable data-efficient learning and deal with the issue of poor generalization in unseen domains, boundary functions and source terms. The authors propose approximating the Green’s function for a given problem domain allowing to compute solutions for arbitrary source and boundary functions by numerically evaluating the integral. They demonstrate superior generalization and speed-ups compared to other state-of-the-art approaches on the MCB dataset for Poisson equation.", "questions": "1) Figures: I assume the authors have used the same color-bar for all the error plots in Figure 1. A color-bar for the true solution (column 2) and errors (columns 3-6) would be helpful to the reader. The same applies to Figures 3 and 4 in the main text and Figure 6 in Appendix D. \n2) Lines 41-42: I do not think it is quite true that neural PDE solvers generally match the accuracy of classical numerical methods. Optimization difficulties can often make it very hard for the neural PDE solvers to reach the accuracy of numerical solvers. I think this sentence needs to be rephrased. Could you elaborate your view? \n3) The authors compare the FEM run-time with the Neural Green's function evaluation/inference time, not the training time. I think it is perfectly fine to do so to stress the advantage during inference. However, it is also important to mention the training times and the fact that the major computational load for neural operators is offloaded to the training phase. Hence, the speed-ups during inference do not come without a price. Moreover, the operator frameworks are trained on data and thus have an advantage that they have more information, compared to solving a PDE for any particular instance with an FEM solver, which does not leverage any data from other simulations.  Could you elaborate your view and why you do not mention the training times?", "limitations": "Not mentioned at all, no separate limitation section.", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 2, "originality": 2, "strengths_and_weaknesses": "Strengths:\n1) Approach: The central idea of incorporating Green’s functions into the framework design is quite interesting. Neural green's function is agnostic to source and boundary functions used for training by design and seems to generalize well across diverse geometries. \n2) Empirical validation for the Poisson equation: For the Poisson equation, the baseline comparisons with other state-of-the-art architectures and approaches, different source and boundary terms, and different irregular shapes are strong. \n3) Writing: The paper is well-written, and is relatively easy to follow. \n\nWeaknesses:\n1) Reproducibility: The code and data are unavailable. The details mentioned in the paper are certainly not enough to reproduce the results. For instance, the architectural details (number of layers, width) are missing, along with many other parameters such as number of interior and boundary points, activation functions, etc. \n2) Lack of examples: The study is restricted to a single PDE - Poisson equation. While the Poisson equation is explored in depth, an additional PDE like the Helmholtz equation (for instance, from reference [30] in the paper) would strengthen the results further. \n3) Limitations: The limitations of the proposed approach are not mentioned at all. I believe it is quite important to mention the limitations and key directions for future work that the authors feel are important. \n4) Weak baselines: the FEM solver is \"off the shelf\" (l.222), so not necessarily state of the art. Comparing the inference time of a (trained) neural operator in specific examples to an arbitary (i.e., not specifically designed) FEM solver is not reasonable. Poisson's equation is very well known in the scientific community and can be solved extremely efficiently (e.g. using multigrid). An off-the-shelf FEM solver probably cannot compare to a solver designed for this problem. It is ok to use it to generate the training data, but claiming superior performance in terms of inference time is not reasonable (but a very common mistake in physics-informed ML).", "was_revised": true, "final_justification": "I appreciate the rebuttal of the authors, and working on a new example. It is very difficult to judge how the updated version (or the code) look like, and the additional statement on the solution being an \"initial guess\" for subsequent numerical solvers also cannot be evaluated in short time, so I will not reconsider my score.", "rebuttal": "Dear Reviewer KHba,\n\nThank you for recognizing the integration of Green’s functions into our framework as a key strength. We also appreciate your positive remarks on the clarity of our writing. Below, we address your comments and summarize our responses accordingly:\n\n**Implementation Details for Reproducibility**\n\nWe promise to release the full codebase (including model training and dataset preparation) along with the datasets upon acceptance to support reproducibility and facilitate future research. The network architecture consists of an initial MLP that maps input features (XYZ coordinates and, if applicable, sampled function values) to latent representations. This is followed by eight Transolver blocks and a final MLP that decodes the latent features into either the solution function or, in our case, features used for Green’s function approximation. We use a consistent feature dimension of 128 across all layers.\n\n**Lack of Examples**\n\nWe demonstrate that the proposed framework can be extended to other PDEs, such as the biharmonic equation, as long as their differential operators admit an eigendecomposition. To verify this empirically, we prepared a toy setup of the biharmonic equation and compared the performance of our model against Transolver, the strongest baseline in our experiments. The detailed experimental setup and results are described below.\n\nWe consider a triangular mesh obtained by discretizing the unit square [0, 1] x [0, 1]. Similarly to our experiment setup involving shapes from the MCB-B dataset, we define a template function of form\n$$\nf(x, y, z) = A(x^4 - 6x^2 y^2 +y^4) + B(x^4 - 6x^2z^2 + z^4) + C(y^4 - 6y^2z^2 +z^4).\n$$\nWe use this template function—whose Bilaplacian is zero—to generate our training and test samples by varying the coefficients $A$, $B$, and $C$. For the training set, all coefficients are sampled from $\\mathcal{U}[-1, 1]$, while for the test set, they are sampled from $\\mathcal{U}[1, 2]$. We generated 100 samples each for training and testing. Along with the corresponding triangular mesh, these samples are used to train both Neural Green’s Function and Transolver. During training, we observed that Transolver required a larger number of epochs to reach comparable accuracy on the training set. Accordingly, we trained our model for 40 epochs, while the Transolver was trained for 200 epochs.\n\nThe following table summarizes the training and test accuracies of both models after convergence. While Transolver and our model achieve comparable accuracy on the training set, our approach—being agnostic to the specific boundary and source functions used during training—demonstrates stronger generalization to unseen problem instances. We will incorporate these results in future updates and revise the Methods section accordingly.\n\n| Method              | Transolver | Ours   |\n|---------------------|------------|--------|\n| Relative $\\mathcal{L}_2$ (Train) | 0.0252     | 0.0102 |\n| Relative $\\mathcal{L}_2$ (Test)  | 0.3373     | 0.0093 |\n\n**Limitations**\n\nWhile we have demonstrated that our framework extends beyond the originally considered Poisson’s equation, we acknowledge that its applicability could be further improved by supporting a broader range of boundary conditions. We will include a discussion of these limitations in the upcoming revision.\n\n**Clarification on L41-42.** We agree with your point and appreciate your thoughtful feedback, which allows us to clarify and correct our statement. Similar to other neural operator approaches, the primary goal of Neural Green’s Function is to provide an efficient approximation of numerical solvers—potentially useful during the early stages of exploration, before resorting to more precise but more computationally demanding methods. We will revise the relevant text accordingly in the upcoming version.\n\n**Training Time of Learning-Based Solvers** We agree with your point that the acceleration provided by learning-based solvers comes at the cost of training time, making it an important factor to consider when evaluating such methods. In our experiments on the Screws and Bolts and Gear classes from the MCB-B dataset—each containing point clouds with average sizes of 8.6K and 19.8K points, respectively—both our method and Transolver require approximately 15 and 30 minutes, respectively, to process 3,200 examples per training epoch. Since all models are trained for 40 epochs in our experiments, the total training time amounts to approximately 10 hours for the *Screws and Bolts* category and 20 hours for the *Gear* category.\n\n**Clarification on Baseline Setup**\n\nWe would like to clarify the FEM baseline setup, particularly regarding the runtime evaluation. As detailed in Appendix C, the baseline was constructed using a state-of-the-art tetrahedralization algorithm [1] and a GPU-accelerated linear solver based on Cholesky decomposition [2], to ensure fair comparison. We will include a more detailed description of the FEM baseline in the main paper in a future revision.\n\n**Color-Bar in Figures.** Thank you for notifying us. As you correctly pointed out, the color bars for the ground-truth solutions (second column) and the error maps (third to sixth columns) in Fig. 1 are shared within each example (i.e., per row). The same convention applies to other visualizations, such as Figs. 3, 4, and 6. We will revise the figures to make this clearer in the upcoming revision.\n\n**Reference**\n\n[1] Fast Tetrahedral Meshing in the Wild (ACM ToG 2020): https://dl.acm.org/doi/10.1145/3386569.3392385  \n[2] cholespy: An easily integrable Cholesky solver on CPU and GPU: https://github.com/rgl-epfl/cholespy"}, {"review_id": "vZ56UUxKiY", "reviewer": "Reviewer_FSUU", "summary": "This paper offers a neural network based approach to solving Poissons equation with an approach that leverages the computational form of classical linear algebra solution. The main insight is the Shur inverse, the inverse of the interior of the Laplacian, analogous to the discrete Green's function, is approximated with a neural network. This is only possible because of the linearity of the Poisson problem.\nFollowing this the authors make the key insight that the Greens function is only geometry and not boundary value or source term dependent. Given this insight the problem is reduced to predicting Eigen-value scaled Eigen-vectors, a component from the decomposition of the Greens function.\n\nTo ensure the method is fully mesh independent this predicted quantity is also used to predict M and KLS.T (the surface to boundary sub-matrix of the Laplacian). \n\nExperiments are performed on 5 complex, non-convex geometries with Neural Green beating Transovler, the next best method, by >10% error reduction, with errors against an FEM approach in the 16-34% range but with a 350 times speed up.", "questions": "What kind of neural network is used to predict $\\Phi$? In the theme of adding \"stronger priors\", would it be possible to make a stronger inductive bias for example roto-translational equvariance in this network?\n\nDid the authors try with any exotic feature engineering, ie signed distance function from geometry instead of positional coordinate?\n\nI understand equation (7) inspires the dot products between features, but in (6) G involves the inverse of K.L.K^T. Then how is it justified to use non-verted K/K^T multiplication?\n\nIf I could share for idea for the authors to be able to leverage their Poisson specialist method to other PDEs is taking an approach similar to\nZeng et al - \"PhyMPGN: Physics-encoded Message Passing Graph Network for spatiotemporal PDE systems\", where the Laplace Block calves out the Laplacian solve, deferring other dynamics contributions from parallel NN architectures, GNN (or a NERF if requiring mesh free property).", "limitations": "Yes", "rating": 4, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 2, "originality": 2, "strengths_and_weaknesses": "Formulating the solution operator closely to the linear alegra computational structure is a strong mathematically grounded approach. Allowing to use a stronger inductive bias, utilising the linearity of Poissons equation and the stencil of the numerical method, replacing the computationally complex parts (matrix decomposition) with NN acceleration.\n\nThe work is clear, well presented and self contained. The benchmark datasets of the Poisson solves on the MCB dataset is a valuable contribution to the community.\n\nHowever I have some concerns to the value of the work. A 13.9% improvement is not huge, especially since the errors against an FEM approach still in the 16-34% range. So it's a 1-2% actual error improvement. Moreover the Transolver baseline can be applied to PDEs other than the linear Poisson equation.\n\nRuntime comparison - the paper compares against FEM runtime (which is good) but don't include the runtime of the other ML baselines.\nThere might be be a sitution where Transolver has a 22.1% error against FEM but a x1000 speedup versus Neural Greens that has a 18.9% error but a x330 speedup.\n\nThe authors critise the transformer style baselines saying \"but they require large-scale training datasets to achieve generalization.\" \nbut I couldn't see any analysis that they overcome the need to for large datasets or any evidence that there method generalizes better than the baselines beyond the 200/20 train/test split.\n\nAlthough the method is novel, the paper is clear, well explained, my prevailing concerns about the benefits compared to existing neural solvers leave my rating at weak reject unless the authors can convince me on the evaluation and usefulness of this work.", "was_revised": true, "final_justification": "The authors clarifications, new experiments and including solve times have addressed my concerns.", "rebuttal": "Dear Reviewer FSUU,\n\nWe sincerely appreciate your detailed review and positive feedback, especially your recognition of our work as a “strong, mathematically grounded approach” presented in a clear and self-contained manner. We have reviewed your questions and summarized our responses below.\n\n**Concerns on Effectiveness & Applicability of Neural Green’s Function**\n\n**Effectiveness.** We would like to highlight that similar magnitudes of absolute improvement (on the order of 1E-3 to 1E-2) were also reported by Transolver, one of the most powerful neural operators to date, on the ShapeNet-CFD dataset, yet were interpreted as meaningful thanks to their substantial relative gains (approximately 20–30%). Following this, we reported relative performance improvements to better reflect the practical significance of our results.\n\n**Application to Other PDEs and Unique Capabilities.** We agree that other neural operators, such as Transolver, are capable of handling a broader range of PDEs. However, we believe that Neural Green’s Function offers distinct advantages due to its design, which closely mirrors the structure of analytical solution operators. In particular, it is independent of the specific test functions used during dataset generation and relies solely on geometry: a property exclusive to ours that we find crucial for generalizing to unseen source and boundary functions. To better support this claim and further demonstrate the generalizability of our framework, we conducted a toy experiment on Poisson’s equation using a fixed geometry with an extended set of boundary function examples.\n\nSpecifically, we consider a triangular mesh obtained by discretizing the unit square [0, 1] x [0, 1]. The problem instances were generated using the same template function as in our MCB-B experiments, with the following form used as the boundary condition:\n$$\nf(x, y, z) = A (x^3 - 3xy^2) + B(y^3 -3x^2y) + (x^2 - z^2).\n$$\n\nWe generated 100 training and 100 test examples by sampling tuples $(A, B)$. For the training set, all coefficients are sampled from $\\mathcal{U}[-1, 1]$, while for the test set, they are sampled from $\\mathcal{U}[1,2]$. Each sample, along with the associated triangular mesh, was used to train both Neural Green’s Function and Transolver. During training, we observed that Transolver required a larger number of epochs to reach comparable accuracy on the training set. Accordingly, we trained our model for 40 epochs, while the Transolver was trained for 200 epochs.\n\nAfter training, we evaluated both models on the test set by measuring the relative $\\mathcal{L}_2$ errors, as summarized in the table below. We note that both Transolver and Neural Green’s Function converge on the training set, as evidenced by their low training errors. However, Transolver shows a decline in performance on unseen examples compared to the training set, which may be attributed to the increased difficulty of generalizing to unseen boundary functions. This is due to the joint conditioning over both geometry and boundary function values. In contrast, our method remains agnostic to specific boundary functions and relies solely on geometric features, achieving similar accuracy on the test set.\n\n| Method              | Transolver | Ours   |\n|---------------------|------------|--------|\n| Relative $\\mathcal{L}_2$ (Train) | 0.0529     | 0.0139 |\n| Relative $\\mathcal{L}_2$ (Test)  | 0.3716     | 0.0124 |\n\n**Application to Other PDEs.** We demonstrate that the proposed framework can be extended to other PDEs, such as the biharmonic equation, as long as their differential operators admit an eigendecomposition. Following the previous experiment setup, we constructed a toy example in which both Transolver and our method were tasked with learning the solution operator for the biharmonic equation on a 2D domain. For the biharmonic cases, we consider the following template function\n$$\nf(x, y, z) = A(x^4 - 6x^2 y^2 +y^4) + B(x^4 - 6x^2z^2 + z^4) + C(y^4 - 6y^2z^2 +z^4),\n$$\nwhose bilaplacian ise zero. Following the setup of the previous Poisson experiment, we generated 100 training and 100 test samples by varying the coefficients $A$, $B$, and $C$ of the boundary functions. For the training set, coefficients were sampled from $\\mathcal{U}[-1, 1]$, while for the test set, they were sampled from $\\mathcal{U}[1, 2]$. Both Transolver and our method were trained on the training split—Transolver for 200 epochs and our model for 40 epochs—and evaluated on the test set after convergence. The relative $\\mathcal{L}_2$ errors are summarized in the table below. Our Neural Green’s Function demonstrates consistent generalizability in this new PDE class as well, in contrast to Transolver, whose network is jointly conditioned on domain geometry and function examples.\n\n| Method            | Transolver | Ours   |\n|-------------------|------------|--------|\n| Relative $\\mathcal{L}_2$ (Train) | 0.0252     | 0.0102 |\n| Relative $\\mathcal{L}_2$ (Test)  | 0.3373     | 0.0093 |\n\n**Runtime Comparison**\n\nWe agree that the runtime of neural solvers is an important factor for a thorough analysis. Accordingly, we measured the runtimes of all learning-based baselines on the test set used in our main experiments. As shown in the table below, these methods achieve runtimes on the order of fractions of a second, including both interior point querying and the network forward pass. While Transolver achieves the best runtime performance, please note that UPT and Neural Green's Function are already 280 times faster than a numerical solver, enabling new applications such as real-time solution preview.\n\n| Method    | FEM | Transolver | UPT  | Ours |\n|-----------|------|------------|------|------|\n| Time (s)  | 28.02 | 0.07       | 0.12 | 0.10 |\n\n**Clarifications**\n\n**Network Architecture.** As noted in Section 5.1, our network is based on the Transolver architecture; however, we acknowledge that the provided details were insufficient. Specifically, we use the same network as the Transolver baseline, with only minor modifications—namely, the addition of linear layers to transform the extracted features into the quantities used in Eqn. 10 of the main paper. (e.g., $\\Phi_{\\theta}$, $\\Psi_{\\theta}$, and $\\mathbf{M}_{\\theta}$). We will include these details in the upcoming revision.\n\n**Additional Feature Engineering.** Thank you for your suggestion. While we have not yet experimented with incorporating additional geometric information beyond raw point clouds, we agree that leveraging more sophisticated features (e.g, SDFs) could further improve model performance. We appreciate your input and will explore this direction in future work. \n\n**Use of Non-Inverted K/K^T Multiplications.** In our framework, the solution operator $\\mathbf{G}$ is directly approximated by $\\mathbf{G}_{\\theta}$, which is parameterized as dot products between learned features and can be directly used in Eqn. 10 to compute solutions, without matrix inversion. The multiplication of $\\mathbf{K}$ and its transpose in Eqn. 8 serves to index the features corresponding to interior points, ensuring that the resulting matrix has the appropriate shape in $\\mathbb{R}^{(N_v - N_b) \\times (N_v - N_b)}$.\n\nLastly, we appreciate your recommendation of an inspiring work (PhyMPGN), which seems highly relevant for advancing our approach."}]

Author rebuttal:

Meta-review: This paper presents the Neural Green's Function (NGF), a neural operator for linear Poisson PDEs with Dirichlet boundary conditions. The authors demonstrate that NGF enables data-efficient learning, improves generalization across unseen domains and boundary/source functions, and achieves significant speed-ups compared to state-of-the-art approaches.

The majority of reviewers recommended acceptance, highlighting the paper's strong mathematical grounding, clear exposition, and compelling empirical results. The authors also effectively showcase NGF's utility beyond the Poisson equation, making a strong case for its value.

One reviewer (KHba) raised several concerns that I believe the authors have adequately addressed in their rebuttal. They have committed to releasing their code to ensure **reproducibility**, added a new experiment on the biharmonic equation to **broaden the examples**, provided sufficient justification for their **choice of baselines**, and presented promising preliminary results for **initializing numerical solvers**.

After weighing the paper's clear strengths against its minor weaknesses, I recommend **acceptance**. This principled approach, backed by strong results, will be a valuable contribution to the community and will likely inspire further work. I strongly encourage the authors to incorporate the new experiments and valuable discussions from their rebuttal into the final version of the paper.
Decision: Accept (poster)