Neural functional a posteriori error estimates

We appreciate the review and thank the reviewer for providing comments on the article.

We understand that in parts the presented work may be hard to follow. This is, however, subjective, since we can see that pN9X, and SLC5 rank presentation as good and excellent, CFkx ranks it as fair and finally, the present reviewer h1aP ranks it as poor. Given that we have a full spectrum of opinions.

The introduction provides minimum information about astra -- only that it is a loss function with certain benefits like more robustness over residuals and variational losses. In my view, there should be a high-level description of (a) what motivates this specific loss function, and (b) what exactly is this loss function. There should also be reasoning behind it at a high level. It is difficult to grasp the concept and also get interested otherwise.

There is a part of the paper that provides high-level details on the Astral loss. We tried to convey this information in 2.1, where the rationale behind using the functional error estimate as a loss function is given (certain favorable properties are listed).

I think the paper can also be more kind to the readers by defining certain terminology before using them, e.g., majorants, posterior error, priori error, etc.

We agree that these terms are not properly defined. We will make sure to provide brief comments when appropriate. More specifically, we will comment on the terms majorant and minorant right after equation (2) where they first appeared. Similarly, we will comment on a priori and a posteriori errors in the introduction where these terms appeared first.

I found section 2.3 to be difficult to parse. The paper states several equations without explaining them in sufficient depth. For example, in equation (9), I did not understand where the definition of astral is from. The paper states "It is possible to derive" but at least in the appendix, these derivations should be shown. Equation 12 is also similar. I could not understand the derivation from equation 9 to equation 12

We agree that the presentation is sketchy here. The main reason why is the limit on the length of the paper. We provided references in the text sufficient for the reproduction of equations but we agree that it is better to make paper as self-sufficient as possible. Shortly we will provide a revised version of the manuscript with two additional appendices: the first one with the derivation of upper bound (9), and the second one with the explanation of how to arrive at (12) from (9).

There should also be sufficient reasoning "why", for example, equation 12 is a good measure of the predicted solution.

We agree that this is an important question and thank the reviewer for bringing this up. Indeed, the manuscript does not explain properly why the energy norm (and its upper bound) is a good measure of accuracy. In some sense this material is standard, but we will review it properly. For details we refer to the answer https://openreview.net/forum?id=z62Xc88jgF&noteId=SQ0cDlqod2, where we address this question in some detail.

We want to thank the reviewer for taking part in the discussion of the manuscript.

We agree that the approach is limited to the PDEs where the functional upper bound can be derived, and this is the main downside. However, for many practically relevant equations, the functional error estimate exists. A non-exhaustive list of these equations is given at the end of 2.3 of the present manuscript.

Here, we would like to mention that error majorants are available for the Maxwell equation and a wide class of variational problems with convex functions.

Besides that, we provide additional examples of how to apply functional error estimate: $D=2$ elliptic equation in the L-shaped, $D=1+1$ convection-diffusion equation. More details are available in https://openreview.net/forum?id=z62Xc88jgF&noteId=Wy4TDErd0T.

As for the code, we decided not to share the GitHub repository because direct sharing is prohibited by rules. On the other hand, anonymization of the repository is tedious. The Google Colab notebooks that we share are anonymous and can be immediately used to reproduce the results since they run in the browser. We will make the code available after the review stage.

Thanks for releasing the code. I am revising my rating up.

"Constructing the the upper bound on the solution error is non-trivial. "

• No. The upper bound is trivially obtained for elliptic problems if one uses functional-type a posteriori error estimation. For Poisson problems, inserting the flux variable $\\|\nabla u - \nabla u_h\\|^2$ and integration by parts you are good to go. The difficult part is to the lower bound locally, in the sense that the local error can be also bounded below by the error estimator (local version). Check Repin's book chapter 2.2.2.2 for example, and also chapter 5 about the local efficiency. For residual-based a posteriori error estimation approaches (dated back to 60s which the authors cited none), the difficulty reverses (local efficiency is trivial but the upper bound needs the discretization's stability and a quasi-interpolant).

"My intuition is that the error certificate $w_U$ can be thought of as an explanation of how the solution emerges"

• $\\|u - u_h\\|_V$ is not computable, because $u$ is unknown (unless in a toy model). If $V$ is Hilbert, then this error is a functional in $V'$, $w_U$ can be viewed as an approximation to the Riesz representation of this functional in a computable norm $L^2$.

We would like to thank the reviewer for the time spent reading the paper and for the feedback.

Typos:

• Indeed, the approximate solution should be $\widetilde{u}$. We will fix that in the revision.
• “deep learning is to” --> “deep learning to”

As for the question. Indeed, the solution may be good be the certificate produced by the surrogate model is not good, i.e., they lead to an upper bound that is much larger than the energy norm. In our experiments, however, we did not observe this behavior when input data is in distribution. Out of distribution, this may happen, but the situation is hard to analyze. In this unfortunate case when certificates are useless, the practitioner has two choices:

1. Recompute the solution with a more reliable method.
2. Optimize the loss with a fixed solution to find better certificates.

So the certificates are useful only in the situation when the predicted upper bound is smaller than the desired tolerance. The situation in machine learning and statistics when the method provides a measure of certainty for the obtained result (e.g., vote counting, document sorting). When the measure is good enough, no intervention is needed. When the measure is bad, the sample is processed with a more reliable method (e.g., by a human being).

Constructing the the upper bound on the solution error is non-trivial and requires expert knowledge, but this provides a way of incorporating domain knowledge into the model.

Indeed, we agree that the construction of the upper bound is not straightforward and requires expert knowledge. We recognize this as a weakness of the approach and provide some discussion in these lines in https://openreview.net/forum?id=z62Xc88jgF&noteId=L8Rln3shGY. This discussion includes the list of equations where the upper bound is known and the additional example of the convection-diffusion equation that we introduce to showcase that not only elliptic equations can be treated with functional a posteriori error estimate (see https://openreview.net/forum?id=z62Xc88jgF&noteId=Wy4TDErd0T).

I have one of the most positive reviews for this paper because I think it is well-written and brings new ideas to the ML community. I admit that it is out of my area of expertise and the other reviewers likely have more knowledge of this area than me, but I think the authors have responded well to the criticisms and comments, so I see no reason to change my positive review.

In particular, the public comments from an expert in this area were critical of the paper's novelty. (https://openreview.net/forum?id=z62Xc88jgF&noteId=1MwG3D88ca) The responses from the authors seem thorough and make sense to me, but I admit I do not have the depth of expertise to be confident in my judgement.

We would like to thank the reviewer for bringing up an important question. We agree that the relation between errors computed in energy and the $L_2$ norm is an important topic not covered in the article. There are two parts to consider: the theoretical one — how the losses are related, and the practical one — what are typical magnitudes of errors in each norm when networks are trained with Astral and residual losses.

The theoretical part is as follows. Suppose the error is $e$, the diffusion coefficient is $a$, and the reaction term is $b^2$. In this case, we can find

$|||e|||^2 = \int_{\Gamma} dx \sum_{ij} \frac{\partial e}{\partial x_{i}}\frac{\partial e}{\partial x_{j}}a_{ij} + ||be||^2 \geq \lambda_{\min}^2 ||e||^2 + ||be^2|| \geq \lambda_{\min}^2 ||e||^2 + \inf_{x} b(x)^2||e^2||,$

where $\lambda_{\min}$ is a minimal eigenvalue of the elliptic problem $-\frac{\partial}{\partial x_{i}} a_{ij} \frac{\partial}{\partial x_{j}}u(x) = \lambda u(x)$ defined on the same domain and with the same boundary conditions as the original elliptic problem. From the expression above, we obtain the bound

||e||\_{2}^{2} \\leq \\frac{1}{\\lambda_{\\min}^2 + \\inf_{x} \\left(b(x)\right)^2} |||e|||.

This bound can be simplified and made more practical (e.g., by removing minimal eigenvalue that may be hard to compute) by various means but the bound above is already sufficient to claim that when error is sufficiently small in the energy norm, this is also the case for the standard $L_2$ norm.

To answer the practical part, we conducted the set of additional experiments described in the reply https://openreview.net/forum?id=z62Xc88jgF&noteId=Wy4TDErd0T. There, we compared relative errors and energy norms for two additional test cases ($D=2$ elliptic equation in L-shaped domain, $D=1+1$ convection-diffusion equation) for Astral and residual losses.

We would like to also add that, indeed, solutions of PDE can be considered in various spaces. Most typically, however, not a space $L_2$ is used but Sobolev spaces to properly define weak solutions (see [E] Part II). Errors in natural Sobolev norms are sometimes used in deep learning literature. For example, a recent paper with FNO extension uses this norm [K]. Other metrics besides relative $L_2$ norms are also used, e.g., $L_1$ norm divided by $L_{\infty}$ [B], relative RMSE of energy spectrum (turbulence) [C], etc

The main reason we use energy norms is because it is a natural norm for elliptic problems. Namely, the energy norm will weigh the error according to the importance of spaces spanned by eigenvalues of the elliptic operator in question (eigenspaces with larger eigenvalues receive more weight). The energy norm is a de facto standard for elliptic PDEs. For example, FEM error estimate is done with error energy norm [SB, for example, 1.5.2 and 1.5.3]; the convergence theory for Conjugate Gradient method (Krylov subspace method used to solve a linear system with symmetric positive-definite matrices including the ones originated from elliptic equations) is naturally formulated with energy norm [S].

Finally, we would like to point out that we do not claim that the solution obtained with Astral loss is better. Our main claim is that Astral loss provides a natural way to construct an error majorant that provides a tight bound on the error. Note that the residual does not have this property. It is easy to see that the residual can be arbitrarily bad for the solution that is arbitrarily close to the exact one. Namely, we can consider a small high-frequency perturbation $\\delta(x) = \\epsilon\\sin(2\\pi L x)$ to the exact solution. This perturbation will introduce an error of the order $\epsilon$ and residual of order $\left(2 \\pi L\right)^2 \epsilon$. Now one can take $L$ large enough to have an arbitrarily large residual.

We hope that this discussion clarifies the relation between $L_2$ and energy norms and explains our choice of metric. We will incorporate this material in the corrected version of the paper that will appear shortly.

[E] — Evans LC. Partial differential equations. American Mathematical Society; 2022 Mar 22.

[K] — https://arxiv.org/abs/2310.00120v1

[B] — https://arxiv.org/abs/2202.03376

[C]— https://arxiv.org/abs/2112.01527

[SB] — Szabó B, Babuška I. Finite Element Analysis: Method, Verification and Validation.

[S] — Shewchuk JR. An introduction to the conjugate gradient method without the agonizing pain.

3. Errors.

   Without a doubt, the manuscript contains some portion of typos. We would like to thank Shuhao Cao for spotting a few. We will cover all of them one by one:

   The approach in Gupta et al. (2021) is not transformer-based.

   It is not stated in our paper that the approach in Gupta is transformer-based. The whole sentence reads “This interest led to the construction of a plethora of novel efficient architectures including Fourier Neural Operator (FNO) Li et al. (2020), Deep Operator Network (DeelONet) Lu et al. (2019), wavelet- and transformer-based approaches Gupta et al. (2021), Hao et al. (2023).” This means Gupta et al. (2021) provide a wavelet-based approach (WNO) and Hao et al. (2023) — a transformer based.

   $a_{ij}\geq c > 0$ is not the sufficient condition that equation (6) is elliptic.

   This is a classic linear algebra shortcut for comparing positive-definite matrices. So what is written in plain English is “$a$ is a positive definite matrix with eigenvalues uniformly bound from below by constant $c$ that is greater than zero”. One can find that we are aware of how to properly define elliptic PDE by, for example, consulting the appendix. In equation (17) for one of the test problems, we generated positive definite matrices with uniformly bound eigenvalues using Cholesky decomposition.

   PDE is usually defined on open domains, not closed ones.

   We did not define PDE on the closed domain. Part $x\in\Gamma = [0, 1]^{D}$ is a definition of the domain itself.

   It is "Cauchy–Schwarz", not "Cauchy–Schwartz"

   We will properly spell Schwarz in the revision.

   equation (10) is wrong

   There is a minor typo in the equation. It read $\\sqrt{2\left(J(\widetilde{u}) - J(\widetilde{u} + w)\right)} \leq |||\widetilde{u} - u_{\text{exact}}|||^2$ and it should be $\\sqrt{2\left(J(\widetilde{u}) - J(\widetilde{u} + w)\right)} \leq |||\widetilde{u} - u_{\text{exact}}|||$. This equation provides a lower bound and is not used directly to obtain any results of the paper. This equation will likely be removed in the next revision to make space for other explanations.

   In summary, we would like to thank the reviewer for vigilance. In general, typos are very minor and unimportant. Still, it is good to spot as much of them as possible. We will also make careful proofreading at the later stages to make sure to remove the rest of the typos.

4. Benchmarks.

   We have two comments on the benchmarks proposed by Shuhao Cao:

   1. We performed a few additional tests: an elliptic equation with corner singularity (L-shape domain) and a posteriori functional error estimate for the convection-diffusion equation. Please, see https://openreview.net/forum?id=z62Xc88jgF&noteId=Wy4TDErd0T for more details.
   2. We do not find the standard FEM a posteriori error benchmarks helpful. There are two reasons for that:

      1. Neural PDE solvers and classical methods differ significantly.

         In classical methods, one can directly control the degree of approximation by increasing the number of degrees of freedom (e.g., the number of elements or points in physical space). Because of that classical benchmark aims to show how the situation changes for less smooth functions or domains with singularities. Often the focus is on the asymptotic properties (does the estimate converge) and on the adaptive sampling of grid points. All that is pointless for physics-informed neural networks and neural operators. The latter one does not show convergence when the grid is refined (see Figure 1 in [2]). The best result for the physics-informed neural network is often $10^{-2}$ relative error in $L_2$ norm with no sign of improvement for larger networks or finer grids (see [3] for SOTA results table 3; also table 3 from [4]). This absence of convergence makes classical benchmarks unsuitable. Indeed as shown in https://openreview.net/forum?id=z62Xc88jgF&noteId=Wy4TDErd0T nothing essentially new is observed for the L-shaped domain.

      2. We use a functional a posteriori error estimate that was already tested for standard problems. We merely incorporate the known results into the deep learning approach. It is not clear why standard tests are warranted in this case.

Overall, we would like to thank Shuhao Cao for spending time on reading the paper and making valuable suggestions.

[1] — https://arxiv.org/abs/1911.02109

[2] — https://arxiv.org/abs/2309.15325

[3] — https://arxiv.org/abs/2308.08468

[4] — https://arxiv.org/abs/2306.08827

We would like to thank Shuhao Cao for providing public comment with suggestions. We will treat this comment as an additional review given the background and keen interest of the author.

There are several issues raised by Shuhao Cao:

1. Paper [1] is claimed to contain similar results.

   We respectfully disagree that this contribution contains Astral loss defined in the present publication. To help other reviewers make up their minds, we provide a list of differences between our work and [1]:

   1. In [1] authors consider a neural network that takes the set of points $\mathbb{R}^{n}$ as an input and produces $\mathbb{R}^{n}$ solution values. We consider physics-informed neural networks and neural operators, i.e., our formulation is continuous.
   2. Equation (3.6) pointed out by Shuhao Cao is a discrete version of the loss function that can not be directly compared with (5) — continuous Astral loss defined for abstract PDE with known functional error estimate. It is more appropriate to compare FOSL functional (2.6) with concrete Astral loss for the elliptic equation given in (9). Undoubtedly, they are related since both of them reformulate elliptic problems in terms of fluxes (usual practice for a posteriori estimate) but there are a few significant differences:
      1. The claim that they are the same is exaggerated. For example, the terms with fluxes are different. Besides, Astral loss contains additional parameter $\beta$ that enters the optimization problem and the upper bound on Friedrichs's constant $C$.
      2. The FOSL loss is a posteriori error estimate and Astral is an error majorant. This means we know Astral always bound error from above — the property that we highlight through the text of the paper. For the FOSL loss, this is not the case.
   3. In paper [1] the results are given only for ordinary differential equations (boundary value problems) in a non-parametric setting. In our paper, we consider parametric and non-parametric PDEs.

2. Missing references.

   More specifically Shuhao Cao claims

   There are lots of missing references here on the newer development since Repin. Too bad it seems that none of the reviewers have ever worked in the area of a posteriori error estimation for traditional discretization methods (FEM or FVM or dG).

   We mostly reference Repin since he provided pioneering contributions to the field of a posteriori functional error estimate. We acknowledge that the literature on the a posteriori error estimate for classical discretization methods is vast and reach. However, we see no need to introduce additional references since they are largely irrelevant to our work. The functional error estimate for the elliptic equation derived by Repin remains correct (since it is supported with formal mathematical proof) even though there is “newer development” by other researchers.

   We also want to highlight that the functional a posteriori error estimate is markedly different from other approaches since it does not rely on discretization and smoothness properties. This is the main reason it is uniquely suitable for deep learning applications. It is very unfortunately this approach is not widely known by deep learning practitioners. We hope that our paper will improve the situation.

[1] — https://arxiv.org/abs/1911.02109

We can propose a few plausible explanations for that:

1. Astral loss is merely an upper bound on the error energy norm, not the energy norm itself. If the neural network predicts poor certificates, it can affect the quality of the solution since the upper bound will be loose.
2. Error norm contains derivatives, so the smoother the solution, the better. Residual loss heavily penalizes smoothness since it uses a strong form, i.e., it has a second derivative with respect to the network input. (This also means, it does not work for problems where only a weak form is available.) On the other hand, Astral loss contains only first derivatives, so the obtained solution can be less smooth.

Thank you for the additional experiments. Could you please comment on why the residual loss achieved better results for the L-shaped domain? As the theory suggests, I would expect Astral to always give better energy norms.