Toward Efficient Kernel-Based Solvers for Nonlinear PDEs

C2: Beyond this issue, I also do not understand why the choice of tensorial kernel makes any sense except for computational convenience

R2: As we have emphasized throughout the paper, our motivation for introducing the product kernel is rooted in reducing computational costs and enabling the use of a large number of collocation points --- strictly from a computational standpoint. We do not claim any additional implications beyond computational convenience. It is surprising to see the reviewer’s strong focus on points that are tangential to our objectives, indicating a misunderstanding of our motivation. Nevertheless, the product kernel corresponds to a tensor product structure in the latent feature space, consistent with the tensor product approach commonly employed in numerical methods [1].

[1] ARNOLD, D. N., BOFFI, D., and BONIZZONI, F. (2012). Tensor product finite element differential forms and their approximation properties. arXiv preprint arXiv:1212.6559.

C3: Why solving the minimization problem in (6) is better than solving (2a) needs some clarifications? Can you provide numerical evidence in terms of accuracy and efficiency? The equation numbers here correspond to the equations labelled in the comments above.

R3: We did try with direct optimization of the coefficients $\mathbf{\alpha}$ instead of using the approach in Eq. (6). However, we observed that this led to a solution error increased by 1 -- 2 orders of magnitude. Our empirical findings suggest that each coefficient $\mathbf{\alpha}$ globally influences the solution approximation, to which both the boundary condition fit and the PDE residuals are highly sensitive. As a result, small adjustments of the coefficients can cause significant perturbations, leading to an imbalance between these components. This makes the optimization process much more challenging than directly estimating the local solution values at the collocation points. Below, we present results from solving the Burgers' equation with $\nu=0.02$ to illustrate this effect.

C4: Can you clarify what DAKS is and what its relationship to their current work?

R4: As highlighted in Lines 353-355, DAKS (Chen et al., 2021a) is the prior work that motivates our approach. DAKS is based on the nested optimization framework introduced in Eq. (2) of Section 2. It incorporates a set of derivative values — more generally, the linear operators within the PDE, evaluated at the collocation points, denoted as $z^j_m$ in (2). DAKS then uses kernel interpolation to construct the solution approximation, as shown in Eq. (3), where the Gram matrix consists of sub-blocks computed by applying linear operators to the kernels (see Eqs. (4) and (5)). All relevant details have been provided in Section 2

We thank the reviewer for their time. However, the reviewer seems to completely misunderstand the motivation and contribution of our work. The reviewer is even unclear what the major baseline method DAKS (Chen et al., 2021a) is, which directly motivates our work, and has been explained and highlighted in numerous places in our paper, e.g., Line 34-45 and whole Section 2.

C1: I just cannot see any new idea with this formulation. BTW, if the PDE is nonlinear, I believe RBF will result in a minimization problem as well

R1: It appears the reviewer has misunderstood our motivation and contribution, resulting in a superficial assessment. The objective of our work is not to propose a new framework to replace kernel ridge regression. As highlighted throughout the paper, our aim is to develop an efficient kernel-based approach for PDE solving that (1) is computationally efficient and convenient to implement, and (2) provides convergence guarantees along with rate analysis. Evaluating a method without understanding its objectives is unproductive. Following the reviewer’s logic, the baseline method, DAKS (Chen et al., 2021a), could also be dismissed as merely a variation of kernel ridge regression, lacking "new ideas". By this reasoning, numerous neural network papers could similarly be deemed unoriginal just because they employ neural network formulations.

While the reviewer proposed equations relevant to solving linear PDEs, our work focuses on developing a nonlinear PDE solver (see our test benchmarks). The context and results are entirely different, as nonlinear cases are much more challenging. For instance, in the nonlinear setting, it is generally not possible to achieve the optimal approximation in the interpolation form specified in Eq. (6) in the review comments.

The real challenge is: for nonlinear PDEs, if one applies kernel ridge regression directly, justifying optimality and convergence becomes very difficult. The reviewer might be aware of the the theoretical foundation of kernel ridge regression. When we solve an optimal recovery problem framed as minimizing $\left\lVert u\right\rVert_{RKHS}$ s.t. $u(z_j) = y_j (j=1 \ldots M)$, the representation theory guarantees that the optimal solution takes an interpolation form, $u(x) = \kappa(x, \mathbf{Z}) \mathbf{\alpha}$. Only with this form can we express the RKHS norm of $u$ is $\mathbf{\alpha} \mathbf{K} \mathbf{\alpha}$. Based on this result, the optimal recovery problem can be reformulated as a soft-constraint optimization problem, resulting in the well-known kernel ridge regression: $\min_\mathbf{\alpha} \frac{1}{M}\sum_j u(\mathbf{x}_j - y_j)^2 + \lambda\mathbf{\alpha}^\top \mathbf{K} \mathbf{\alpha}$.

However, if you impose nonlinear constraints, such as $u^2(z_j) =y_j$, the optimal solution no longer takes the simple linear form, and applying kernel ridge regression loses the foundation!

To address this, the prior work DAKS (Chen et al., 2021a) introduces a nested optimization formulation (see Section 2 of our paper for a detailed explanation). This approach ensures that at the inner level, the problem aligns with the standard optimal recovery form, allowing for a claim that, given derivative values $z^j_m$, the optimal solution for $u$ still takes the standard linear/interpolation form. However, this formulation requires augmenting the RKHS with derivative operators over the kernels, which inevitably introduces kernel derivative blocks in the kernel matrix.

Our contribution, in contrast, is to maintain the favorable interpolation form without derivative augmentation (see Eq. (6) in our paper), allowing for a simpler implementation and more efficient computations while reconstructing theoretical guarantees!!. Achieving this required moving beyond the classical optimal recovery framework, as we constrain $u$ to a reduced search space within the RKHS. We have rigorously demonstrated that, even in this reduced space, convergence to the true PDE solution is achievable under mild regularity assumptions commonly employed in PDE convergence analysis. Moreover, our convergence rate remains comparable to that of DAKS. For the detailed proof, see Sections 4, A, and B, where we present a technical and highly nontrivial analysis.

I admitted that I am very confused reading this paper since the presentation does not help improving my understanding either. While the presentation states that the proposed work is trying to improve upon DAKS, but then it is very unclear to me why one wants to regress to a set of features that include C_{ij} corresponds to derivatives of the kernel as in (3)-(5) in the manuscript. Effectively, the PDE constraints readily provide matrix $A = K''+\alpha K$. At the same time, such a question does not seem to be fair since I am not reviewing DAKS. Second, the proposed new algorithm in this paper now seems to abandon the needs of regressing over C_{ij}, which is how it should be done in the first place. But then I am bothered with why would one devise an algorithm to beat previous method that is unclear. This stems me to look at what is actually being proposed if one solves a simple linear PDE, which helps me clarify what this paper seems to be about.

Responding to C2, I am still not convinced why one should do the tensorial kernels. There must be cases when such idea is advantageous and/or not advantageous beyond computational costs. I understand your motivation but the current answer seems to be based on empirical evident of a few examples.

Responding to C3, this comparison is not complete. The fact that minimizing over $\eta$ is better than $\alpha$ requires a lot of explanation. There are many factors that can affect the solutions (including optimization scheme, initial conditions, etc) and whether this phenomenon occurs only on 1 example is unclear.

We thank reviewer for their time, and we would like to reply with the following.

W1: I did not like the trade-off of making the kernel just a constant smaller. ... The product kernel trick was a nice idea, but I am convinced that it can also be applied in Chen et al., 2021a.

In Chen et al., although each block of the full Gram matrix can be converted in kronecker product if product kernel is used. However, the full Gram matrix cannot be decomposed into kronecker product. Thus, we break the big Gram matrix into smaller ones and exploit the kronecker structure for efficient computation.

W2: In the experiments setup they say SKS is run using the ADAM optimizer for 1e6 epochs, but in the table 9 in the appendix, where they report time per iteration and total, the ratio is not 1e6.

For reported runtime, total time is calculated with early stopping criterion. We stopped the optimization, if the performance stopped improving for 1000 iterations. This is the reason why our method total time is not 1e6 * per-iter time.

W3: I also found it very surprising that the time per iteration is so small, especially compared to those of PINNs

We believe our method time per iteration is surprisingly small compared to PINNs is due to the following:

1. Our method is based on linear operations, where we did not calculate gradient explicitly. But PINNs are based on autograd.
2. Optimization in our method is simpler and more convenient for modern libraries. Our method only needs to optimize parameter \eta, where PINNs need to optimize with computation graphs on derivative terms.
3. Our method is non-parametric, the only tunable parameters are the lengthscale parameters. We further improved our computational efficiency by fixing the kernel matrix terms at first iteration.

Q1: Why learn $\eta$, instead of $K^{-1}\eta$?

Great question, we actually tried this approach! But the performance is just subpar, we believe it is due to $K^{-1}\eta$ is hard to optimize. Below, we present results from solving the Burgers' equation with $\nu=0.02$ to illustrate this effect.

Q2: In the numerical experiments, the collocation points should be the same for SKS and DAKS, as I am convinced the higher errors of DAKS arise from random sampling: with random sampling, there may be "holes" in the domain without collocation points, where the error is bigger.

We actually tested with both random sampling and grid sampling and random sampling is constant better than grid sampling, thus, we used 5 different seeds and reported the average.

Below, we present results of burgers $\nu=0.02$ with two difference sampling methods.