Neural Green’s Functions

Dear Reviewer FSUU,

We 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.

Concerns on Effectiveness & Applicability of Neural Green’s Function

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.

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.

Specifically, 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:

$$f(x, y, z) = A (x^3 - 3xy^2) + B(y^3 -3x^2y) + (x^2 - z^2).$$

We 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.

After 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.

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

$$f(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),$$

whose 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.

Runtime Comparison

We 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.

Clarifications

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.

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.

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)}$.

Lastly, we appreciate your recommendation of an inspiring work (PhyMPGN), which seems highly relevant for advancing our approach.

Dear Reviewer KHba,

Thank 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:

Implementation Details for Reproducibility

We 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.

Lack of Examples

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. 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.

We 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

$$f(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).$$

We 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.

The 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.

Limitations

While 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.

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.

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.

Clarification on Baseline Setup

We 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.

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.

Reference

[1] Fast Tetrahedral Meshing in the Wild (ACM ToG 2020): https://dl.acm.org/doi/10.1145/3386569.3392385
[2] cholespy: An easily integrable Cholesky solver on CPU and GPU: https://github.com/rgl-epfl/cholespy

Dear Reviewer fcRy,

We 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.

Extension to Other PDE Problems

We 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.

We 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

$$f(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).$$

We 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.

After 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.

Extension to Shapes with Multiple Components and Stability

We 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.

Discussion of Numerical Stability

Since 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.

Dimensionality of Feature Dimension

We 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.

Sampling Method for Interior Query Points

For 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.

Failure Cases of Neural Green’s Function

We 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.

Limitations

We 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.

Dear Reviewer SXcQ,

Thank 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:

Generalization Across Shape Categories and Boundary Conditions

We 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.

Training Time

Using 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.

Clarification of L111-112

We 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.

Comparisons Against Boundary Element/Integral Methods and Monte-Carlo Methods

Thank 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.

Starting 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.

Although 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.

Impact of Choosing Query Points

Our 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).

Table 1. Statistics of Point Cloud Sizes

Composing Combinations of Domain, Boundary Function, and Source Function

By 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.

Reference

[1] Transolver: A Fast Transformer Solver for PDEs on General Geometries: https://arxiv.org/abs/2402.02366
[2] Universal Physics Transformers: A Framework For Efficiently Scaling Neural Operators: https://arxiv.org/abs/2402.12365