Refined generalization analysis of the Deep Ritz Method and Physics-Informed Neural Networks

We sincerely appreciate your time in reviewing our manuscript and your valuable insights. Below, we address your questions point by point.

Q1: Extension of theories to more complex high-dimensional problems.

A1: In this work, we analyze two scenarios: (1) when the solutions of PDEs reside in Barron spaces, and (2) when they belong to Sobolev spaces. For Sobolev spaces, certain constants in our results demonstrate exponential dependence on dimensionality, which limits their applicability to high-dimensional problems. This limitation is inevitable because the complexity of Sobolev spaces grows dramatically with increasing dimensions, leading to exponential scaling of certain constants in neural network approximation errors. In contrast, for Barron spaces, the constants in the generalization error exhibit at most polynomial dependence on dimensionality, making our results meaningful for high-dimensional settings.

To establish fundamental theoretical insights, this work primarily analyzes the settings where solutions in Barron or Sobolev spaces. The extension to more complex spaces is an important direction for future research. We conjecture that our methodology can be generalized to these settings with appropriate modifications.

Q2: Extension to other types of boundary conditions and nonlinear PDEs.

A2: Regarding the extension to other boundary conditions, such as the following mixed boundary conditions: $-\Delta u+Vu =f \ in \ \Omega, \ u+\beta\frac{\partial u}{\partial n}=g \ on \ \partial \Omega,$ where $\beta\neq 0$. A direct analysis shows that the solution $u^{* }$ of this equation satisfies that $u^{* }=\mathop{\arg\min} _{u \in H^1(\Omega) } \mathcal{E}(u):=\mathop{\arg\min} _{u \in H^1(\Omega) } \int _{\Omega} \left(||\nabla u||^2+Vu-2fu \right)dx+ \frac{1}{\beta} \int _{\partial \Omega} \left(u^2-2gu \right)ds.$ Moreover, for any $u\in H^1(\Omega)$, we have $\mathcal{E}(u)-\mathcal{E}(u^{* })\lesssim ||u-u^{* }|| _{H^1(\Omega)}^2\lesssim \mathcal{E}(u)-\mathcal{E}(u^{* }),$ which is similar to the strong convexity property required in our analysis for Neumann boundary conditions (see equation (4) and (6)), suggesting that our method can be naturally extended to mixed boundary conditions. For further details, please refer to Section D.2 (line 2851), where we discuss other boundary conditions.

Regarding the extension to nonlinear PDEs within the PINNs framework, our method can also generalize to a broad class of nonlinear PDEs. For instance, let us consider the following equation: $\mathcal{D}(u)=f(u) \ in \ \Omega, \ u=g \ on \ \partial \Omega,$ where $\mathcal{D}$ is a linear differential operator and $f$ is the Lipschitz nonlinear term. Here, for simplicity, we denote the loss function as $L(u)=||\mathcal{D}(u)-f(u)|| _{L^2(\Omega)}^2+||u-g|| _{L^2(\partial \Omega)}^2.$ Assume that $u^{* }$ is the true solution, then we can deduce that

$

L(u)&=||\mathcal{D}(u)-f(u)-\mathcal{D}(u^{* })-f(u^{* })|| _{L^2(\Omega)}^2+||u-u^{* }|| _{L^2(\partial \Omega)}^2 \\ &\lesssim ||\mathcal{D}(u)-\mathcal{D}(u^{* })|| _{L^2(\Omega)}^2 +||f(u)-f(u^{* })|| _{L^2(\Omega)}^2+ ||u-u^{* }|| _{H^1(\Omega)}^2\\ &\lesssim ||\mathcal{D}(u)-\mathcal{D}(u^{* })|| _{L^2(\Omega)}^2 +||u-u^{* }|| _{L^2(\Omega)}^2+ ||u-u^{* }|| _{H^1(\Omega)}^2.

$

Therefore, this bound is similar to that in the linear case, and our method remains applicable.

Moreover, we agree with the reviewer's suggestion. Extending this framework to broader physical problems (e.g., inverse problems) is indeed an important direction.

Q3: Practical applications.

A3: In this work, we primarily focus on generalization analysis and approximation rates, while optimization aspects are beyond the scope of this study. As you rightly pointed out, the computational resource and time requirements of neural network-based PDE solvers indeed limit their applicability to complex physical systems and large-scale problems. Developing computationally efficient optimization methods for these solvers with guaranteed accuracy remains a key focus of our ongoing research.

Q4: Issues in viewing PINNs as a multi-task learning problem.

A4: From a theoretical perspective, formulating PINNs as a multi-task learning problem allows us to derive tighter generalization bounds through the use of local Rademacher complexity in the multi-task setting. In the experimental perspective, as you mentioned, interactions between different tasks may indeed lead to certain issues, since the loss functions contain multiple different additive terms that can disagree and yield conflicting update directions. Some recent studies, such as [1] and [2], have developed efficient algorithms to mitigate such conflicts.

References:

[1]: Config: Towards conflict-free training of physics informed neural networks. ICLR 2025.

[2]: Dual cone gradient descent for training physics-informed neural networks. NeurIPS 2024.

Thank you for taking the time to review our article and for your insightful comments. Let us address your concerns point by point.

Q1: Missing references.

A1: We appreciate the reviewer for pointing out this important literature. It has established rigorous analysis of PICNN on the sphere, bridging the gap in understanding PINNs on manifolds. In the revised version, we will cite it with proper discussion.

Q2: Experiments.

A2: Regarding why we did not include experiments, like most theoretical papers, we have chosen to omit optimization error. In future work, we plan to account for optimization error and conduct experiments to validate our theory. For experimental validation of the Deep Ritz Method (DRM) for high-dimensional PDEs in the setting of Barron spaces, the experiments in [1] demonstrate that even in 100-dimensional cases, the DRM achieves a relative error of less than $20$% for Poisson equation and static Schrödinger equation. Moreover, when comparing experimental results with theoretical bounds, [1] observed that the convergence rates in their generalization error may not be sharp. Therefore, their experiments can, to some extent, support our findings, as our theoretical setting aligns with theirs—while we provide better generalization bounds.

Q3: Assumptions for the solutions.

A3: As the article primarily focuses on theoretical aspects, we–like most existing works–initially consider simpler cases. Investigating the setting where solutions belong to more complex spaces remains an important direction for future research.

Q4: Extension to other PDEs.

A4: In our analysis for PINNs, we require that the expected loss function of PINNs can be controlled by certain Sobolev norms between the exact and approximate solutions. For time-dependent PDEs, such as $\partial_t u=\Delta u+f(u),$ where $u^{* }$ is the true solution and $f$ is the nonlinear term assumed to be Lipschitz. For simplicity we only consider the interior term, then the loss function is $L(u)=||\partial _t u-\Delta u-f(u)|| _{L^2(0,T;L^2(\Omega))}^2.$ Then we have

$

\begin{aligned} L(u)&=||\partial _t u-\Delta u-f(u)-\partial _t u^{* }+\Delta u^{* }+f(u^{* })|| _{L^2(0,T;L^2(\Omega))}^2\\ &\lesssim ||\partial _t u-\partial _t u^{* }|| _{L^2(0,T;L^2(\Omega))}^2+||\Delta u-\Delta u^{* }|| _{L^2(0,T;L^2(\Omega))}^2+||u-u^{* } || _{L^2(0,T;L^2(\Omega))}^2. \end{aligned}

$

Then combining our method with neural network approximation results for spatiotemporal Sobolev spaces can yield similar generalization bounds.

Q5: Comparison with traditional methods.

A5: The error bounds of finite element method depend on the maximum mesh diameter $h$ and the computational complexity scales as $O(h^{-d})$ in high dimensions. Moreover, the error constants grow significantly with increasing dimension $d$, leading to the curse of dimensionality.

In this work, under the assumption that the solutions belong to Barron spaces, the relevant constants demonstrate at most polynomial dependence on the dimension, making the generalization bounds meaningful even in high-dimensional settings.

Q6: The use of $\text{ReLU}^k$.

A6: The choice of activation functions is just to demonstrate that our method can achieve better results. Our method also works for other activations. To illustrate this, we provide a concrete example using tanh for the static Schrödinger equation.

For the DRM, the generalization bound takes the form $\frac{\alpha}{n}+\epsilon _{app}^2,$ where $n$ is the sample size, $\alpha$ represents the neural network complexity and $\epsilon _{app}$ is the approximation error. Theorem 5.1 of [2] shows that for the solution $u _{S}^{* }\in W^{k,\infty}$ of the static Schrödinger equation, there exists a two-hidden-layer tanh network with width $O(N^d)$ achieving an approximation error $O(N^{-(k-1)})$. In this case, we have $\alpha=O(N^d)$, then taking a proper $N$ yields $\frac{\alpha}{n}+\epsilon _{app}^2=\frac{N^d}{n}+\frac{1}{N^{2(k-1)}}=O(n^{-\frac{2k-2}{d+2k-2}}).$ This result exactly coincides with the derived rate in Theorem 2.5(2) of our work.

Q7: Optimality.

A7: Some studies have shown that neural network-based estimators can achieve the minimax optimal rates for regression problems. These upper and lower bounds are estimated under the $L^2$ norm. However, for PINNs, different PDEs may require distinct norms to measure the discrepancy between approximate and exact solutions, which significantly differs from the framework of regression problems. Therefore, whether the improved generalization bounds for PINNs derived in this work are indeed optimal still requires further validation in future studies.

References:

[1]: A Priori Generalization Analysis of the Deep Ritz Method for Solving High Dimensional Elliptic Partial Differential Equations. COLT 2021.

[2]: On the approximation of functions by tanh neural networks. Neural Networks 2021.

We are grateful for your thorough review and thoughtful suggestions. Below we provide detailed responses to each of your comments.

Q1: Missing experiments.

A1: We agree that numerical validation would enhance the work. About experimental validation, the experiments in [1] have demonstrated that even in 100-dimensional cases, the Deep Ritz Method (DRM) attains a relative error of less than 20% for the Poisson equation and static Schrödinger equation in the setting of Barron space. Furthermore, when comparing experimental results with theoretical bounds, [1] observed that the convergence rates in their generalization bounds may not be sharp. Consequently, the experimental findings in [1] provide partial support for our results, as our theoretical framework aligns with theirs—while offering improved generalization bounds.

Regarding why we did not include experiments: Like most theoretical papers, we have chosen to omit optimization error. In future work, we plan to account for optimization error and conduct experiments to validate our findings.

Q2: Consider $H^2$ conforming networks for Deep Ritz Method (DRM) .

A2: We totally agree with you. In fact, if we assume that the solutions belong to Barron space $\mathcal{B}^s(\Omega)$ with $s\geq 3$, then $\text{ReLU}^2$ activation functions (i.e. $H^2$ conforming networks) could be employed for DRM, and our theoretical framework Theorem 2.4&2.5 still remains valid.

In this work, both DRM and PINNs are analyzed under the weakest solution regularity assumption -- specifically, we assume the solutions reside in relatively low-order Barron spaces (e.g. $\mathcal{B}^3(\Omega)$ in Theorem 3.4). This fundamental assumption governs our choice of activation functions throughout the paper. For example, for solutions in $\mathcal{B}^s(\Omega)$ with $s\geq 4$, $\text{ReLU}^3$ activations would become feasible for PINNs, and our theoretical framework Theorem 3.4 remains valid.

Q3: Regarding the positioning of this work.

A3: We can understand the reviewer’s concern about whether our purely theoretical error analysis fits ICML’s scope, especially given its focus on theoretical contributions without numerical simulations. However, we believe this work aligns with ICML’s mission to advance machine learning theory, as it addresses certain challenges arising in physics-informed learning where commonly used ML tools fail or are insufficient. Whereas papers in mathematical journals often rely on numerical integration with well-established error estimates, our work tackles the more complex setting of Monte Carlo methods for computing high-dimensional integrals arising in neural network loss functions---a setting that requires different theoretical tools. Specifically, we contribute by (1) developing a new theoretical framework for physics-informed learning and (2) deriving meaningful error bounds for Barron spaces in the overparameterized regime, advancing the theoretical foundations of machine learning for scientific computing domain.

Q4: Missing references.

A4: We sincerely appreciate the reviewer for pointing out this important reference. It indeed provides valuable insights, particularly for the second-order elliptic equations where Lemma 1 parallels our Lemma C.11. Additionally, this reference provides a thorough analysis for other types of PDEs, which broadens the scope. In the revised manuscript, we will cite this reference and incorporate this discussion on page 8 to offer readers a more comprehensive perspective.

Q5: Extension to more general activation functions.

A5: While our current analysis primarily focuses on $\text{ReLU}^k$ activation functions, the framework can indeed be extended to more general activations. Here we outline a concrete example using tanh activations for the static Schrödinger equation (Theorem 2.5(2)):

For the Deep Ritz method, the generalization bound under this setting takes the form: $\frac{\alpha}{n}+\epsilon _{app}^2,$ where $n$ is the sample size, $\alpha$ represents the neural network complexity and $\epsilon _{app}$ is the approximation error in the $H^1$ norm. When employing tanh activations, Theorem 5.1 of [2] demonstrates that for the solution $u _{S}^{* }\in W^{k,\infty}$ of the static Schrödinger equation, there exists a two-hidden-layer tanh network with width at most $\mathcal{O}(N^d)$ achieving an approximation error $\epsilon _{app}=\mathcal{O}(\frac{1}{N^{k-1}})$. In this setting, we have $\alpha=\mathcal{O}(N^d)$. Then taking a proper $N$ yields: $\frac{\alpha}{n}+\epsilon _{app}^2=\frac{N^d}{n}+\frac{1}{N^{2(k-1)}}=\mathcal{O}(n^{-\frac{2k-2}{d+2k-2}}).$ This matches the convergence rate derived in Theorem 2.5(2) of our work.

References:

[1]: A Priori Generalization Analysis of the Deep Ritz Method for Solving High Dimensional Elliptic Partial Differential Equations. COLT 2021.

[2]: On the approximation of functions by tanh neural networks. Neural Networks 2021.

We sincerely appreciate your time in reviewing our manuscript and your valuable insights. Below we provide a point-by-point response to your concerns.

Q1: Experiments.

A1: We agree that experiments would further strengthen our work. In future work, we plan to account for optimization error and conduct experiments to validate our results. For experimental validation of DRM in the setting of Barron space, the experiments in [4] show that even in 100-dimensional cases, the DRM achieves a relative error of less than $20$%. Moreover, when comparing experimental results with theoretical bounds, [4] found that the convergence rates in their generalization error may not be sharp. Therefore, their experiments can, to some extent, justify our findings, as our theoretical setting aligns with theirs—but we provide better generalization bounds.

Q2: References.

A2: We sincerely appreciate the reviewer for pointing out these important references. In the revised version, we will follow the format of [1] and incorporate all suggested references [2-8] with proper discussion.

Q3: Lower bounds.

A3: What [1] has achieved better than our work is that they also derived lower bounds for both DRM and PINNs. Their results show that the bound for DRM is not minimax optimal, whereas that for PINNs is minimax optimal. However, the metric used in [1] to evaluate PINNs is the $H^2$ norm, which requires strong convexity assumption on PDEs and neural network functions to belong to $H_0^1$. Such assumptions appear too stringent.

Recent studies have shown that neural network-based estimators can achieve minimax optimal rates for regression problems. These bounds are estimated under the $L^2$ norm. However, for PINNs, different PDEs require distinct norms to measure the discrepancy between the empirical and true solutions, which differs significantly from the regression framework. Therefore, whether the bounds we derived for PINNs and DRM are truly optimal still requires further investigation in future research.

Q4: Comparison of proof strategies with [1].

A4: For the DRM, we also analyze the Barron space setting and derive novel generalization bounds that remain meaningful even in overparameterized regimes. This advancement surpasses the capabilities of [1]'s approach (see Proposition 2.7). For the Poisson equation, the variational form $\mathcal{E} _{P}(u)$ is not equal to the expectation of its corresponding empirical part $\mathcal{E} _{n,P}(u)$. This makes the local Rademacher complexity (LRC, [3]) and the method in [1] infeasible. Morover, for example, the core lemma (Lemma B.4) in [1] requires the condition (from [3]) that $R _n({f\in \mathcal{F}: \mathbb{E}[f]\leq r}) \leq \phi(r),$ where $R _n$ is the empirical Rademacher complexity. Then in [1], an appropriate functions class $\mathcal{F}$ can be chosen such that $\mathbb{E}[f] = \mathcal{E} _{P}(u)-\mathcal{E} _{P}(u^{* })$. By doing so, the strong convexity can be used. However, for the Poisson equation, there does not exist a function class $\mathcal{F}$ such that $\mathbb{E}[f]=\mathcal{E} _{P}(u)-\mathcal{E} _{P}(u^{* })$. In contrast, this work develops a novel peeling technique to establish fast rates. The static Schrödinger equation shares similar strong convexity properties with certain classical problems like bounded-noise regression [2]. Although approaches from [1], [2] and our Poisson equation method remain applicable, they are much complicated. To address this, we instead develop a novel error decomposition method that allows direct application of the results in [3] (LRC), yielding better generalization bounds through a more concise proof framework.

For PINNs, [1] only considered the case of the static Schrödinger equation and assumed that the neural network function class is a subset of $H_0^1$. This assumption makes PINNs contain only interior terms, thereby giving them a strong convexity property similar to that of the DRM. Thus, the approach for PINNs in [1] is identical to that for the DRM. In contrast, we consider more general PDEs and treat PINNs as a multi-task learning (MTL) problem. The key difference from the DRM is that for PINNs, we only require a non-exact form of oracle inequality, eliminating the need for the strong convexity. Then, by using LRC under MTL (where we also derive an improved Talagrand-type concentration inequality for MTL), we obtain sharper generalization bounds. Moreover, this approach can also be extended to other types of PDEs and neural-network-based PDE-solving methods similar to PINNs.

References:

[1] Machine learning for elliptic pdes: Fast rate generalization bound, neural scaling law and minimax optimality. ICLR,2022.

[2] Deep neural networks for estimation and inference. Econometrica,2021.

[3] Local Rademacher complexities. AoS,2005.

[4] A Priori Generalization Analysis of the Deep Ritz Method for Solving High Dimensional Elliptic Partial Differential Equations. COLT,2021.