Solving High Frequency and Multi-Scale PDEs with Gaussian Processes

C4: "What would we find if we included a non-machine-learning-based solver in the experiments and compared the runtime to the accuracy of all these methods? I would like to see this comparison. I am not saying that the proposed method must outperform traditional solvers, which have been studied and optimised for a long time, but that the context is critical."

R4: Great suggestion! We actually have included the result of a Spectral method based on Fourier bases plus least mean square (denoted by "Spectral method") in our tables. We followed your advice to test with three other non-ML solvers. The first one is finite difference (discretion resolution = 400, no less than the number of collocation points used in our method), implemented by PyPDE library. The second one is Chebfun that uses Chebshev polynomials. The third one is spectral Galerkin method (with legendre and chebyshev polynomials, 48 quadrature nodes), implemented by shenfun library. We list the results on solving the 1D Poisson equations in the paper. We highligh our method and the best of the other methods.

We can see that finite difference though gives reasonable results, its accuracy is consistently inferior to our method. The solution accuracy of Chebfun is consistently high. Spectral Galerkin, however, seems to fail in all the cases. The three methods run much faster than our method and PINN. This is reasonable, because the ML solvers use optimization to approximate the solution rather than interpolation. We will supplement the results for all the test PDEs in our paper.

C5: "This could be personal taste, but I think the paper might benefit from a slightly more rigorous notation surrounding Equation (14), for instance, by clarifying which joint probability density Equation (14) is the logarithm of. I am imagining something like mentioning "

R5: This is a great suggestion. We will add the clarification of the joint probability (as suggested) in Eq4 to make the whole probabilistic modeling structure (learning objective) more clear.

Thanks for your insightful and constructive comments. Here are our responses. C: comments; R: responses

C1: "The manuscript misses a range of closely related work about solving partial differential equations with Gaussian processes"

R1: Thanks much for providing these excellent references. We do agree these works are related and important. We will cite and discuss them in our paper to appreciate their contributions.

C2: "I am surprised the submission never mentions traditional approaches to solving multi-scale problems, such as those based on numerical homogenisation (a starting point could be [4]). I know that the paper targets the ICLR community and that large parts of this community may be more interested in machine-learning-based PDE solvers than traditional PDE solvers. Still, existing approaches should be credited by (at least) mentioning their existence. (This relates to point 4. below)"

R2: We appreciate the reviewer suggesting and providing references for traditional approaches. We do agree. We will add references and credit those methods (including [4]) in our paper.

C3: "I am also surprised that the experiments only discuss the precision of the approximations rather than the work required to achieve this precision. The experiments in Appendix C suggest that the per-iteration runtime of the proposed algorithm is comparable to that of a physics-informed neural network. However, Figure 9 indicates that tens of thousands of iterations are needed. What is each algorithm's overall runtime (training time) to achieve the precisions in Tables 1 and 2?"

R3: Great observation and questions. Actually, the number of iterations needed to reach a decent accuracy depends on the problem. That is why we report the average per-iteration time. In most cases, our method needs <100K iterations to achieve the accuracy reported in the tables, e.g., 10K for $u_1$ and 20K for $u_3$ for 1D Poisson. In the most challenging cases, including $u_5$ for 1D Poisson and $u_7$ for 2D Poisson, it needs more iterations (>=200K). This is also comparable to the number of iteration required in PINN based methods, e.g., in [1], the PINN is trained with 10K ADAM iterations and 50K L-BFGS iterations. Thanks for the great suggestion. We will supplement the detailed total running time for each test case in our paper.

[1] Shibo Li, Michael Penwarden, Yiming Xu, Conor Tillinghast, Akil Narayan, Robert Kirby, and Shandian Zhe, “Meta Learning of Interface Conditions for Multi-Domain Physics-Informed Neural Networks”, The 40th International Conference on Machine Learning (ICML), 2023.

Dear Reviewer J8Ew,

We have followed your suggestions to update the new version of our manuscript, where new results, like running non-ML solvers (including chebfun), total running time, and polished statements are added and highlighted in blue.

Since the End of author/reviewer discussions is just in one day, may we know if our response addresses your main concerns? We are happy to engage in more discussion and paper improvements.

Thank you again for reviewing our paper!

Authors

C4: "The authors offer very little conclusion and the main limitation section from the Appendix should appear in the main text. One of the key limitation of the method is that it is limited to simple geometry, where one can use a Kronecker product structure. However, this severely limits the applicability of the method and questions its advantages over spectral methods."

R4: Due to the space limit, we leave the limitation section in Appendix. We will re-organize the paper to merge the limitation and conclusion section, and put it in the main paper.

We believe our method is NOT limited to simple geometry. Even if the domain is not a grid, we can create a grid that covers the domain, and use the GP conditional mean (see Eq4 and Eq13) to predict the boundary values and fit the boundary conditions via the GP likelihood $\text{log} \mathcal{N}(\mathbf{g}|\mathbf{u}_b, \tau_1^{-1}\mathbf{I})$ (see Eq14). We then maximize the log joint probability to solve the problem. Note that the computation efficiency (via Kronecker product plus tensor algebra) remains the same.

C5: "Is there any concrete advantage of using the Student t distribution over the Gaussian distribution?..."

R5: Great question. As explained in Section 3 (see the text under Eq5), the student $t$ distribution has a fat-tailed density, and hence can more robustly capture "the (potentially many) minor frequencies". By contrast, the Gaussian distribution density has a very thin tail, and can be much more sensitive to minor frequencies, especially when they are not ignorable.

However, we agree that our most test cases comprise only one or just a few principle frequencies, and hence the results do not exhibit a strong contrast of the two choices. But two examples worth attention. One is the 1D Poisson with solution $u_4 = x\cdot sin(200 x)$, and the other is 1D Poisson with solution $u_5 = sin(500x) - 2(x-0.5)^2$. In both cases, the solution contains many minor frequencies (due to non-trigonometric terms). We can see that our method with student $t$ mixture (GP-HM-StM) achieves much smaller error (one magnitude smaller) than with Gaussian mixture (GP-HM-GM).

We will add more such test cases to highlight their difference.

C6: "The experiments in Section 6 should provide the source terms (at least in Appendix) to allow for reproducibility."

R6: Thanks for the suggestion. As explained in Section 6, specifically, the text under Eq17, Eq18 and Eq19, we follow existing literature to first craft the solution, e.g., $u_1 = \sin(100x)$ and $u_2=\sin(x) + 0.1\sin(20x) + 0.05\cos(100x)$, and then apply the differentiation operator to obtain the source term $f$. We will supplement differential results (e.g., $f$) in the paper.

C7: "Section 6: Timings should be reported in main text as it's one of the main advantage of the proposed method."

R7: Thanks for the suggestion. We will add the time table to the main paper.

Thanks for the many valuable questions and comments. Here are our clarification and responses. C: comments; R: responses

C1: "One of the weaknesses of the paper is that the difference between prior GP-based works in the literature is unclear. The authors should clearly state the differences between their work and the existing ones, in particular the papers by Chen et al (JCP, 2021) and Raissi et al (JCP, 2017), cited in this work. I believe that the two differences are: (1) The use of student-t distribution and (2) using a Kronecker product grid to speed-up computations. In the current version, the title and abstract of the paper might be misleading to readers in the sense that this paper is not the first one to consider solving a PDE with GP."

R1: In fact, we did explicitly emphasize the key difference with Chen et al (JCP, 2021). Please see the bottom of page 4 and top of page 5. In Section Related Work, we have also discussed the challenge of the SE/Matern kernel used by Chen et al (JCP, 2021) in solving the high-frequency/multi-scale PDEs (4th line, 2nd paragraph), which motivates our design of mixture kernel (a 3rd critical difference). We agree your summary about the differences as well. We believe the key difference of Raissi et al (JCP, 2017) has also been covered by our paper --- see Section Releated Work, the 1st sentence of 2nd paragraph --- it only solves linear PDEs with noisy measurement of source terms. Our work (as well as Chen et al (JCP, 2021)) is more general in that it applies to both linear and nonlinear PDEs. We will highlight the difference in a stronger manner.

Second, we never claim/imply that our paper is "the first one to consider solving a PDE with GP". In fact, In Introduction Section (2nd paragraph, Page 2), we explicitly claim that "we resort to an alternative arising ML solver framework, Gaussian processes (GP)" with (Chen et al., 2021) cited. We will emphasize more on this point.

C2: "I am very surprised by the Spectral method experiment performed by the authors. A quick experiment with Chebfun (https://www.chebfun.org/) allows me to solve the 1D Poisson equation to 10 digits accuracy in a fraction of a second."

R2: Our spectral method was implemented based on Fourier bases to be consistent with kernel construction in our method for a fair comparison. The performance is indeed inferior (as shown in the tables). This is to show that introducing a mixture model in the frequency domain (rather than the temporal domain), with which to construct the GP, is much more effective.

Second, we do agree there can be more effective traditional solvers, such as Chebfun (but not all the numerical methods work better for this problem; please see R4 to Reviewer J8Ew below). In many forward problems, machine-learning (ML) based PDE solvers, like GP or PINN, have not been able to perform equally well (in accuracy/efficiency) as traditional solvers. However, as explained in the 1st paragraph of Introduction, ML-based solvers "are simple to implement, and can solve inverse problems efficiently and conveniently". In addition, GP models can easily, robustly handle noisy measurements and/or source terms, offering a principled and convenient framework to quantify and reason under uncertainty. We therefore view developing ML-based solvers as a promising direction, like the concurrent blooming efforts in this direction, even in several aspects the ML-based solvers are now still behind traditional numerical solvers that have been developed for decades.

C3: "The authors state at the end of Section 5 that their work is "the first to [...] use the Kronecker product structure to efficiently solve PDEs". I believe this sentence is erroneous. It is completely natural when using spectral methods to exploit the Kronecker product structure of the domain (see e.g. Fortunato & Townsend, 2019)."

R3: Thanks for the comment. We do agree that our statement is not precise. We mean to say GP plus Kronecker product structure (rather than Kronekcer product structure only). We will revise our statement to be accurate and clear. Please see R2 to Reviewer GK3t above for the detailed clarification.

"Hence, the authors first reported a high error using spectral methods for one experiment but updated it to 10^{-10} error in the response to Reviewer J8Ew."

We believe this has already been clarified in our paper and response R2.

First, spectral method contains a large number of instances. Chebfun is just one instance, which adopts Chebshev polynomials as bases. One shouldn't view Chebfun is the only choice of spectral methods. You can have many other choices. In fact, different bases can end up with different performance. In our paper, we used the Fourier bases (the most straightforward choice), and it does show inferior performance. As explained in R2, we believe putting such a baseline is important to highlight the advantage of designing a mixture model in the frequency space. As reported in R4 to J8Ew, we have also tried spectral Galerkin methods, which completely failed in our ``simple'' cases.

Second, we would love to come out more test cases and more challenging cases. But even for the current test cases, finite difference --- which is another commonly used and fundamental numerical approaches --- have already been surpassed by our methods. Though Chebfun works very well in 1D cases, it cannot solve 2D stationary problem, like 2D Poisson and 2D Allen-Cahn in our paper; please see the updated draft. Similarly, can you write a Chebfun solver to solve a 2D Darcy-flow? The usage of Chebfun is limited, even in the set of cases tested in our paper.

Third, circling back to the points we have made for ML based solvers. One biggest advantage of ML based solvers we believe is that they are really convenient to implement and use --- they are much more general in that you don't need to change the modeling and computational framework to handle different PDEs. After 70+ years' development, traditional methods contain many useful tools/tricks. If you pick up the right tool/trick, it indeed works well; but how to pick up and combine various tools for different PDEs to make them work is always a headache. Of course, ML solvers. which start to get attention for only a few years, are far from mature, in both practice and theory. We therefore view such immaturity as a sign of new promising direction, like many researchers do, rather an excuse to suppress the direction.

"Concerning the complexity of the geometry, the authors reply that one can always embed the domain into a square with a uniform grid. I agree with this comment but the number of points required to approximate the domain geometry would likely grow extremely fast, which would result in a much less efficient representation than using a meshing algorithm combined with a finite element solver. "

We would like to point out that, even if the statement is true (though rigorous proof is unkown), including a large number of points on the boundary will not be an issue in that various stochastic mini-batch optimization can efficiently handle those points, just like training neural networks (or PINNs) with a huge number of data points. There is not any memory or computational bottleneck. Beside, it is very easy to leverage the modern distributed ML platform, e.g., [1], to further speed-up, which is much cheaper than HPC/super computers.

[1] Li, M., Andersen, D. G., Park, J. W., Smola, A. J., Ahmed, A., Josifovski, V., ... & Su, B. Y. (2014). Scaling distributed machine learning with the parameter server. In 11th USENIX Symposium on operating systems design and implementation (OSDI 14) (pp. 583-598).

Thanks for your valuable comments and suggestions. We will polish our paper and improve the presentation to highlight the contributions in a more concise and easy-to-understand way, especially for readers who are not experts in Gaussian process, kernel methods and PDE solving. We will also follow your suggestions to move figure 2 and 3 up to better articulate the problem we aim to solve, and add additional analysis and examples to illustrate the learning behaviors of random Fourier PINNs.

We thank the reviewer for their careful and constructive comments.

C1: "... Härkönen et al. (2023) construct a specific kernel for linear PDEs, which as they write explicitly at the end of Section 4.1 recovers the spectral mixture kernel as a special case"

R1: Thanks for the excellent reference. We will modify our claim. Actually, our work is done roughly at the same time; this is not the first submission (and we did not release our work on ArXiv). Anyway, we will remove the statement, cite and discuss this nice work.

C2: "...I disagree with this being the first instance of Kronecker product structure being used to efficiently solved PDEs..."

R2: Thanks for the valuable comments. Our statement is indeed not precise, and we will remove the claim. From the context, we mean to say, we are using GP + Kronecker product structure --- taking advantage of the corresponding computational advantage --- to solve PDEs (rather the Kronecker product itself). The computational advantage is fulfilled by combining the properties of Kronecker product of kernel matrices and tensor algebra, which we believe is rarely used in the domain of PDE solving. In fact, we have acknowledged tensor-product structures are widely used in PDE solving; see the text under Eq11.

We never claimed that our work invents/discovers the computation advantage of GP+Kronecker product. As emphasized in related work, this advantage has been known in the GP community; our contribution is to enable this computational advantage in PDE solving.

We do agree (Wang et al. 2021) have also proposed the product kernel, but the paper is mainly about theoretical analysis (e.g., the sample path property) and does not fulfill efficient computation. Thanks this nice paper --- it provides additional strong justification of using a product kernel; we will cite and discuss about it. Thus we believe, the contributions of our work are still orthogonal to existing work in a large degree. We will make a more detailed, careful discussion and clarification.

C3: "...the vanilla GP baseline with a Matern kernel. Is this using a product kernel?..."

R3: Great question. We did use product kernels for vanilla GP baselines. We use the same computational approach (Kronecker product + tensor algebra) to fulfill efficient training and prediction. We will clarify it in the paper.

C4: More reference on GP based PDE solvers

R4: We appreciate the reviewer for providing these great references. We will cite and discuss them in our work.

C5: "What about the marginal uncertainty output by the Gaussian process?..."

R5: Great question! We show the result on 1D Poisson (with solution $u_1$ and $u_3$) in Fig. 1 and Fig. 2. We compared our method with student-t mixture (GP-HM-StM) and GP with Matern kernel (GP-matern), in terms of log variance (left), log absolute error (middle), and prediction (right). We can see that the marginal uncertainty indeed reflects the error. In both figures, GP-Matern exhibits much larger errors in almost all the locations, and its corresponding variance is much larger than GP-HM-StM as well. In addition, both GP-Matern and GP-HM-StM show smaller variance on the boundary, and meanwhile their boundary errors are very small. It reflects both methods fit the boundary well. We will supplement the result and discussion in our paper.

C6: "What's the impact of the number of components of the mixture kernel on the approximated solution? ..."

R6: Great question! We have followed your suggestion to give the result on 1D Poisson (with solution $u_1$ and $u_3$). We varied the component number from {3, 10, 30, 50, 80}. See Here. As we can see, using two few or too many components can both increase optimization challenges (too few can make it harder to move to the right frequency). The best choice is in between. We will supplement the results and discussion in our paper.

C7: "Do the kernel matrices you define have Toeplitz structure? Stationary kernels on a 1D grid would have. This could further accelerate the required computations (see Section 3.2. of Wilson et al. 2015)."

R7: Great question. Our method does not require the grid must be evenly-space in each dimension. In fact, we can randomly sample locations in each dimension and then generate the grid (see text above Eq12). Hence the kernel matrices do not have the Toeplitz structure in general. However, if we specify a regular (evenly-spaced) grid, the Toeplitz structure will show up and we can certainly leverage it to further accelerate the computation. Thanks for suggesting this great idea, and we will add it into our discussion.