We appreciate the comment on the readability and simulation results. In response to your inquiries, we would like to provide detailed explanations below:

1. General Contributions: (Please refer to General Response for more details)

Thanks for recognizing the technical contributions of our paper. As for our main contribution, our goal is to introduce the idea of combining the Neural Field method with the Monte Carlo method using Control Variates. This contribution is very general as it can be applied to any spatial integration setting where the domain is structured. Unlike previous neural control variates methods where the networks are restricted to specific types of architectures, our methods impose minimal restrictions on the choice of network architectures. This allows users to choose more powerful and expressive network architectures, potentially leading to better numerical integration tools in the future. We believe that our contributions add knowledge to the field of numerical integration and can be beneficial to users with various applications in mind.

2. Wall-Time-Comparison:

Since our current implementation is based on Jax, which is highly under-optimized for both WoS and Neural networks, wall-time comparison might not be very reliable. It also depends on the hardware where we run the experiments in, as well as what kind of PDE and neural networks we used for the experiment. As a result, we prioritized to present equal sample comparison, as done by many papers (Sawhney et al. 2022).

Per the reviewer’s request, we’ve added the of our wall-time breakdown of our methods in comparison with the baselines, alongside our discussion of wall-time performance in the limitation sessions (Table 1 and L262-L265). Right now our proposed method takes more time compared to WoS to reach a certain MSE threshold in solving the 2D Poisson equation. We believe that the current profiling can be greatly improved if we investigate different ways to make neural network inference faster, such as implementing customized kernels. We also hypothesize that our method’s advantage will be more obvious when applying to solve more complicated PDEs, where WoS baseline will require more time per evaluation of integrands $f$ or per sampling from the domain.

3. Optimizing a network $g$ or just $c$ with a fixed, approximated $g$:

If we only optimize $c$ with a fixed $g$, then the error between the shape of $f$ and the shape of g cannot be fixed. On the other hand, optimizing the network $g$ allows us to match f not only in scale but also in shape. As a result, optimizing a network $g$ can be more powerful to reach lower variance.

Let’s examine the following example. Let $f(x) = x^2$ and $g(x) = x$. No matter how we optimize the $c$ we cannot reduce the variance of a Monte Carlo estimator for $\int_0^1 f(x) - cg(x) dx$ to zero. With our methods, we can define $G(x) = a_0 + a_1 x + a_2 x^2 + a_3 x^3$, where $a_0, \dots, a_3$ are optimizable parameters. We can see that $g(x) = G’(x) = a_1 + 2 a_2x + 3a_3 x^2$. One way to minimize the variance is to set $g(x) = f(x)$, namely $a_1 = 0, a_2=0, a_3 = 1/3$. This way our method can achieve zero variance, which is not possible if we only optimize with a fixed $g$ of choice: $g(x) = x$. In other words, our methods allows the data to pick the appropriate choice of $g$, which can be more efficient comparing to prior methods where $g$ is chosen via some (potentially wrong) heuristics.

4. Plotting issues:

The actual result has a few slight bumps accumulating for 10k steps, you can zoom in to see the small bumps.

Thank you for your thoughtful review! We would like to address some of your concerns as follows:

1. Convergence Rate of our method:

The convergence rate of our proposed method is also O(1/N). This is shown in Figure 2 Convergence Plot and Section 5.3 (L233-L235)

2. Application Scenarios:

We have discussed this in the second paragraph of the introduction section(L21- L27). Applications include solving PDEs and physics-based rendering, etc. These applications greatly impact various industries including movie, gaming, and mechanical engineering. Specifically, in this paper, we demonstrate our method’s ability to produce more accurate and efficient PDE solvers for Poisson and Laplace equations, which are among the most fundamental PDEs.

3. Limited Numerical Experiments: (Please refer to the general response for more details)

Per the reviewer's requests, we have provided more numerical experiment results. The additional results are consistent with those we provided within the paper and these results also support our claims that the proposed neural control variates estimator is unbiased and also low-variance. The additional results further suggest that our method can generalize to different 2D/3D domains and boundary conditions.

4. Lack of theoretical guarantees of variance:

We would like to clarify the misunderstandings of the variance guarantee. Our method inherits the merits of Monte Carlo methods - unbiased and with predictable variance. As shown in the paper shows that our method’s variance can be proved to decrease in O(1/N) (Equation 15 and L233-L235)

Furthermore, our network is directly optimized to reduce the variance, which we’ve shown in Equation 16. These advantages come with very few restrictions on the network architectures and parameters except for letting the network be n-th time differentiable, which we believe greatly expands the choices of functions for control variates.

Thank you for pointing out all the typos and unclear notations. We have fixed them in the main paper and marked the change with a cyan color. We hope our revision can improve the mathematical rigorousness and clarity of our method section. We would also like to address the following questions:

1. Computational complexity when d is large:

Thank you for asking this question. In our original paper, we discussed the computational complexity for different d: “Computing the integration requires evaluating the network approximating the anti-derivative $2^d$ times, with $d$ being the dimension of the space we are integrating in.” ( L271-L274)

(Please refer to the General Response for details.)

We've discussed how our method scaled with $d$ naively in the limitation in the paper. Empirically, our method could also be applied to higher dimensional problems with efforts currently: Our PDEs’ path integration domain is infinite dimensions. We demonstrate the applicability of our method by breaking down the infinite-dimensional integration domain into smaller dimensions. In other words, even though our techniques could not naively be applied to higher dimensional problems directly, as long as the problem could be broken down into smaller dimensions, our method would work out. In the meantime, we believe extending our method to be naively applicable in higher dimensions is an interesting future work problem. We’ve added this to the discussion in our current paper. (L274)

2. Axiomatic Unbiased:

Thank you for pointing out this language issue. Our method is unbiased by construction. The experiments' results further confirm this as our MSE decays at the rate of O(1/N) where $N$ is the number of samples. We’ve corrected the language issue in the paper in L55-L56.

We greatly appreciate the insightful comments and suggestions provided by the reviewer. In response to your questions, we provide detailed explanations below:

1. Explain the importance of solving PDEs:

Thank you for this suggestion! We originally motivated our applications in the introduction section (L21-L27). Per reviewer’s request, we also include a motivation in the result section to improve the flow (L189-L192 in the paper).

2. Scalability: (Please refer to the General Response for more details.)

Thank you for providing the relevant papers. As noted in the general response, our experiment is a demonstration that our method can be applied to solve integration problems in higher dimensions. We believe that similar principle might be applicable to other applications where high-dimensional integrations can be broken down into a set of lower-dimensional integrations. We believe that exploring how to apply the principle proposed by our method to different applications will be great future work research. We’ve also included discussion in our revision of the paper. (L271-274)

3. Using parametric function instead of NN:

The parametric function can be used. Some specific parametric functions (such as a weighted sum of polynomial kernels) have been studied by prior works (Salaün et al. 2022) (L349- L350). In the meantime, our methods open a wider variety of network architectural choices for control variates. We achieve this by proposing a more general way to construct the control variates pairs (i.e. the G and g such that G = \int g), placing minimal restrictions on the network architecture. This allows us to use not only differentiable parametric functions but also other powerful universal approximators such as deep neural networks with nonlinearities.

4. Bayesian Inference:

We agree that extending our methods to Bayesian inference can be a very interesting research direction. One of the key insights of this paper is that when we need to construct a pair of functions for control variates, namely $G$ and $g$ such that $G = \int g$, we can first define a function to approximate $G$ and take its derivative fo find $g$. We believe that key insights like this can be potentially extended to other problems that might benefit from control variates, including Bayesian inference.

5. Result aggregating different $n$:

In the 3D Laplace experiment setting (Fig.3 in the paper and provided above), we have provided an example with different sample sizes for our method. As the sample size became larger, our method results were closer to ground truth. Our error is decreasing in the rate of O(1/N)