Toward Efficient Kernel-Based Solvers for Nonlinear PDEs

Thanks for your valuable and constructive comments.

The method assumes a structured grid, which limits its generalization to unstructured meshes and other complex discretization routines.

R1: We appreciate your insightful comments. Indeed, we agree that our method relies on a structured grid and does not directly apply to more complex discretizations. A simple yet effective workaround — which we have validated (see Section 6 and Appendix Section C.3) — is to employ a virtual grid that covers the irregular domain.

In future work, we plan to explore two directions to better generalize our method. First, we will investigate domain decomposition and hierarchical grid construction to adaptively adjust the grid resolution across local regions while preserving computational efficiency. Second, we aim to learn a mapping from unstructured mesh points to a latent grid, where our efficient kernel solvers can be applied. This approach is inspired by recent work such as GeoFNO[1] in the neural operator literature.

We will include this discussion in the paper.

[1] Li, Z., Huang, D. Z., Liu, B., & Anandkumar, A. (2023). Fourier neural operator with learned deformations for pdes on general geometries. Journal of Machine Learning Research, 24(388), 1-26.

The hyperparameter selection process (e.g., kernel length scales) is not discussed in detail.

R2: Thank you for your great suggestion. Given the wide range of hyperparameters --- including $\alpha$ and $\beta$ (see our response R3 to Reviewer fQAE), as well as kernel length-scales and nugget values (see Lines 284–297 right column) --- performing a full grid search would be prohibitively expensive. Therefore, we adopt a hybrid strategy. We begin with a random search to identify a promising set of hyperparameters. Then, we perform a grid search over $\alpha$ and $\beta$, keeping the other parameters fixed. Once $\alpha$ and $\beta$ are selected, we fix them and conduct a grid search over the remaining hyperparameters, including the nugget and length-scales. We will include a detailed description of this procedure in the paper.

I'm also curious how the runtime of the Kronecker based method compares with the naive full matrix method.

R2: Great question. Below, we provide the runtime of our model using the naive full matrix computation (i.e., without exploiting the Kronecker product structure). As shown, the per-iteration runtime with naive matrix operations is consistently around 100× slower compared to our method, which leverages Kronecker product properties. Furthermore, when the number of collocation points increases to 22,500 for the Burgers’ equation and 43,200 for the Allen–Cahn equation, the naive approach exceeds available memory, resulting in out-of-memory (OOM) errors — rendering it infeasible to run. We will include these results in our paper.

We thank the reviewer for the valuable and constructive comments.

C1: Can you provide more details on the virtual grid on an irregular domain, which is important for the application of your method?

R1: Great question. Our virtual grid is chose as the smallest rectangular grid that fully covers the irregular domain. Specifically, for the nonlinear elliptic PDE (see Appendix C.3 and Figure 4), the domain is a circle inscribed within $[0, 1]\times [0, 1]$, and we hence set the virtual grid on $[0, 1]\times [0, 1]$ with a grid size of $50 \times 50$. For the Allen-Cahn equation (see Appendix Figure 4), the domain is a triangle with vertices at (0, 0), (1,0) and (0.5, 1). To fully cover this domain, we again place the virtual grid on $[0, 1] \times [0, 1]$, using a grid size of $50 \times 50$. We will include these details in our paper.

C2: Can you provide the running time comparison with the finite difference method?

R2: Great suggestion. Below, we provide the time comparison results in seconds (for the Allen-chan equation with $a=15$).

As we can see, the finite difference (FD) is computationally much more costly per iteration compared to our method (SKS). This might be due to the expensive cost of computing the inverse Jacobian during the root finding procedure. However, FD typically converges within just a few dozen iterations — substantially faster than the stochastic optimization used in both SKS and PINN — leading to a much lower total runtime. We will include this comparison in the paper.

C3: In future work, can the proposed method be extended to the neural operator setting for further acceleration?

R3: We appreciate the reviewer for bringing up this excellent idea — we completely agree. In fact, we have already begun exploring this direction in our ongoing work. One line of effort involves replacing the trunk network in the Deep Operator Network (DeepONet) with Gaussian process bases. This enables us to adopt a similar computational strategy to accelerate both training and prediction, while also providing a natural framework for uncertainty quantification. We are also extending our approach to deep kernel-based operator learning, where our method can further accelerate function transformations in the latent channel space. We look forward to continuing our exploration in this promising direction.

Thanks for your valuable and insightful comments.

The introduction lacks references. Only the direct relevant work is cited (2 papers)

R1: Thank you for the helpful suggestion. We will include additional references in the introduction to provide broader context and better situate our work within the existing literature from the outset.

no standard deviations

R2: Great comment. In our experiments, we observed that despite using stochastic optimization (ADAM), our method (SKS) consistently converge to the same solution. So is PINN. We ran multiple trials on each PDE and found that the standard deviation of the error is negligible. For example, below we show the standard deviation of the $L^2$ error for our method when solving Burgers' equation ($\nu = 0.01$) and the nonlinear elliptic equation. Given the extremely low variance, we chose to omit standard deviation values from the main results (Note that the finite difference method is deterministic and does not exhibit variance). We will clarify this point in our paper.

Burgers($\nu = 0.01$):

Nonlinear elliptic:

I think that there is a problem with the formulation of (7). As it is, $\epsilon$ has no impact on the problem thus I hardly see how it can be equivalent to (6)... A typical way to do this is to square this quantity again.

R2: Thank you for your insightful question. Here is our clarification.

First, $\epsilon$ plays an important theoretical role (though in practice, it can often be set to zero). Our convergence proof and rate estimate (see Lemma 4.2 and Appendix Section A) are established by systematically varying $\epsilon$. Specifically, we set $\epsilon = C_0 h^{2\tau}$ and vary $\epsilon$ by adjusting the fill-in distance $h$. Each choice of $\epsilon$ defines a distinct instance of problem (6), with its own solution. We prove that as $\epsilon \to 0$ (i.e., $h \to 0$), the solution of (6) converges to the ground-truth solution of the PDE (see Appendix Section A for details).

Second, to show that solving (7) (with appropriately chosen $\alpha$ and $\beta$) is equivalent to solving (6), we use the Lagrangian formulation of (6). It is a standard result that constrained optimization problems can be reformulated as mini-max problems over the Lagrangian. The Lagrangian is a linear combination of the objective and the constraints — it does not involve squaring the constraints. From this perspective, it is straightforward to show that when $\alpha$ and $\beta$ are selected as the optimal dual variables in the mini-max problem, solving (7) yields the optimal $u$, and thus is equivalent to solving (6). This equivalence holds for any $\epsilon$, not just in the limit as $\epsilon \to 0$. The full proof is provided in Appendix Section B.

Finally, we agree that an alternative approach --- minimizing the objective plus the squared constraints (i.e., a penalty method) --- is also a good and viable idea. However, this method is primarily justified for equality constraints (i.e., $\epsilon = 0$) (While it is possible to design penalty terms for inequality constraints, they typically introduce non-differentiable terms.) Moreover, equivalence to the original constrained problem is only guaranteed as as the penalty weights tend to infinity. Therefore, we believe our formulation in (7) offers stronger theoretical guarantees and greater practical flexibility when approximating solutions of (6). We will include a more detailed discussion of this alternative approach in our paper.

Could you include a discussion about the tuning of two regularization parameters (7)?

R3: Thank you for the great suggestion — we completely agree. In our experiments, we selected $\alpha$ and $\beta$ in Equation (7) from a wide range: $[10^{-2}, 10^{-1}, 1, 10, 10^2, 10^3, 10^4,10^5, 10^6, 10^7, 10^8, 10^{10}, 10^{12}, 10^{14}, 10^{15}, 10^{20}]$, jointly with other hyperparameters, including the kernel length scales and nugget terms (see their ranges in Lines 284-297 right). To efficiently tune hyperparameters, we first performed a random search to identify a promising group of hyperparameters. We then fixed all other hyperparameters and conducted a grid search over $\alpha$ and $\beta$. Finally, we fixed $\alpha$ and $\beta$ and performed a grid search over the remaining hyperparameters. We will include this discussion in the paper.