Posterior Concentration of Bayesian Physics-Informed Neural Networks for Elliptic PDEs

Problem Statement
This paper addresses the statistical uncertainty quantification of Physics-Informed Neural Networks (PINNs) when solving a general class of elliptic partial differential equations (PDEs) with non-homogeneous Dirichlet boundary conditions. While standard PINNs produce point estimators by minimizing a residual-based loss function, they fail to quantify the uncertainty arising from noisy or limited data. The authors investigate the posterior concentration rate of a Bayesian formulation of PINNs, specifically analyzing how the posterior distribution over neural network solutions contracts around the exact solution as the sample size increases. The data model involves independent, noisy measurements of both the source term (interior domain) and the boundary values, representing an indirect observation problem where the solution itself is not directly observed.

Methodology
The authors adopt a Bayesian framework where the solution to the elliptic PDE is treated as a random variable.

• Model: They consider the elliptic equation $-\text{div}(A\nabla u) + V u = f$ in a domain $\Omega$ with boundary condition $u = g$ on $\partial\Omega$. The coefficients $A$ and $V$ are assumed smooth, and the source term $f$ and boundary data $g$ are observed with additive Gaussian noise.
• Neural Network Architecture: The solution space is approximated by Deep Neural Networks (DNNs) using $\sigma_3$ activation functions (powers of ReLU). The network class includes a $C^2$-cutoff function to ensure uniform boundedness, which is critical for controlling statistical complexity.
• Prior Construction: A spike-and-slab prior is placed on the neural network weights. This prior models sparsity by assigning an activation probability to each neuron and includes an exponential prior on the amplitude of the weights, allowing for unbounded weights. The prior is designed to be adaptive to the unknown smoothness of the true solution.
• Loss Function: The PINN loss is defined as the sum of the $L^2$-norm of the PDE residual in the interior and the $L^2$-norm of the boundary residual. The empirical loss is constructed from the noisy observations.
• Theoretical Framework: The analysis relies on the general theory of posterior contraction rates for nonparametric Bayesian inference (Ghosal & Van der Vaart). The proof involves constructing a "sieve" (a subset of the function space with controlled complexity) and verifying three key conditions:
  1. Entropy Condition: Bounding the covering number of the sieve.
  2. Prior Mass Condition: Ensuring the prior assigns sufficient mass to a neighborhood of the true solution.
  3. Generalization Gap: Controlling the difference between the empirical loss and the population loss using Rademacher complexity arguments.

Key Contributions

1. Posterior Contraction Rates: The paper establishes that the posterior distribution of Bayesian PINNs concentrates around the exact solution at a near-minimax optimal rate. Specifically, if the true solution is $\beta$-Hölder smooth, the contraction rate is $\epsilon_n \asymp n^{-\frac{\beta-2}{d+2(\beta-2)}} (\log n)^{1/2}$ with respect to the norm induced by the PINN loss.
2. Rate Adaptivity: The proposed Bayesian PINN is rate-adaptive. It achieves the optimal contraction rate without requiring prior knowledge of the smoothness level $\beta$ of the exact solution. This is achieved through the specific construction of the spike-and-slab prior.
3. Minimax Lower Bound: The authors derive a minimax lower bound for the estimation error under the given data model (noisy interior and boundary measurements). They show that the derived upper bound for the Bayesian PINN matches this lower bound (up to logarithmic factors) when the number of interior and boundary measurements are equal, confirming the near-optimality of the method.
4. Handling Non-Homogeneous Boundaries: Unlike previous works that often assume homogeneous boundary conditions (where the solution vanishes on the boundary), this work handles non-homogeneous Dirichlet conditions, which is more representative of practical applications but statistically more challenging.

Results

• Theorem 3.1: Proves that for a $\beta$-Hölder smooth solution, the posterior probability of the event that the PINN loss exceeds $M_n^2 \epsilon_n^2$ converges to zero in probability as $n \to \infty$.
• Corollary 3.2: Establishes that the posterior mean also achieves the same convergence rate in the PINN loss.
• Remark 3.3: The authors clarify that the contraction in the PINN loss implies contraction in the $H^{1/2}(\Omega)$ norm. However, due to the use of an $L^2$-penalty on the boundary, they cannot currently guarantee convergence in stronger norms (e.g., $H^s$ for $s > 1/2$) without stronger boundary penalties or different variational formulations.
• Theorem 3.5: Provides the minimax lower bound, showing that any estimator (including PINNs) suffers from the curse of dimensionality under the assumption of Hölder continuity, consistent with existing literature.

Significance and Claims
The paper claims to provide the first rigorous statistical guarantees for the uncertainty quantification of PDEs via Bayesian PINNs in a general elliptic setting. By bridging the gap between the empirical success of Bayesian PINNs and theoretical understanding, the work demonstrates that these methods are not only flexible but also statistically sound, achieving near-optimal rates of convergence. The authors emphasize that their results hold for indirect observation models (noisy source and boundary data) rather than direct solution observations, which represents a harder statistical problem.

Limitations and Future Directions
The authors modestly note several limitations:

• The contraction rates are established in the $H^{1/2}$ norm (induced by the PINN loss) rather than stronger norms like $H^2$, which is a known challenge for non-homogeneous boundary conditions with $L^2$ penalties.
• The current proof techniques for posterior contraction do not directly extend to cases with highly imbalanced numbers of interior and boundary measurements, although the lower bound suggests optimal rates might still hold.
• The analysis is restricted to Hölder smooth solutions; extending these results to irregular solutions (e.g., Besov spaces) where adaptive methods are typically required remains an open question.