remarktheorem remark exampletheorem example

# On the estimation rate of Bayesian PINN for inverse problems

Yi Sun Department of Mathematics and Statistics, Boston University, Boston, MA
<EMAIL_ADDRESS>Debarghya Mukherjee Department of Mathematics and Statistics,
Boston University, Boston, MA<EMAIL_ADDRESS>Yves Atchadé Department of
Mathematics and Statistics, Boston University, Boston, MA<EMAIL_ADDRESS>
###### Abstract

Solving partial differential equations (PDEs) and their inverse problems using
Physics-informed neural networks (PINNs) is a rapidly growing approach in the
physics and machine learning community. Although several architectures exist
for PINNs that work remarkably in practice, our theoretical understanding of
their performances is somewhat limited. In this work, we study the behavior of
a Bayesian PINN estimator of the solution of a PDE from $n$ independent noisy
measurement of the solution. We focus on a class of equations that are linear
in their parameters (with unknown coefficients $\theta_{\star}$). We show that
when the partial differential equation admits a classical solution (say
$u_{\star}$), differentiable to order $\beta$, the mean square error of the
Bayesian posterior mean is at least of order $n^{-2\beta/(2\beta+d)}$.
Furthermore, we establish a convergence rate of the linear coefficients of
$\theta_{\star}$ depending on the order of the underlying differential
operator. Last but not least, our theoretical results are validated through
extensive simulations.

###### keywords:

partial differential equations, inverse problems, physics-informed neural
network, nonparametric regression, Bayesian deep learning, posterior
contraction

##  1 Introduction

We consider a class of inverse problems that consists in estimating the
parameters of a partial differential equation from noisy measurements of the
solution. Specifically, let $\Omega\subset\mathbb{R}^{m}$ be a bounded open
domain with a smooth boundary equipped with a probability measure $\nu$. For
$j\geq 1$, let $L^{2}(\Omega,\mathbb{R}^{j},\nu)$ denote the $L^{2}$-space of
$\mathbb{R}^{j}$-valued functions on $\Omega$, i.e.

$L^{2}(\Omega,\mathbb{R}^{j},\nu)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\left\\{f:\Omega\mapsto{\mathbb{R}}^{j}:\int_{\Omega}\|f(x)\|_{2}^{2}\
\nu(\mathrm{d}x)<\infty\right\\}\,.$

We consider a differential equation with parameter $\theta\in\mathbb{R}^{d}$
that aims to find a smooth function $u\in L^{2}(\Omega,\mathbb{R},\nu)$
satisfying

(1.1) $\mathcal{H}_{0}u+\theta^{\textsf{T}}\mathcal{H}_{1}u=f,\;\mbox{ with
initial/boundary condition }\;\;\mathcal{B}u={\bf 0},$

for some known function $f\in L^{2}(\Omega,\mathbb{R},\nu)$, and a pair of
differential operators $\mathcal{H}_{0}:\;L^{2}(\Omega,\mathbb{R},\nu)\to
L^{2}(\Omega,\mathbb{R},\nu)$ and
$\mathcal{H}_{1}:\;L^{2}(\Omega,\mathbb{R},\nu)\to
L^{2}(\Omega,\mathbb{R}^{d},\nu)$. The operator $\mathcal{H}_{0}$ is typically
(although not necessarily) the ”time” partial derivative, and the operator
$\mathcal{B}:\;L^{2}(\Omega,\mathbb{R},\nu)\to
L^{2}(\Omega,\mathbb{R}^{d_{1}},\nu)$ imposes the boundary/initial conditions.
A large class of differential equations (linear and nonlinear) that are linear
in their parameters can be written in this form. A classical example of such
an inverse problem arises from the heat equation as elaborated below:

######  Example  (Heat equation).

Given $L>0,T>0$, and parameter $\theta\in\mathbb{R}$, consider the heat
equation $u_{t}-\theta u_{xx}=0$ on $(0,L)$, with time domain $(0,T)$, initial
condition $u(0,\cdot)=g(\cdot)$, and boundary condition
$u(\cdot,0)=u(\cdot,L)=0$. Here $u_{t}$ (resp. $u_{xx}$) denotes the partial
derivative of $u$ with respect to $t$ (resp. the second order partial
derivative of $u$ with respect to $x$). With $\Omega=(0,T)\times(0,L)$, this
equation can be written as $\mathcal{H}_{0}u-\theta\mathcal{H}_{1}u=f$, where
$\mathcal{H}_{0}u=u_{t}$, $\mathcal{H}_{1}u=u_{xx}$, and $f=0$. The
boundary/initial condition operator is given by

$\mathcal{B}u=\left(\begin{array}[]{c}u_{\\{t=0\\}}-g\\\ u_{\\{x=0\\}}\\\
u_{\\{x=L\\}}\end{array}\right).$

where $u_{\\{t=0\\}}$ is the map $u_{\\{t=0\\}}(x)=u(0,x)$, with similar
definition for $u_{\\{x=0\\}}$ and $u_{\\{x=L\\}}$. Here $d=1$, $d_{1}=3$.

Let $(\theta_{\star},u_{\star})$ be a tuple that satisfies the PDE (1.1). We
assume to observe some noisy measurements of $u_{\star}$ as noted below:

######  H 1.

We have $n$ i.i.d. random locations and observations $({\bf
s}_{i},Y_{i})\in\Omega\times\mathbb{R}$, where

$Y_{i}\;|\;{\bf s}_{i}\stackrel{{\scriptstyle
ind.}}{{\sim}}\mathbf{N}\left(u_{\star}({\bf s}_{i}),\sigma^{2}\right),$

where $u_{\star}:\Omega\to\mathbb{R}$ is the unique solution of the pde (1.1)
with $\theta=\theta_{\star}$, for some unknown vector
$\theta_{\star}\in\mathbb{R}^{d}$. Here, $\textbf{N}(m,v^{2})$ denotes the
univariate Gaussian distribution with mean $m$ and variance $v^{2}$. The
variance parameter $\sigma^{2}$ is assumed to be known. Throughout we write
$\mathbb{P}$ for the joint distribution of $({\bf s}_{1},Y_{1}),\ldots,({\bf
s}_{n},Y_{n})$, and $\mathbb{P}_{n}$ for their corresponding empirical
measure.

Our goal is to estimate $(\theta_{\star},u_{\star})$ using the noisy
measurements. Over the last few years, physics-informed neural network (PINN)
has taken the numerical pde literature by storm (Raissi et al. (2019); Yu et
al. (2018); Sirignano and Spiliopoulos (2018); Sun et al. (2020), and Cuomo et
al. (2022) for an extensive review of the literature). In the setting of H1,
the approach consists of estimating $u_{\star}$ by regression while explicitly
using the information that the true function $u_{\star}$ is the solution of a
PDE. In this paper, as typically done in the PINN literature, we resort to a
sieve-based approach for estimating $u_{\star}$, i.e., we consider a sequence
of increasingly complex parametric models. Hence let
$\mathcal{F}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\\{u_{W},\;W\in\mathbb{R}^{q}\\}$
be a function class, where $u_{W}:\;\Omega\to\mathbb{R}$ is a function with
parameter $W\in\mathbb{R}^{q}$. Due to their superior empirical performances,
$u_{W}$ is typically taken as a deep neural network, and we follow that
practice, although our results can be applied more broadly. Given the data
$\mathcal{D}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\\{({\bf
s}_{i},Y_{i}),\;1\leq i\leq n\\}$, PINN in its frequentist formulation
estimates jointly $\theta$ and $W$ by minimizing the loss function

(1.2) $\frac{1}{2\sigma^{2}}\sum_{i=1}^{n}\left(Y_{i}-u_{W}({\bf
s}_{i})\right)^{2}+\frac{\lambda}{2}\left\\{\alpha_{1}|f-\mathcal{H}_{0}u_{W}-\theta^{\textsf{T}}\mathcal{H}_{1}u_{W}|_{2}^{2}+\alpha_{2}|\mathcal{B}u_{W}|_{2}^{2}\right\\},$

for some regularization parameters $\lambda\geq 0$, and
$\alpha_{1},\alpha_{2}\geq 0$, where $|\cdot|_{2}$ denotes the $L^{2}$-norm on
the appropriate function space. In this work, we approach the problem from a
Bayesian perspective. Hence, starting from a standard Gaussian prior for
$\theta\in\mathbb{R}^{d}$, and a sparsity inducing prior density $\Pi_{0}$ for
$W\in\mathbb{R}^{q}$ (see Section 1.5 below for our choice of $\Pi_{0}$), we
define the informative PINN prior distribution for $(\theta,W)$ as the
probability measure on $\mathbb{R}^{d}\times\mathbb{R}^{q}$ with density
proportional to

(1.3) $(\theta,W)\mapsto\Pi_{0}(W){\sf
exp}\left(-\frac{\lambda}{2}\left\\{\alpha_{1}|f-\mathcal{H}_{0}u_{W}-\theta^{\textsf{T}}\mathcal{H}_{1}u_{W}|_{2}^{2}+\alpha_{2}|\mathcal{B}u_{W}|_{2}^{2}\right\\}-\frac{1}{2}\|\theta\|_{2}^{2}\right).$

For appropriately large choices of $\lambda$, any realization $(\theta,W)$
from the PINN prior distribution (1.3) produces $(\theta,u_{W})$ that is
biased toward solving the pde (1.1). Rigorous general results (not specific to
PINN) of this flavor can be found, for instance, in Hwang (1980). When
$\theta$ is known, realizations from the conditional distribution of $W$ given
$\theta$ in (1.3) produce $u_{W}$ that approximately solve the pde (1.1) for
the given $\theta$. This corresponds to the initial PINN methodology of Raissi
et al. (2019).

For some convenience in the analysis, we will use the prior (1.3) with
$\alpha_{1}=1$, and $\alpha_{2}=0$. Since the parameter $\theta$ does not
appear in the boundary condition, the case $\alpha_{2}>0$ (although useful in
practice) does not induce any fundamentally new behavior, and at the expense
of more involved notations, our analysis can be easily modified to handle
$\alpha_{2}>0$. Given the observed data $\mathcal{D}$, the data generating
model postulated in (H1), and the PINN prior distribution (1.3), we thus
consider the posterior distribution

(1.4) $\Pi(\mathrm{d}\theta,\mathrm{d}W|\mathcal{D})\propto\\\ {\sf
exp}\left(-\frac{1}{2\sigma^{2}}\sum_{i=1}^{n}\left(Y_{i}-u_{W}({\bf
s}_{i})\right)^{2}-\frac{\lambda}{2}|f-\mathcal{H}_{0}u_{W}-\theta^{\textsf{T}}\mathcal{H}_{1}u_{W}|_{2}^{2}-\frac{\|\theta\|_{2}^{2}}{2}\right)\Pi_{0}(\mathrm{d}W)\mathrm{d}\theta.$

###  1.1 Main contributions

We study the behavior of the posterior distribution (1.4) as $n\to\infty$ (and
$q\to\infty$). We consider the case where the operator
$\mathcal{H}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}(\mathcal{H}_{0},\mathcal{H}_{1})$
is a differential operator that involves derivatives up to order $\tau$, for
some $\tau>0$. Under some additional regularity conditions, we establish in
Theorem 3.1 that when the true pde solution $u_{\star}$ is a classical
solution that possesses derivatives up to order $\beta\geq\tau$, the mean
square error of the posterior mean of $\theta$ under (1.4) satisfies (up to
log terms that we ignore),

(1.5)
$\mathbb{E}\left[\left\|\int_{\mathbb{R}^{d}}\vartheta\Pi^{(\theta)}(\mathrm{d}\vartheta|\mathcal{D})-\theta_{\star}\right\|_{2}^{2}\right]\lesssim
n^{-2(\beta-\tau)/(m+2\beta)},$

where $\Pi^{(\theta)}$ is the $\theta$-marginal of (1.4). In fact we obtain
this rate by showing that the 2-Wasserstein distance between $\Pi^{(\theta)}$
(the posterior distribution of $\theta$) and the Gaussian distribution
$\mathbf{N}(\theta_{\star},\lambda^{-1}\Sigma_{\star})$ converges to zero at
the rate given above (see Section 2 for the definition of $\Sigma_{\star}$).
Our simulation results suggest that, in general, this convergence does not
hold in the total variation metric. We also show in Theorem 2.2 that the pde
solution $u_{\star}$ is estimated at a rate at least as fast as the
nonparametric minimax optimal rate:

$\mathbb{E}\left[\int_{\mathbb{R}^{q}}\;|u_{W}-u_{\star}|_{2}^{2}\Pi^{(W)}(\mathrm{d}W|\mathcal{D})\right]\;\lesssim
n^{-2\beta/(m+2\beta)},$

where $\Pi^{(W)}$ denotes the $W$-marginal of (1.4). We conjecture that the
PINN posterior distribution actually converges faster than the nonparametric
minimax rate. Indeed, in a limiting case of an infinitely strong PINN prior,
we construct an estimator which, although computationally intractable in
general, achieves the parametric rate $n^{-1/2}$ (Theorem 2.3).

###  1.2 Related work

Despite their popularity, the theoretical aspects of PINN-related methods
remain under-studied. Most of the existing theoretical literature on PINN
deals with the forward problem of estimating $u$ for a given $\theta$, when
the $L^{2}$-norm $|\cdot|_{2}$ is replaced by a Monte Carlo approximation
(Mishra and Molinaro (2022); Shin et al. (2023); Lu et al. (2022)). The
objective in this line of work is to derive the mean square error of the PINN
solution as a function of the Monte Carlo sample size. More related to our
work, Mishra and Molinaro (2021) considers a PINN inverse problem where the
parameter of interest appears in the boundary condition. They work mainly in a
noiseless (or small-magnitude noise) data regime. Furthermore, they assume a
continuous space observation model, which is substantially different from the
framework considered here.

Inverse problems as related to differential and partial differential equations
have a long history that predates the current PINN literature (Ramsay et al.
(2007); Xun et al. (2013); McGoff et al. (2015)). In particular, we note the
close similarity between the PINN methodology and the approach by Xun et al.
(2013). The estimation rate in pde-driven inverse problems has been
extensively studied in statistics in recent years (Spokoiny (2020); Giordano
and Kekkonen (2020); Nickl et al. (2020); Monard et al. (2021); de Hoop et al.
(2023)). However, none of these results are directly applicable to our setting
or PINN more generally, since in this literature, the forward map is typically
assumed known and is not parametrized as in PINN.

A Bayesian version of PINN was proposed and studied in (Yang et al. (2021))
from an empirical viewpoint, where the authors noted a certain robustness of
Bayesian PINN to measurement noise. However, no theoretical analysis is
proposed.

A related version of the inverse problem considered in this work assumes that
$\theta$ is high-dimensional and sparse (Both et al. (2021); Chen et al.
(2020)). The goal is then to estimate $\theta$ under a sparsity assumption.
These physics discovery problems can be viewed as high-dimensional versions of
the problem investigated here. Thus, our work serves as a stepping-stone
toward a theoretical understanding of high-dimensional physics discovery
inverse problems.

###  1.3 Outline of the paper

We close the introduction with a description of the function class
$\mathcal{F}=\\{u_{W},\;W\in\mathbb{R}^{q}\\}$, the prior distribution
$\Pi_{0}$ for $W$, and a description of the main notations that we use
throughout the paper. Section 2 deals with the estimation rate for
$u_{\star}$. A Bernstein-von Mises theorem for the marginal distribution of
$\theta$ is discussed in Section 3. Section 4 presents our numerical
illustrations of the theoretical findings.

###  1.4 The function class $\mathcal{F}=\\{u_{W}\\}$

Although our results apply widely, we follow the PINN literature and focus on
the case where the function class $\mathcal{F}$ is constructed using neural
networks. Let $\xi>0$ be the depth of the network. Let
$(p_{\xi},\ldots,p_{0})$ be a sequence of integers representing the sizes of
the layers of the network, with $p_{0}=m$, and $p_{\xi}=1$. For $1\leq
i\leq\xi$, let $\mathsf{A}_{i}:\;\mathbb{R}^{p_{i}}\to\mathbb{R}^{p_{i}}$ be a
component-wise application of a 1-Lipschitz function
$\mathsf{a}_{i}:\mathbb{R}\to\mathbb{R}$. For $B\in\mathbb{R}^{p_{i}\times
p_{i-1}}$, and $b\in\mathbb{R}^{p_{i}}$, we set

(1.6) $\Psi_{B,b}^{(i)}({\bf
z})\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\mathsf{A}_{i}(B{\bf
z}+b),\;\;\;{\bf z}\in\mathbb{R}^{p_{i-1}}.$

We consider functions $u_{W}$ of the form

(1.7) $u_{W}({\bf
x})=\Psi_{W_{\xi},w_{\xi}}^{(\xi)}\circ\cdots\circ\Psi_{W_{1},w_{1}}^{(1)}({\bf
x}),\;\;\;{\bf x}\in\mathbb{R}^{m},$

with parameter $W=(W_{\xi},w_{\xi},\ldots,W_{1},w_{1})$, where
$W_{i}\in\mathbb{R}^{p_{i}\times p_{i-1}}$, and $w_{i}\in\mathbb{R}^{p_{i}}$,
and where $f\circ g$ is the composition of $f$ with $g$. For convenience, and
by vectorization we view $W$ as an element of $\mathbb{R}^{q}$, where

$q\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\sum_{i=1}^{\xi}p_{i}(1+p_{i-1}).$

We assume that the activation functions $\mathsf{a}_{\ell}$ are of class
$\mathcal{C}^{\infty}$, which implies that $u_{W}$ is of class
$\mathcal{C}^{\infty}$. The activation function $\mathsf{a}_{\xi}$ in the last
layer is typically taken as the identity function.

###  1.5 The prior distribution $\Pi_{0}$

We use sparsity to control the complexity of the function class
$\\{u_{W},\;W\in\mathbb{R}^{q}\\}$. We use a version of the spike-and-slab
prior distribution taken from Atchade and Bhattacharyya (2019). We set
$\mathcal{S}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\\{0,1\\}^{q}$. To
construct the prior $\Pi_{0}$, we first define a pair of random vectors
$(\Lambda,W)\in\mathcal{S}\times\mathbb{R}^{q}$ as follows. Let
$\Lambda_{j}\stackrel{{\scriptstyle
i.i.d}}{{\sim}}\textbf{Ber}((1+q^{\mathsf{u}+1})^{-1})$, $1\leq j\leq q$, for
some sparsity parameter $\mathsf{u}\geq 1$, where $\mathbf{Ber}(\alpha)$
denotes the Bernoulli distribution with parameter $\alpha\in(0,1)$. Given
$\Lambda\in\mathcal{S}$, the components of $W$ are independent with joint
density on $\mathbb{R}^{q}$ given by

$W\mapsto\prod_{j:\;\Lambda_{j}=1}\;\sqrt{\frac{1}{2\pi}}e^{-\frac{1}{2}W^{2}_{j}}\prod_{j:\;\Lambda_{j}=0}\sqrt{\frac{\rho_{0}}{2\pi}}e^{-\frac{\rho_{0}}{2}W_{j}^{2}},$

for some given parameter $\rho_{0}>1$. We denote $\bar{\Pi}_{0}$ the joint
distribution of $(\Lambda,W)$, and we let $\Pi_{0}$ (our prior on
$\mathbb{R}^{q}$) be the distribution of $W\odot\Lambda$, where for
$a,b\in\mathbb{R}^{q}$, $a\odot b$ denotes the component-wise product of $a$
and $b$. By construction, for any measurable function
$f:\mathbb{R}^{q}\to\mathbb{R}$, we have

$\Pi_{0}(f)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\int_{\mathbb{R}^{q}}f(v)\Pi_{0}(\mathrm{d}v)=\int_{\mathcal{S}\times\mathbb{R}^{q}}f(\Lambda\odot
W)\bar{\Pi}_{0}(\mathrm{d}\Lambda,\mathrm{d}W).$

###  1.6 Notation

For $\alpha\geq 1$, $\|\cdot\|_{\alpha}$ denotes the $\ell^{\alpha}$-norm on
finite-dimensional Euclidean spaces ($\mathbb{R}^{d}$, $\mathbb{R}^{n}$,
$\mathbb{R}^{d}$, etc, – which space should be clear from the context), and we
use $a^{\textsf{T}}b$ to denote the inner product between two finite
dimensional vector $a,b$. As usual, we allow $\alpha=+\infty$ (resp.
$\alpha=0$), by defining
$\|a\|_{\infty}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\max_{1\leq i\leq
d}|a_{i}|$ (resp. $|a|_{0}$ is the number of non-zero components of $a$).
Given a matrix $A$, $\lambda_{\textsf{min}}(A)$ (resp. $\|A\|_{\textsf{op}}$)
denotes its smallest (resp. largest) singular value.

Let $\Omega\subset\mathbb{R}^{m}$ be as above with a probability measure $\nu$
(typically the uniform measure on $\Omega$). Given an integer $j\geq 1$,
$L^{2}(\Omega,\mathbb{R}^{j},\nu)$ denotes the Hilbert $L^{2}$-space of
$\mathbb{R}^{j}$-valued functions with inner product

$\left\langle
f_{1},f_{2}\right\rangle\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\int_{\Omega}f_{1}({\bf
x})^{\textsf{T}}f_{2}({\bf x})\nu(\mathrm{d}{\bf x}),\;\;\mbox{ and norm
}\;\;|f|_{2}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\sqrt{\left\langle
f,f\right\rangle}.$

For $f:\;\Omega\to\mathbb{R}^{j}$, we set
$|f|_{\infty}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\sup_{x\in\Omega}\|f(x)\|_{\infty}$.

We use multi-index derivatives: given ${\bf k}=(k_{1},\ldots,k_{m})$, where
$k_{i}\geq 0$ is an integer,

$D^{\bf k}u(x)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\frac{\partial^{|{\bf
k}|}u(x)}{\partial x_{1}^{k_{1}}\ldots\partial x_{m}^{k_{m}}},$

where $|{\bf k}|=\sum_{i}k_{i}$. Given $\beta>1$, we set

$|u|_{\mathcal{C}^{\beta}}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\sum_{{\bf\alpha}:\;|{\bf\alpha}|\leq\lfloor\beta\rfloor}\;|D^{\bf\alpha}u|_{\infty}+\sum_{\alpha:\;|\alpha|=\lfloor\beta\rfloor}\;\sup_{x\neq
y}\frac{|D^{\bf\alpha}u(y)-D^{\bf\alpha}u(x)|}{\|y-x\|_{2}^{\beta-\lfloor\beta\rfloor}},$

where $\lfloor\beta\rfloor$ is the largest integer strictly smaller than
$\beta$. $|u|_{\mathcal{C}^{\beta}}$ is the Holder norm of $u$. We set
$\mathcal{C}^{\beta}(\Omega)$ (resp. $\mathcal{C}^{\beta}(\Omega,b)$) as the
set of all function $u:\;\Omega\to\mathbb{R}$ with finite Holder norm (resp.
with Holder norm bounded by $b$).

Throughout the paper, we use $C_{0},C_{1},C_{2},{\it etc.}$ to denote absolute
constants and $C$ to denote a generic absolute constant that depends on
$C_{0},C_{1},{\it etc.}$. The actual value of $C$ may change from one
appearance to the next. Similarly, $c_{0},c_{1},{\it etc.}$ denote problem-
dependent constants (that do not depend on the sample size $n$ or the model
size $q$). Specifically, these constants typically depend on the noise
parameter $\sigma^{2}$, the dimension $d$, and the true parameter
$\theta_{\star}$. We will also use $c$ to denote a generic constant that
depends on $c_{0},c_{1},{\it etc.}$, and that we do not track. The actual
value of $c$ may change from one appearance to the next.

##  2 Estimation rate of $u_{\star}$

We maintain throughout the basic assumption that the operator
$\mathcal{H}=(\mathcal{H}_{0},\mathcal{H}_{1})$ is a differential operator of
order $\tau>0$, in the sense that for all $\beta>\tau$, and for
$u,v\in\mathcal{C}^{\beta}(\Omega)$,

(2.1) $\left|\mathcal{H}u-\mathcal{H}v\right|_{2}\leq C_{0}\max_{{\bf
k}:\;|{\bf k}|\leq\tau}\;\left|D^{\bf k}u-D^{\bf k}v\right|_{\infty},$

for some absolute constant $C_{0}$.

Because the pde equation is linear in $\theta$, and the log-prior is quadratic
in $\theta$, it is straightforward to integrate out $\theta$ from the
posterior distribution (1.4). We first introduce some appropriate notations to
do this. For $W\in\mathbb{R}^{q}$, we define
$\bar{f}_{W}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}f-\mathcal{H}_{0}u_{W}$.
Let $\Phi_{W}\in\mathbb{R}^{d}$ be the vector with $j$-th component given by

$(\Phi_{W})_{j}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\left\langle\bar{f}_{W},(\mathcal{H}_{1}u_{W})_{j}\right\rangle,$

and where $(\mathcal{H}_{1}u_{W})_{j}$ denotes the $j$-component of
$\mathcal{H}_{1}u_{W}$. We define $\Phi_{\star}\in\mathbb{R}^{d}$ similarly,
replacing $u_{W}$ by $u_{\star}$. Let $\Sigma_{W}\in\mathbb{R}^{d\times d}$
given by

$(\Sigma_{W})_{j,k}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\left\langle(\mathcal{H}_{1}u_{W})_{j},(\mathcal{H}_{1}u_{W})_{k}\right\rangle\;+\;\frac{1}{\lambda}\textbf{1}_{\\{j=k\\}},\;\;1\leq
j,k\leq d.$

We define $\Sigma_{\star}$ similarly by replacing $u_{W}$ by $u_{\star}$.
Given $u\in L^{2}(\Omega,\mathbb{R},\nu)$ we define

(2.2) ${\cal
J}(u)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\min_{\theta\in\mathbb{R}^{d}}\;\left[|f-\mathcal{H}_{0}u-\theta^{\textsf{T}}\mathcal{H}_{1}u|_{2}^{2}+\frac{1}{\lambda}\|\theta\|_{2}^{2}\right].$

It is easy to check that ${\cal
J}(u_{W})=|\bar{f}_{W}|_{2}^{2}-\Phi_{W}^{\textsf{T}}\Sigma_{W}^{-1}\Phi_{W}$,
and we can integrate out $\theta$ from the posterior distribution (1.4) to
obtain the marginal distribution of $W$ given by

(2.3) $\Pi(W|\mathcal{D})\propto\frac{\Pi_{0}(W)}{\sqrt{{\sf
det}(\Sigma_{W})}}{\sf
exp}\left(-\frac{1}{2\sigma^{2}}\sum_{i=1}^{n}(Y_{i}-u_{W}({\bf
s}_{i}))^{2}-\frac{\lambda}{2}{\cal J}(u_{W})\right).$

At times, when the distinction is needed we will write $\Pi^{(W)}$ to denote
(2.3). We make the following basic assumption on the function class
$\mathcal{F}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\\{u_{W},\;W\in\mathbb{R}^{q}\\}$.

######  H 2.

  1. 1.

For all $W_{0}\in\mathbb{R}^{q}$, and $r>0$, we can find $L_{W_{0},r}\geq 1$
such that if $W_{1},W_{2}\in\mathbb{R}^{q}$ satisfy
$\max(\|W_{1}-W_{0}\|_{2},\|W_{2}-W_{0}\|_{2})\leq r$, then it holds

$\left|u_{W_{1}}-u_{W_{2}}\right|_{\infty}\leq
L_{W_{0},r}\;\|W_{1}-W_{2}\|_{2}.$

  2. 2.

For all $W\in\mathbb{R}^{q}$, $u_{W}$ has derivatives to the order $\tau$, and
there exists $\kappa>0$ and a constant $c_{1}$ such that for all
$W_{1},W_{2}\in\mathbb{R}^{q}$, such that
$|u_{W_{1}}-u_{W_{2}}|_{\infty}\leq\varepsilon$, we have

$\max_{{\bf k}:\;|{\bf k}|\leq\tau}\;\left|D^{\bf k}u_{W_{1}}-D^{\bf
k}u_{W_{2}}\right|_{\infty}\leq c_{1}\varepsilon^{\kappa}.$

######  Remark .

The local Lipschitz condition imposed in H2-(1) is known to hold for most deep
neural net architectures with Lipschitz activation functions (see e.g.
Proposition 6 of Taheri et al. (2021)). Specifically, if $u_{W}$ is a deep
neural network as constructed in Section 1.4 with 1-Lipschitz activation
functions, then with $W_{1},W_{2}\in\mathbb{R}^{q}$ satisfying
$\max(\|W_{1}-W_{0}\|_{2},\|W_{2}-W_{0}\|_{2})\leq r$, Proposition 6 of Taheri
et al. (2021) and some easy calculations yield

$|u_{W_{1}}-u_{W_{2}}|_{\infty}\leq L_{W_{0},r}\|W_{1}-W_{2}\|_{2},$

with

(2.4) $L_{W_{0},r}=\left(\sup_{{\bf x}\in\Omega}\|{\bf
x}\|_{2}\right)\;\sqrt{\xi}\;\left(1+\frac{(\|W_{0}\|_{2}+r)^{2}}{\xi}\right)^{\xi},$

where $\xi$ is the depth of the DL function.

######  Remark .

Since we aim to solve pdes, it is natural to require the function $u_{W}$ to
be smooth. H2-(2) further requires a Holder-type assumption on the
differential operator $D^{{\bf k}}$ over the function class $\mathcal{F}$.
This assumption is less common in the literature, and we check its validity in
Theorem 2.4.

Let $\epsilon:\;{\mathbb{N}}\to\mathbb{R}$ be a nonincreasing function such
that for all $s\geq 1$

(2.5)
$\min\left\\{|u_{W}-u_{\star}|_{\infty},\;W\in\mathbb{R}^{q}:\;\|W\|_{0}\leq
s\right\\}\leq\epsilon(s).$

$\epsilon(s)$ defines the approximation skills of the function class
$\\{u_{W}\\}$ at the sparsity level $s$, and is a nonincreasing function of
$s$. Of importance are sparsity levels $s>1$ at which the approximation error
$\epsilon(s)$ matches the statistical error:

(2.6) $\epsilon(s)\leq\sigma\sqrt{\frac{\mathsf{u}s\log(q)}{n}}.$

We make the following assumption.

######  H 3.

There exists $\epsilon:\;{\mathbb{N}}\to\mathbb{R}$ such that (2.5) holds.
Given $n$ and $q$, let $s_{0}\geq 2$ be the smallest integer satisfying (2.6),
and $\epsilon_{0}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\epsilon(s_{0})$.
We assume that there exists $W_{0}\in\mathbb{R}^{q}$, with $\|W_{0}\|_{0}\leq
s_{0}$ that satisfies $|u_{W_{0}}-u_{\star}|_{\infty}\leq\epsilon_{0}$, and

(2.7)
$\min\left(\lambda_{\textsf{min}}(\Sigma_{W_{0}}),\lambda_{\textsf{min}}(\Sigma_{\star})\right)\geq
C_{2},$

for some absolute constant $C_{2}$. Furthermore, we assume that
$\lim_{n\to\infty}\epsilon_{0}=0$, and

(2.8)
$(\|W_{0}\|_{\infty}+1)^{2}+\log\left(Ls_{0}^{1/2}/\epsilon_{0}\right)\leq
C_{3}\mathsf{u}\log(q),$

for some absolute constant $C_{3}$, where
$L\stackrel{{\scriptstyle\mathrm{def}}}{{=}}L_{W_{0},1}$ is as in H2.

######  Remark .

We recall that we aim to estimate $u_{\star}$ by the sieve method, i.e., we
approximate $u_{\star}$ by a collection of functions
$\mathcal{F}\equiv\mathcal{F}_{n}$, where the _size_ of the function class
typically grows with the sample size $n$. The assumption that $\epsilon_{0}\to
0$ as $n\to\infty$ in H3 highlights the relationship between the sample size
$n$ and the model size: given $n$, we assume that the model $\mathcal{F}$ is
chosen such that its approximation error at some sparsity level $s$ can match
(or be smaller than) the statistical error at sparsity level $s$.

There has been a flurry of research activities in recent years to derive
precise estimates on $\epsilon(s)$ for various DL function classes, using
various norms (Yarotsky (2017); Schmidt-Hieber (2020); Lu et al. (2021);
Belomestny et al. (2022)). For example for piecewise polynomial activation
functions, and $u_{W}$ as defined in Section 1.4, it is known (see e.g.
Yarotsky and Zhevnerchuk (2020); Lu et al. (2021)) that at the sparsity level
$s$, and depth $\xi=O(\log(s))$, $\mathcal{F}$ can approximate a function
$u_{\star}\in\mathcal{C}^{\beta}(\Omega)$ with precision $(1/s)^{\beta/m}$.
This yields

$\epsilon(s)\lesssim(1/s)^{\beta/m}.$

In that case, solving for $s_{0}$ in (2.6) yields

$s_{0}\sim\left(\frac{n}{\log(q)}\right)^{\frac{m}{m+2\beta}},\;\mbox{ and
}\;\;\epsilon_{0}\sim\left(\frac{\log(q)}{n}\right)^{\frac{\beta}{m+2\beta}}.$

In the high-dimensional setting considered here where $q$ is typically much
larger than $n$, (2.8) is typically satisfied.

Ultimately, we need some assumptions on the stability of the pde operator
$\mathcal{H}$ as presented below:

######  H 4.

With $\kappa$ as in H2-(2), there exists an absolute constant $C_{4}$ such
that for all $\varepsilon>0$, if $|u_{W}-u_{\star}|_{2}\leq\varepsilon$, then
$|\mathcal{H}u_{W}-\mathcal{H}u_{\star}|_{2}\leq C_{4}\varepsilon^{\kappa}$.

######  Remark .

H4 is similar to H2-(2) and both can be checked in a similar way if
$u_{\star}$ is a classical solution. Indeed, if the pde solution $u_{\star}$
is a classical solution, and $u_{\star}\in\mathcal{C}^{\beta}(\Omega)$ for
some $\beta>\tau$, then by (2.1), and Theorem 2.4 below, H4 holds with
$\kappa=1-\tau/\beta$.

Given $s\geq 0$, and constant $C$, we set
$\tau_{s}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\sqrt{(2+\mathsf{u})(1+s)\log(q)}$,
and

$\mathcal{F}_{s}=\mathcal{F}_{s,C}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\left\\{u_{W}\in\mathcal{F}:\;\|W\|_{0}\leq
s,\;\|W\|_{\infty}\leq\tau_{s},\;\;{\cal
J}(u_{W})\leq\frac{Cn\epsilon_{0}^{2}}{\lambda\sigma^{2}}\right\\},$

where ${\cal J}$ is as defined in (2.2). We show below that the PINN posterior
distribution puts high probabilities on the sets $\mathcal{F}_{s}$.
Specifically, the following theorem is proved as Lemma A.5 in Section A.2.

######  Theorem 2.1.

Assume H1-H4. Let Let $\epsilon_{0}$, $W_{0}$ be as in H3. We can find an
absolute constant $C$ such that with $s=C\|W_{0}\|_{0}$, it holds,

$\mathbb{E}\left[\Pi\left(\mathcal{F}_{s,C}|\mathcal{D}\right)\right]\geq
1-C_{0}e^{-\frac{n\epsilon_{0}^{2}}{2\sigma^{2}}},$

for some absolute constant $C_{0}$.

To highlight this point, we note that if a function $u\in
L^{2}(\Omega,\mathbb{R},\nu)$ solves the pde equation (1.1) (without the
initial/boundary condition) for some parameter $\theta$, say, then

${\cal
J}(u)=\frac{1}{\lambda}\theta^{\textsf{T}}\left(\mathbbm{1}-\frac{1}{\lambda}\Sigma_{u}^{-1}\right)\theta\leq\frac{\|\theta\|_{2}^{2}}{\lambda}.$

Hence with the size of $\theta$ remaining bounded, and for $\lambda$ large,
functions that approximately solve the pde equation (1.1) satisfies ${\cal
J}(u)\lesssim C/\lambda$ for some constant $C$. Theorem 2.1 thus implies that
for $\lambda$ large, the PINN posterior distribution inherits the inductive
bias of the prior and restricts the search of a solution for the nonparametric
regression problem on approximate solutions of the pde (1.1). The rate of
convergence of PINN, therefore, depends mainly on the complexity of the
function class $\mathcal{F}_{s}$, as measured, for instance, by its covering
number. The next assumption models the metric entropy of $\mathcal{F}_{s}$.

######  H 5.

Given $s\geq 2$, let $3\leq b<\infty$ be such that

(2.9) $\sup_{u\in\mathcal{F}_{s}}\;|u-u_{\star}|_{\infty}\leq b.$

There exists $V_{1}=V_{1}(s)$ and $V_{2}=V_{2}(s)\geq 6b$ such that for all
$0<\varepsilon\leq\tau_{s}$,

(2.10) $\log\mathcal{N}(\varepsilon,\mathcal{F}_{s},\|\cdot\|_{\infty})\leq
V_{1}\log\left(\frac{V_{2}}{\varepsilon}\right),$

where $\mathcal{N}(\varepsilon,\mathcal{F}_{s},\|\cdot\|_{\infty})$ denotes
the $\varepsilon$-covering number of $\mathcal{F}_{s}$ in the $L_{\infty}$
norm.

The following result establishes the rate of convergence of the PINN posterior
and is proved in Section A. We assume that the prior parameter $\lambda>0$ is
taken such that

(2.11) $\|\theta_{\star}\|_{2}^{2}\lambda\leq
C_{5}n\epsilon_{0}^{2(1-\kappa)}$

for some absolute constant $C_{5}$. We also impose the following mild
technical condition: there exists an absolute constant $C_{6}$ such that

(2.12)
$d\left(\log\left(n\|\Sigma_{\star}\|_{\textsf{op}}\right)+\|\theta_{\star}\|_{\infty}^{2}\right)\leq\frac{C_{6}n\epsilon_{0}^{2}}{2\sigma^{2}}.$

######  Theorem 2.2.

Assume H1-H5, (2.11), and (2.12). Then we can find absolute constants $C$ and
$M\geq 1$ such that with

(2.13) $s=Cs_{0},\;\mbox{ and
}\textsf{r}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}2(b+\sigma)\sqrt{\frac{V_{1}\log(V_{2}\sqrt{n})}{n}},$

where $V_{1}=V_{1}(s)$, $V_{2}=V_{2}(s)$ are as defined in H5, it holds for
all $n$ large enough,

$\mathbb{E}\left[\Pi\left(\left\\{W\in\mathbb{R}^{q}:\;|u_{W}-u_{\star}|_{2}>M\textsf{r}\right\\}\;|\;\mathcal{D}\right)\right]\leq
c\left(e^{-n\epsilon_{0}^{2}/(2\sigma^{2})}+e^{-V_{1}\log(V_{2}\sqrt{n})/C}\right),$

where $c$ is a constant that depends only on $c_{1},c_{2}$ and $\sigma^{2}$.

###  2.1 Minimax nonparametric rate

Bounding the metric entropy of the function class $\mathcal{F}_{s}$ in a way
that leverages the pde information ${\cal J}(u_{W})\leq
Cn\epsilon_{0}^{2}/(\lambda\sigma^{2})$ has proved challenging. However, we
note that $\mathcal{F}_{s}\subseteq\\{u_{W}\in\mathcal{F}:\;\|W\|_{0}\leq
s,\;\|W\|_{\infty}\leq\tau_{s}\\}$, and the metric entropy of this latter set
is straightforward to bound using H2-(1). See also (Schmidt-Hieber (2020); Lu
et al. (2021)). Indeed, for all $W,W^{\prime}$ in the ball
$\\{W\in\mathbb{R}^{q}:\;\|W\|_{0}\leq s,\;\|W\|_{\infty}\leq\tau_{s}\\}$,
$\max(\|W\|_{2},\|W^{\prime}\|_{2})\leq s^{1/2}\tau_{s}$. Hence by H2-(1),
$|u_{W}-u_{W^{\prime}}|_{\infty}\leq Ls^{1/2}\|W-W^{\prime}\|_{\infty}$, where
$L=L_{{\bf 0},s^{1/2}\tau_{s}}$. The size of all the
$(\varepsilon^{\prime}=\varepsilon/Ls^{1/2})$-cover of all the
$\tau_{s}$-balls of all the $s$-sparse subspaces of $\mathbb{R}^{q}$ in the
infinity-norm is bounded from above by

${q\choose s}\left(1+\frac{2r_{s}}{\varepsilon^{\prime}}\right)^{s}.$

It follows that

$\log\mathcal{N}(\varepsilon,\mathcal{F}_{s},\|\cdot\|_{\infty})\leq
s\log(q)+s\log\left(1+\frac{2s^{1/2}\tau_{s}L}{\varepsilon}\right),$

and H5 holds with

$V_{1}(s)=s,\;\;\mbox{ and }\;\;V_{2}(s)=qr_{s}\left(1+2s^{1/2}L\right)\leq
c\sqrt{\xi}q\left(1+s^{2}\log(q)\right)^{\xi+1/2},$

for some constant $c$ that depends only on $\mathsf{u}$, and
$\sup_{x\in\Omega}\|{\bf x}\|_{2}$, where $\xi$ denotes the depth of the
neural network class. Hence, with $s=Cs_{0}$,

$V_{1}(s)=s=Cs_{0},\;\mbox{ and
}\;\;\log(\sqrt{n}V_{2}(s))\lesssim\log(n)+\xi\log(q),$

where $\xi$ is the depth of the neural network. Furthermore, we have seen in
Remark 2 that if $u_{\star}$ is a classical solution and
$u_{\star}\in\mathcal{C}^{\beta}(\Omega)$, then for piecewise polynomial
activation at depth $\xi=O(\log(p))$, we have

$s_{0}\sim\left(\frac{n}{\log(q)}\right)^{\frac{m}{m+2\beta}}.$

As a result, up to log terms that we ignore, we obtain the estimation rate

$\textsf{r}\lesssim\sqrt{\frac{s_{0}}{n}}\lesssim\left(\frac{1}{n}\right)^{\frac{\beta}{m+2\beta}}.$

Hence, Bayesian PINN estimates $u_{\star}$ at least at the nonparametric
minimax rate $n^{-\beta/(m+2\beta)}$, that is, at the optimal rate achievable
when $u_{\star}$ is known only to be $\beta$-smooth (Stone (1982)) without any
PDE constraint. We note that our metric entropy calculations of the function
class $\mathcal{F}_{s}$ are rather crude, as they ignore the inductive bias
induced by the PINN prior. Thus, the PINN posterior distributions likely enjoy
a faster convergence rate than the abovementioned nonparametric minimax rate.
To understand this, consider the limiting case when $\lambda=+\infty$, i.e.,
the prior put all the mass on $\\{f\in L^{2}(\Omega,\mathbb{R},\nu):\;{\cal
J}(u)=0\\}$, that is, the set of functions that are solutions of the pde
equation $\mathcal{H}_{0}u+\theta^{\textsf{T}}\mathcal{H}_{1}u=f$ for some
$\theta\in\mathbb{R}^{d}$. In the subsequent section, we show that estimating
$u_{\star}$ at a parametric rate (up to a log factor) is possible in this
limiting case.

###  2.2 Estimation at the parametric rate

We assume that $\theta_{\star}$ belongs to a compact set
$\Theta\subset\mathbb{R}^{d}$, and for each $\theta\in\Theta$ the pde equation
(1.1) (including the initial/boundary condition) admits a unique solution
$u_{\theta}\in\mathcal{C}^{\beta}(\Omega)$ for some given and fixed $\beta$.
Furthermore we assume that the map
$\Gamma:\;\Theta\to\mathcal{C}^{\beta}(\Omega)$ that maps $\theta\to
u_{\theta}$ is Lipschitz. Specifically,

######  H 6.

There exists a constant $L>0$ such that:

$|\Gamma(\theta_{1})-\Gamma(\theta_{2})|_{\infty}\leq
L\|\theta_{1}-\theta_{2}\|_{2}\ \ \forall\ \
\theta_{1},\theta_{2}\in\Theta\,.$

Let
$\Gamma(\Theta)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\\{\Gamma(\theta),\;\theta\in\Theta\\}$.
Under the data generating process assumed in H1, we consider the constrained
empirical risk minimization estimator

(2.14)
$\hat{u}={\arg\min}_{u\in\Gamma(\Theta)}\frac{1}{n}\sum_{i=1}^{n}\left(Y_{i}-u({\bf
s}_{i})\right)^{2}\,.$

The following theorem is proved in Section A.7.

######  Theorem 2.3.

Assume H1, and H6. Then any solution $\hat{u}$ of (2.14) satisfies

$|\hat{u}-u_{\star}|_{2}^{2}=O_{p}\left(\frac{d}{n}\log{\left(\frac{n}{d}\right)}\right)\,.$

It is immediate from the above theorem that
$|\hat{u}-u_{\star}|_{2}=O_{p}(\sqrt{\frac{d}{n}})$ (up to a log factor) which
is the standard parametric rate for estimating a $d$ dimensional parameter.
The estimator $\hat{u}$ is typically not computable since it requires solving
the pde for each $\theta\in\Theta$. But the rate in Theorem 2.3, in light of
Theorem 2.1, suggests that the PINN posterior distribution may actually have a
convergence rate faster than the nonparametric minimax rate. More research is
needed on this issue.

Regarding Assumption 6, one may replace the $L_{2}$ norm in the assumption by
any $L_{p}$ norm without hurting the rate, but to translate the complexity of
$\Theta$ to $\Gamma(\Theta)$, we believe that some sort of smoothness
assumption is necessary. Whether one can do entirely without Assumption 6 is
an intriguing question that we leave for future research.

###  2.3 Checking Assumption H2-(2)

We end this section with a result showing that sufficiently smooth functions
typically satisfy H2-(2). The proof is given in Section A.6. We recall that
$\Omega$ is a bounded open subset of $\mathbb{R}^{m}$, and we write
$\mathbf{B}(x,\alpha)$ to denote the Euclidean ball of $\mathbb{R}^{m}$
centered at $x$ and with radius $\alpha$. We define
$\textsf{int}_{\alpha}(\Omega)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\\{x\in\Omega:\;\mathbf{B}(x,\alpha)\subset\Omega\\}$.
We recall also that $\mathcal{C}^{\beta}(\Omega,b)$ denotes the ball with
radius $b$ in the Holder space $\mathcal{C}^{\beta}(\Omega)$.

######  Theorem 2.4.

Given $\varepsilon>0$, $\beta>1$, integer $\tau\in[0,\beta)$, and $M>0$,
define

$\bar{\alpha}=\left(\frac{\tau\varepsilon\lfloor\beta\rfloor!}{2M(\beta-\tau)}\right)^{1/\beta}.$

For all $u,\tilde{u}\in\mathcal{C}^{\beta}(\Omega,M)$ such that
$\sup_{x\in\Omega}\;|u(x)-\tilde{u}(x)|_{\infty}\leq\varepsilon$, and for all
multi-index ${\bf k}=(k_{1},\ldots,k_{m})$, with $|{\bf k}|\leq\tau$, we have

$\sup_{x\in\mathsf{int}_{\bar{\alpha}}(\Omega)}\;\left|D^{\bf k}u(x)-D^{\bf
k}\tilde{u}(x)\right|_{\infty}\leq
CM^{\frac{\tau}{\beta}}\varepsilon^{\frac{\beta-\tau}{\beta}},$

for some constant $C$ that depends only on $\tau$ and $\beta$.

######  Remark 2.5.

We note that $\bar{\alpha}\downarrow 0$ as $M\uparrow\infty$. Hence if
$u,\tilde{u}\in\mathcal{C}^{\beta}(\bar{\Omega},M_{0})$ for some well-chosen
$\bar{\Omega}\supset\Omega$, and
$\sup_{x\in\bar{\Omega}}\;|u(x)-\tilde{u}(x)|_{\infty}\leq\varepsilon$ then by
taking $M\geq M_{0}$ large enough such that
$\Omega\subset\mathsf{int}_{\bar{\alpha}}(\bar{\Omega})$, we get

$\sup_{x\in\Omega}\;\left|D^{\bf k}u(x)-D^{\bf
k}\tilde{u}(x)\right|_{\infty}\leq
CM^{\frac{\tau}{\beta}}\varepsilon^{\frac{\beta-\tau}{\beta}}.$

##  3 A Bernstein-von-Mises theorem for the $\theta$-marginal

Given $\mu\in\mathbb{R}^{d}$, and $\Lambda\in\mathbb{R}^{d\times d}$ symmetric
and positive definite, we write $\textbf{N}(\mu,\Lambda)(\cdot)$ as the
probability measure of the Gaussian distribution $\textbf{N}(\mu,\Lambda)$.
Given $W\in\mathbb{R}^{q}$, let

(3.1)
$\hat{\theta}_{W}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\Sigma_{W}^{-1}\Phi_{W}.$

It is then easy to see from (1.4) that the conditional distribution of
$\theta$ given $W$ is $\mathbf{N}(\hat{\theta}_{W},\Sigma_{W}^{-1}/\lambda)$.
For clarity’s sake we will write $\Pi^{(\theta)}$ (resp. $\Pi^{(W)}$) to
denote the marginal distribution of $\theta$ (resp. $W$) under (1.4). We have

$\Pi^{(\theta)}(\cdot|\mathcal{D})=\int_{\mathbb{R}^{q}}\left[\mathbf{N}\left(\hat{\theta}_{W},\frac{1}{\lambda}\Sigma_{W}^{-1}\right)(\cdot)\right]\Pi^{(W)}(\mathrm{d}W|\mathcal{D}).$

We set

(3.2)
$\Pi^{(\theta)}_{\star}(\cdot)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\mathbf{N}\left(\theta_{\star},\frac{1}{\lambda}\Sigma_{\star}^{-1}\right)(\cdot).$

We investigate the proximity between $\Pi^{(\theta)}$ and
$\Pi_{\star}^{(\theta)}$ for $n$ large, using the 2-Wasserstein metric that we
denote $\mathsf{W}_{2}$. The following result is established in Section A.5.

######  Theorem 3.1.

Assume H1-H5 and the notations of Theorem 2.2. Then, for all $n$ large enough,
we have

$\mathbb{E}\left[\mathsf{W}_{2}^{2}(\Pi^{(\theta)},\Pi^{(\theta)}_{\star})\right]\leq
c\left(\textsf{r}^{2\kappa}+\lambda
e^{-n\epsilon_{0}^{2}/(2\sigma^{2})}+\lambda
e^{-V_{1}\log(V_{2}\sqrt{n})/C}\right),$

for some constant $c$. In particular,

$\mathbb{E}\left[\left|\int_{\mathbb{R}^{d}}u\Pi^{(\theta)}(\mathrm{d}u|\mathcal{D})-\theta_{\star}\right|^{2}\right]\leq
c\left(\textsf{r}^{2\kappa}+\lambda
e^{-n\epsilon_{0}^{2}/(2\sigma^{2})}+\lambda
e^{-V_{1}\log(V_{2}\sqrt{n})/C}\right).$

For DL functions with sufficiently smooth piecewise polynomial activation
functions, and for $u_{\star}\in\mathcal{C}^{\beta}(\Omega)$, we have seen
above that $\textsf{r}\lesssim n^{-\beta/(m+2\beta)}$, and by Theorem 2.4,
$\kappa=(\beta-\tau)/\beta$. With these we conclude from Theorem 3.1 that the
estimation rate of $\theta_{\star}$ afforded by B-PINN is at least
$n^{-2(\beta-\tau)/(m+2\beta)}$. This estimation rate corresponds to the
optimal plug-in rate. Recall that our pde equation is linear in $\theta$.
Therefore, the bottleneck in estimating $\theta$ is the difficulty of
estimating $u$ and its derivatives up to order $\tau$. As shown above, for a
function $u\in\mathcal{C}^{\beta}(\Omega)$, PINN estimates these derivatives
at the nonparametric minimax rate $n^{-2(\beta-\tau)/(m+2\beta)}$. Hence, the
obtained estimation rate for $\theta_{\star}$. As conjectured above, PINN
likely estimates $u_{\star}$ at a rate faster than the nonparametric minimax
rate. If true, this would translate to a faster estimation rate for
$\theta_{\star}$. More research is needed on this point.

##  4 Numerical illustration

We illustrate the results above using the one-dimensional heat equation of
Example 1. Extensive illustrations of PINN and Bayesian PINN can be found in
the literature (Raissi et al. (2019); Yang et al. (2021); Cuomo et al.
(2022)). We focus here on comparing PINN and an approach that does not
directly exploit the PDE structure of the problem. We consider the heat
equation with $L=\pi$, and $T=1$. Hence $\Omega=(0,1)\times(0,\pi)$. We use
the boundary condition $u(t,0)=u(t,\pi)=0$, and the initial condition
$u(0,x)=\sin(x)$. Given $\theta$ the heat equation has a unique solution
$u_{\theta}(t,x)=\sin(x){\sf exp}(-\theta t)$. We set the true value of
$\theta$ to $\theta_{\star}=0.5$. The observed data $\mathcal{D}$ is generated
as follows: $n$ sensor locations $\\{{\bf s}_{i}=(t_{i},x_{i}),\;1\leq i\leq
n\\}$ are evenly distributed in $\Omega$, yielding measurements

$Y_{i}=u_{\theta_{\star}}({\bf s}_{i})+\sigma\epsilon_{i},\mbox{ where
}\;\epsilon_{i}\sim\mathbf{N}(0,1),$

for some noise parameter $\sigma$ that we control. For the function class
$\mathcal{F}=\\{u_{W}\\}$, $u_{W}:\;\mathbb{R}^{2}\to\mathbb{R}$ is taken to
be a fully connected neural network with depth $4$, width $64$, and tanh
activation function. The resulting posterior distribution is

$\Pi(\mathrm{d}\theta,\mathrm{d}W|\mathcal{D})\propto{\sf
exp}\left(-\frac{1}{2\sigma^{2}}\sum_{i=1}^{n}\left(Y_{i}-u_{W}({\bf
s}_{i})\right)^{2}-\frac{\lambda}{2}\ell_{0}(\theta,W)-\frac{1}{2}\|\theta\|_{2}^{2}\right)\Pi_{0}(\mathrm{d}W)\mathrm{d}\theta,$

where

(4.1)
$\ell_{0}(\theta,W)=|(u_{W})_{t}-\theta(u_{W})_{xx}|_{2}^{2}+|u_{W}(\cdot,0)|_{2}^{2}+|u_{W}(\cdot,\pi)|_{2}^{2}+|u_{W}(0,\cdot)-\sin(\cdot)|_{2}^{2},$

where $|\cdot|_{2}$ denotes the function space $L^{2}$-norm. Observe that
$\lambda=0$ in the above posterior corresponds to the estimation of $W$
without the physics-informed prior (non-PINN estimation). Before presenting
the numerical results, we describe briefly our MCMC sampling method.

###  4.1 Approximation and MCMC sampling

In general, the $L^{2}$-norm in the loss (4.1) is intractable and is typically
replaced by a Monte Carlo or a numerical quadrature approximation. Here we use
Monte Carlo by drawing $N=10,000$ interior points $({\bf
s}_{1}^{(i)},\ldots,{\bf s}^{(i)}_{N})$ uniformly in $\Omega$ to estimate the
$L^{2}$ norm on $\Omega$. For the $L^{2}$ norm along the boundaries we draw
two sets of $B=128$ points $({\bf s}_{1}^{(b_{i})},\ldots,{\bf
s}^{(b_{i})}_{B})$ for $i=1,2$ uniformly on the time and space boundaries
respectively. Hence, we consider the approximate posterior distribution

(4.2) $\widehat{\Pi}(\theta,W|\mathcal{D})\propto\Pi_{0}(W)\\\ {\sf
exp}\left(-\frac{1}{2\sigma^{2}}\sum_{i=1}^{n}\left(y_{i}-u_{W}({\bf
s}_{i})\right)^{2}-\frac{\lambda}{2}\widehat{\ell}_{0}(\theta,W)-\frac{\rho}{2}\|\theta\|_{2}^{2}\right).$

where the loss (4.1) is replaced by

(4.3)
$\widehat{\ell}_{0}(\theta,W)=\frac{1}{N}\sum_{k=1}^{N}\left((u_{W})_{t}({\bf
s}_{k}^{(i)})-\theta\times(u_{W})_{xx}({\bf s}_{k}^{(i)})\right)^{2}\\\
+\frac{1}{B}\sum_{k=1}^{B}\left(u_{W}({\bf
s}^{(b_{1})}_{k},0)\right)^{2}+\frac{1}{B}\sum_{k=1}^{B}\left(u_{W}({\bf
s}^{(b_{1})}_{k},\pi)\right)^{2}+\frac{1}{B}\sum_{k=1}^{B}\left(u_{W}(0,{\bf
s}_{k}^{(b_{2})})-\sin({\bf s}_{k}^{(b_{2})})\right)^{2}.$

Given $W$, and as seen in Section 3, the posterior conditional distribution of
$\theta$ given $W$ is

(4.4)
$\theta|W,\mathcal{D}\sim\mathbf{N}\left(\Sigma_{W}^{-1}\Phi_{W},\frac{1}{\lambda}\Sigma_{W}^{-1}\right)$

where in the particular case of this example

$\Sigma_{W}=\frac{1}{N}\sum_{k=1}^{N}(u_{W})_{xx}(\mathbf{s}_{k}^{(i)})^{2}+\frac{\rho}{\lambda},\;\;\mbox{
and
}\;\;\Phi_{W}=\frac{1}{N}\sum_{k=1}^{N}(u_{W})_{t}(\mathbf{s}_{k}^{(i)})\times(u_{W})_{xx}(\mathbf{s}_{k}^{(i)})$

We sample from (4.2) using the approximate asynchronous sampler of Atchade and
Wang (2023). The algorithm is a data-augmentation Metropolis-with-Gibbs
sampler where the update of the sparsity support $\Lambda$ given $\theta,W$
uses asynchronous sampling, and the update of $W$ given $\theta,\Lambda$ is a
sparse stochastic gradient Langevin dynamics. Then $\theta$ given $W,\Lambda$
is sample from Gaussian
$\mathbf{N}\left(\Sigma_{W_{\Lambda}}^{-1}\Phi_{W_{\Lambda}},\frac{1}{\lambda}\Sigma_{W_{\Lambda}}^{-1}\right)$
using (4.4) where $W_{\Lambda}$ is the component-wise product of $W$ and
$\Lambda$ (a sparse neural network weight). A Pytorch implementation is
available from the GitHub page https://github.com/xliu-522/SA-cSGLD.

Throughout the experiment, we use $\theta^{*}=0.5$, $\lambda=n$ and $\rho=1$.
For each instance of the posterior distribution (4.2) under consideration, we
run the aforementioned MCMC sampler until convergence and keep running another
$200,000$ iterations. We then record every 20th sample, resulting in a total
of $K=10,000$ samples, denoted by $\\{\hat{\theta}^{(k)}\\}_{k=1}^{K}$ . The
mean and standard deviations of the marginal posterior distribution of
$\theta$ are then approximated respectively by

$\mu_{\hat{\theta}}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\frac{1}{K}\sum_{k=1}^{K}\hat{\theta}^{(k)},\;\mbox{
and
}\;\;\sigma_{\hat{\theta}}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\sqrt{\frac{1}{K}\sum_{k=1}^{K}\left(\hat{\theta}^{(k)}-\mu_{\hat{\theta}}\right)^{2}}.$

In this section, assuming that the MCMC sampler has converged, we take the
distribution of the samples $\\{\hat{\theta}^{(k)}\\}_{k=1}^{K}$ to be
$\Pi^{(\theta)}$ and we denote its normal approximation by
$\tilde{\Pi}^{(\theta)}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\mathbf{N}(\mu_{\hat{\theta}},\sigma^{2}_{\hat{\theta}})$.
Since the exact solution $u_{\star}(t,x)=\sin(x){\sf exp}(-0.5t)$ is known, we
can compute the variance of the limiting Gaussian distribution in Bernstein-
von Mises theorem (Theorem 3.1)

$\Sigma_{*}=\int_{0}^{1}\int_{0}^{\pi}\left((u_{\star})_{xx}(t,x)\right)^{2}dxdt+\frac{1}{n}=\frac{\pi(1-\frac{1}{e})}{2}+\frac{1}{n}\approx
0.993+\frac{1}{n}$

Hence the limiting distribution in Theorem 3.1 denoted
$\Pi_{\star}^{(\theta)}$ is

(4.5)
$\Pi_{\star}^{(\theta)}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\mathbf{N}(\theta_{\star},\frac{1}{n}\Sigma_{\star}^{-1}).$

###  4.2 PINN versus non-PINN comparison

In this section, we compare PINN to a two-step approach (denoted by non-PINN)
that does not directly use the PDE structure of $u$. Specifically, in the non-
PINN approach, we fit the same DL model for estimating $u$ but without the
PINN prior. The resulting posterior is given by

$\widehat{\Pi}(W|\mathcal{D})\propto\Pi_{0}(W){\sf
exp}\left(-\frac{1}{2\sigma^{2}}\sum_{i=1}^{n}\left(y_{i}-u_{W}({\bf
s}_{i})\right)^{2}\right).$

We sample from this posterior distribution using the same MCMC algorithm
described above. For each MCMC draw $W$, we then subsequently solve for
$\theta$ by sampling from the Gaussian distribution
$\mathbf{N}\left(\Sigma_{W}^{-1}\Phi_{W},\frac{1}{n}\Sigma_{W}^{-1}\right)$
from (4.4). Hence, in the two-step approach, the PDE information is not used
in the first step for estimating $u$, but is used in the second step to
recover $\theta$ by the linear regression model derived from (1.1). By
contrast, PINN uses the PDE information and estimates $(u,\theta)$ jointly. We
compare the two approaches under varying sample sizes and noise levels, as it
allows us to observe the contribution of the PINN prior $\ell_{0}(\theta,W)$.

Figure 1 provides a comparative analysis of the posterior distribution of
$\theta$ in the two approaches (PINN and non-PINN) at a noise level of 10%
($\sigma\approx 0.025$), across different sample sizes. Both methods exhibit
improved parameter estimates as the sample size $n$ increases, as indicated by
the decreasing bias and interquartile range (IQR). However, training with the
PDE term generally yields more concentrated distributions of $\theta$ with
narrower IQRs, indicating greater stability and reliability in parameter
estimation. This demonstrates that incorporating the PDE term results in more
accurate and robust parameter estimates.

Figure 1: Comparison between boxplots of sampled $\theta$ from training
with/without pde with different sample sizes n at noise level = 10%. True
$\theta^{*}=0.5$ in the dotted line.

To further investigate the role of the PDE term in the recovery of $\theta$,
the root mean square error(rMSE) and the Wasserstein-2 distance ($W_{2}$) are
computed and summarized in Table 1 and Figure 2. The rMSE is defined as
$\sqrt{\frac{1}{K}\sum_{k=1}^{K}(\hat{\theta}^{(k)}-\theta^{*})^{2}}$, whereas
the $W_{2}$ metric is computed using the ot.emd2_1d function from the POT
library by Flamary et al. (2021) with the parameter $p=2$. These metrics
provide a comprehensive comparison of the accuracy (rMSE) and distribution
similarity ($W_{2}$) between sample $\theta$ and the true value $\theta^{*}$.

Figure 2: Comparison between training with/without pde term. The number in each grid is the rMSE between sampled $\theta$ and $\theta^{*}$. The color represents $W_{2}$ distance $W_{2}(\Pi^{(\theta)},\Pi_{*}^{(\theta)})$ Table 1: Comparison of PDE and No-PDE under different noise levels |  | rMSE | $W_{2}(\Pi^{(\theta)},\Pi_{*}^{(\theta)})$
---|---|---|---
Noise Level | n | PDE | No-PDE | PDE | No-PDE
1% | 50 | 0.287 | 0.352 | 0.0379 | 0.0699
1% | 500 | 0.083 | 0.088 | 0.0018 | 0.0024
1% | 1000 | 0.058 | 0.060 | 0.0008 | 0.0010
1% | 5000 | 0.027 | 0.028 | 0.0002 | 0.0003
10% | 500 | 0.103 | 0.140 | 0.0052 | 0.0131
10% | 1000 | 0.080 | 0.116 | 0.0036 | 0.0078
10% | 5000 | 0.031 | 0.040 | 0.0004 | 0.0008
10% | 10000 | 0.023 | 0.031 | 0.0002 | 0.0007
25% | 1000 | 0.094 | 0.377 | 0.0055 | 0.1387
25% | 5000 | 0.042 | 0.104 | 0.0011 | 0.0090
50% | 5000 | 0.129 | 0.474 | 0.0147 | 0.2246

Figure 2 and Table 1 illustrate the performance comparison of models trained
with and without PDE constraints across various noise levels and sample sizes.
As the number of samples ($n$) increases or the noise level decreases, the
values of rMSE and $W_{2}$ distance decrease for both training methods, with
and without the PDE constraint. However, the PDE-constrained models perform
better, demonstrating lower rMSE and $W_{2}$ distances than corresponding
models without PDE constraints. The difference in performance becomes more
noticeable at low sample sizes or high noise levels. This suggests that the
PDE constraint can help the model learn more effectively when there is
insufficient information from the data alone.

###  4.3 Posterior contraction Behavior

From the previous section and Table 1, we observed that the $W_{2}$ distance
between $\Pi^{(\theta)}$ and the limiting distribution
$\Pi_{\star}^{(\theta)}$ given in (4.5) decays significantly as the sample
size $n$ increases, which is consistent with the conclusion of Theorem 3.1. To
further study this contraction behavior, we plot the histogram of
$\Pi^{(\theta)}$ (using samples from the MCMC sampler), its Gaussian
approximation
$\tilde{\Pi}^{(\theta)}=\textbf{N}(\mu_{\hat{\theta}},\sigma^{2}_{\hat{\theta}})$
where $\mu_{\hat{\theta}}$ and $\sigma_{\hat{\theta}}$ are computed from the
MCMC samples, and the limiting distribution $\Pi_{\star}^{(\theta)}$, for
various sample sizes ($n=500$, $1000$, $5000$, $10000$) at a noise level of
10% ($\sigma\approx 0.025$). The perfect match between $\Pi^{(\theta)}$ and
$\tilde{\Pi}^{(\theta)}$ in Figure 3 suggests that the posterior distribution
of $\theta$ is approximately Gaussian. However the persistent discrepancy
between $\Pi^{(\theta)}$ and $\Pi_{\star}^{(\theta)}$ suggests that the
$W_{2}$ convergence to zero between $\Pi^{(\theta)}$ and
$\Pi_{\star}^{(\theta)}$ (as established in Theorem 3.1 and illustrated
above), likely do not hold in total variation. For two probability measures
$P,Q$ with densities $p,q$ respectively, their total variation distance is
$\text{TVD}(P,Q)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\frac{1}{2}\int|p(x)-q(x)|\,dx$.
To further illustrate this point we compute in Table 2 a lower and an upper on
the total variation distance between $\tilde{\Pi}^{(\theta)}$ and
$\Pi_{\star}^{(\theta)}$. We compute these bounds by noting (see Devroye et
al. (2023)) that for any two probability measures $P,Q$ with densities $p,q$
respectively we have

(4.6) $H^{2}(P,Q)\leq\text{TVD}(P,Q)\leq\min\left(1,\sqrt{KL(P\lVert
Q)/2},\sqrt{KL(Q\lVert P)/2}\right),$

where in the univariate Gaussian case where
$P\sim\mathbf{N}(\mu_{1},\sigma_{1}^{2})$ with density $p(x)$ and
$Q\sim\mathbf{N}(\mu_{2},\sigma_{2}^{2})$ with density $q(x)$, the Hellinger
distance $H(P,Q)$ is (see e.g. Pardo (2006))

$H^{2}(P,Q)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\frac{1}{2}\int\left(\sqrt{p(x)}-\sqrt{q(x)}\right)^{2}\,dx=1-\sqrt{\frac{2\sigma_{1}\sigma_{2}}{\sigma_{1}^{2}+\sigma_{2}^{2}}}{\sf
exp}\left(-\frac{(\mu_{1}-\mu_{2})^{2}}{4(\sigma_{1}^{2}+\sigma_{2}^{2})}\right),$

and the $KL$-divergence is (see e.g. Rasmussen and Williams (2005))

$D_{KL}(P\|Q)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\int
p(x)\log\frac{p(x)}{q(x)}\,dx=\log\left(\frac{\sigma_{2}}{\sigma_{1}}\right)+\frac{\sigma_{1}^{2}+(\mu_{1}-\mu_{2})^{2}}{2\sigma_{2}^{2}}-\frac{1}{2}.$

Using these formulas, the third and fourth column of Table 2 shows the left-
hand side and right-hand side of (4.6) respectively. From this table, we see,
for instance, that at $10\%$ noise level, the $W_{2}$ distance decreases with
the sample size, whereas the TVD lower and upper bounds do not. This lack of
TV convergence is because the posterior mean does not converge to
$\theta_{\star}$ fast enough. Indeed, given two univariate normal
distributions $\textbf{N}(\mu_{1,n},\frac{1}{n})$ and
$\textbf{N}(\mu_{2,n},\frac{1}{n})$, as $n\to\infty$, their $W_{2}$ distance
converges to $0$ as soon as $|\mu_{1,n}-\mu_{2,n}|\to 0$, whereas their
Hellinger distance (using the formula above) converges to zero if and only if
$|\mu_{1,n}-\mu_{2,n}|$ converges to zero faster than $1/\sqrt{n}$.

Figure 3: Histograms of sampled $\theta$, denoted $\Pi^{(\theta)}$, in blue with its normal approximation $\tilde{\Pi}^{(\theta)}$ in red line and target Gaussian distribution $\Pi_{*}^{(\theta)}$ in orange with increasing sample size n at noise level 10% Table 2:  $W_{2}$ and upper and lower bound of TVD under different sample sizes and noise levels Noise Level | n | $W_{2}$ | TVD Lower | TVD Upper
---|---|---|---|---
1% | 50 | 0.0377 | 0.149 | 0.482
1% | 500 | 0.0018 | 0.090 | 0.372
1% | 1000 | 0.0008 | 0.084 | 0.360
1% | 5000 | 0.0002 | 0.099 | 0.387
10% | 500 | 0.0051 | 0.176 | 0.515
10% | 1000 | 0.0037 | 0.228 | 0.594
10% | 5000 | 0.0004 | 0.156 | 0.485
10% | 10000 | 0.0002 | 0.173 | 0.512
25% | 1000 | 0.0054 | 0.270 | 0.637
25% | 5000 | 0.0011 | 0.287 | 0.665
50% | 5000 | 0.0148 | 0.616 | 1.000

##  Appendix A Proof of Theorem 2.2

We define

$\Delta_{W}(y,{\bf
s})\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\left(y-u_{W}({\bf
s})\right)^{2},\;\;\Delta_{\star}(y,{\bf
s})\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\left(y-u_{\star}({\bf
s})\right)^{2},\;\;(y,{\bf s})\in\mathbb{R}\times\Omega.$

We recall that $\mathbb{P}_{n}(f)=n^{-1}\sum_{i=1}^{n}f(Y_{i},{\bf s}_{i})$.
It follows from (2.3) that for any set $C\subseteq\mathbb{R}^{q}$, we have

(A.1) $\Pi(C|\mathcal{D})=\frac{\int_{C}{\sf
exp}\left(-\frac{n}{2\sigma^{2}}\mathbb{P}_{n}(\Delta_{W}-\Delta_{\star})+\mathcal{R}_{W}\right)\Pi_{0}(\mathrm{d}W)}{\int_{\mathbb{R}^{q}}{\sf
exp}\left(-\frac{n}{2\sigma^{2}}\mathbb{P}_{n}(\Delta_{W}-\Delta_{\star})+\mathcal{R}_{W}\right)\Pi_{0}(\mathrm{d}W)},$

where

$\mathcal{R}_{W}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\mathcal{R}_{W}^{(1)}+\mathcal{R}_{W}^{(2)},\;\;\mbox{
with
}\;\;\mathcal{R}_{W}^{(1)}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}-\frac{1}{2}\left(\log{\sf
det}(\Sigma_{W})-\log{\sf det}(\Sigma_{\star})\right),\;\;\\\ \;\;\mbox{ and
}\;\;\mathcal{R}_{W}^{(2)}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}-\frac{\lambda}{2}{\cal
J}(u_{W}).$

We break the proof into three parts. First, in Section A.1, we give a lower
bound on the normalizing constant of the posterior distribution as given in
(A.1). Then we show in Section A.2 that the prior distribution $\Pi_{0}$ has a
good inductive bias, as it puts a high probability on $W$ that is sparse. The
third part of the proof in Section A.3 establishes some deviation bounds for
the empirical process of the log-likelihood ratio. Then, we put all the pieces
together in Section A.4.

###  A.1 Lower bound on the normalizing constant

######  Lemma A.1.

Assume H1-H4. Let $\epsilon_{0},s_{0},W_{0}$ be as in H3. Then for all $n$
large enough, with probability at least
$1-2e^{-n\epsilon_{0}^{2}/(2\sigma^{2})}$, we have

$\int_{\mathbb{R}^{q}}{\sf
exp}\left(-\frac{n}{2\sigma^{2}}\mathbb{P}_{n}(\Delta_{W}-\Delta_{\star})+\mathcal{R}_{W}\right)\Pi_{0}(\mathrm{d}W)\geq
C\times{\sf exp}\left(-\frac{Cn\epsilon_{0}^{2}}{\sigma^{2}}\right),$

for some absolute constant $C$.

######  Proof A.2.

By definition of $W_{0}$, $|u_{W_{0}}-u_{\star}|_{\infty}\leq\epsilon_{0}$.
Let $L\stackrel{{\scriptstyle\mathrm{def}}}{{=}}L_{W_{0},1}\geq 1$ be as in
H2. Let $\Lambda_{0}\in\\{0,1\\}^{q}$ denote the sparsity support of $W_{0}$,
and let

(A.2)
$\eta\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\frac{\epsilon_{0}}{Ls_{0}^{1/2}},\;\;\mbox{
and
}\;\;\mathcal{V}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\\{W\in\mathbb{R}^{q}:\;\mathsf{supp}(W)=\Lambda_{0},\;\|W-W_{0}\|_{\infty}\leq\eta\\}.$

Since $\eta s_{0}^{1/2}\leq 1$, for all $W\in\mathcal{V}$, we have

(A.3) $\|W-W_{0}\|_{2}\leq s_{0}^{1/2}\|W-W_{0}\|_{\infty}\leq
s_{0}^{1/2}\eta\leq 1.$

Using (A.3), and the Lipschitz assumption imposed in H2, we deduce that for
all $W\in\mathcal{V}$,

(A.4) $|u_{W}-u_{W_{0}}|_{\infty}\leq L\|W-W_{0}\|_{2}\leq
Ls_{0}^{1/2}\eta\leq\epsilon_{0}.$

From the above display, and appealing to (2.1) and H2-(3) we further deduce
that for all $W\in\mathcal{V}$,

(A.5) $|\mathcal{H}u_{W}-\mathcal{H}u_{W_{0}}|_{2}\leq C_{0}\max_{{\bf
k}:\;|{\bf k}|\leq\tau}\;\left|D^{\bf k}u_{W}-D^{\bf
k}u_{W_{0}}\right|_{\infty}\leq C_{0}c_{1}\epsilon_{0}^{\kappa}.$

In view of H4, we also deduce that for all $W\in\mathcal{V}$,

(A.6) $|u_{W}-u_{\star}|_{\infty}\leq 2\epsilon_{0},\;\;\mbox{ and
}\;\;|\mathcal{H}u_{W}-\mathcal{H}u_{\star}|_{2}\leq
2^{\kappa}C_{4}\epsilon_{0}^{\kappa}.$

#### Step 1

First we show that with probability at least
$1-2e^{-n\epsilon_{0}^{2}/(2\sigma^{2})}$, it holds

(A.7) $\int_{\mathbb{R}^{q}}{\sf
exp}\left(-\frac{n}{2\sigma^{2}}\mathbb{P}_{n}(\Delta_{W}-\Delta_{\star})+\mathcal{R}_{W}\right)\Pi_{0}(\mathrm{d}W)\geq\frac{1}{2}e^{-4n\epsilon_{0}^{2}/\sigma^{2}}\int_{\mathcal{V}}e^{\mathcal{R}_{W}}\Pi_{0}(\mathrm{d}W).$

To establish this, we recall that with $\xi_{i}=(Y_{i}-u_{\star}({\bf
s}_{i}))/\sigma$, we have

$-\frac{n}{2\sigma^{2}}\mathbb{P}_{n}(\Delta_{W}-\Delta_{\star})=-\frac{1}{2\sigma^{2}}\sum_{i=1}^{n}(u_{W}({\bf
s}_{i})-u_{\star}({\bf
s}_{i}))^{2}+\frac{1}{\sigma}\sum_{i=1}^{n}\xi_{i}(u_{W}({\bf
s}_{i})-u_{\star}({\bf s}_{i})).$

Using (A.6), we note that for $W\in\mathcal{V}$,

$\frac{1}{\sigma^{2}}\sum_{i=1}^{n}(u_{W}({\bf s}_{i})-u_{\star}({\bf
s}_{i}))^{2}\leq A,$

where

$A\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\frac{4n\epsilon_{0}^{2}}{\sigma^{2}}.$

By Gaussian tail bounds, for all $W\in\mathcal{V}$,

$\mathbb{P}\left(\frac{1}{\sigma}\sum_{i=1}^{n}\xi_{i}(u_{W}({\bf
s}_{i})-u_{\star}({\bf s}_{i}))>\frac{A}{2}\;|{\bf s}_{1:n}\right)\leq{\sf
exp}\left(-\frac{A^{2}}{8A}\right)\leq e^{-A/8}.$

Therefore, with
$\mathcal{E}_{W}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\\{\mathcal{D}:\;\frac{1}{\sigma}\sum_{i=1}^{n}\xi_{i}(u_{W}({\bf
s}_{i})-u_{\star}({\bf s}_{i}))\leq A/2\\}$,

$\int_{\mathbb{R}^{q}}{\sf
exp}\left(-\frac{n}{2\sigma^{2}}\mathbb{P}_{n}(\Delta_{W}-\Delta_{\star})+\mathcal{R}_{W}\right)\Pi_{0}(\mathrm{d}W)\geq
e^{-A}\int_{\mathcal{V}}\textbf{1}_{\mathcal{E}_{W}}(\mathcal{D})e^{\mathcal{R}_{W}}\Pi_{0}(\mathrm{d}W),$

and with $\mathcal{E}_{W}^{c}$ denoting the complement of $\mathcal{E}_{W}$,
we obtain

$\mathbb{P}\left[\int_{\mathbb{R}^{q}}{\sf
exp}\left(-\frac{n}{2\sigma^{2}}\mathbb{P}_{n}(\Delta_{W}-\Delta_{\star})+\mathcal{R}_{W}\right)\Pi_{0}(\mathrm{d}W)<\frac{e^{-A}}{2}\int_{\mathcal{V}}e^{\mathcal{R}_{W}}\Pi_{0}(\mathrm{d}W)\;|{\bf
s}_{1:n}\right]\\\
\leq\mathbb{P}\left[\int_{\mathcal{V}}\textbf{1}_{\mathcal{E}_{W}}(\mathcal{D})e^{\mathcal{R}_{W}}\Pi_{0}(\mathrm{d}W)<\frac{1}{2}\int_{\mathcal{V}}e^{\mathcal{R}_{W}}\Pi_{0}(\mathrm{d}W)\;|{\bf
x}_{1:n}\right]\\\
=\mathbb{P}\left[\int_{\mathcal{V}}\textbf{1}_{\mathcal{E}_{W}^{c}}(\mathcal{D})e^{\mathcal{R}_{W}}\Pi_{0}(\mathrm{d}W)>\frac{1}{2}\int_{\mathcal{V}}e^{\mathcal{R}_{W}}\Pi_{0}(\mathrm{d}W)\;|\;{\bf
x}_{1:n}\right]\\\
\leq\frac{2}{\int_{\mathcal{V}}e^{\mathcal{R}_{W}}\Pi_{0}(\mathrm{d}W)}\int_{\mathcal{V}}\mathbb{P}\left(\frac{1}{\sigma}\sum_{i=1}^{n}\xi_{i}(u_{W}({\bf
s}_{i})-u_{\star}({\bf s}_{i}))>\frac{A}{2}\;|\;{\bf
s}_{1:n}\right)e^{\mathcal{R}_{W}}\Pi_{0}(\mathrm{d}W)\\\ \leq 2e^{-A/8},$

which is (A.7).

####  Step 2

We now show that

(A.8) $\int_{\mathcal{V}}e^{\mathcal{R}_{W}}\Pi_{0}(\mathrm{d}W)\geq
c\times{\sf exp}\left(-C\frac{n\epsilon_{0}^{2}}{\sigma^{2}}\right).$

Using the definition of $\Sigma_{W}$, we first note that for all $W_{1},W_{2}$

(A.9) $\displaystyle\|\Sigma_{W_{1}}-\Sigma_{W_{2}}\|_{\textsf{op}}$
$\displaystyle\leq$
$\displaystyle\sqrt{2}\left(|\mathcal{H}u_{W_{1}}|_{2}+|\mathcal{H}u_{W_{2}}|_{2}\right)\;|\mathcal{H}u_{W_{1}}-\mathcal{H}u_{W_{2}}|_{2}$
$\displaystyle\leq$
$\displaystyle\sqrt{2}\left(2|\mathcal{H}u_{W_{1}}|_{2}+|\mathcal{H}u_{W_{2}}-\mathcal{H}u_{W_{1}}|_{2}\right)|\mathcal{H}u_{W_{1}}-\mathcal{H}u_{W_{2}}|_{2}.$

We combine this with (A.5) to obtain that for all $W\in\mathcal{V}$,

$\|\Sigma_{W}-\Sigma_{W_{0}}\|_{\textsf{op}}\leq
2C_{0}c_{1}\left(2|\mathcal{H}u_{W_{0}}|_{2}+C_{0}c_{1}\epsilon_{0}^{\kappa}\right)\epsilon_{0}^{\kappa}\leq\frac{C_{2}}{2},$

for all $n$ large enough, since $\epsilon_{0}\to 0$, as $n\to\infty$.
Therefore, Weyl’s inequality, and $\lambda_{\textsf{min}}(\Sigma_{W_{0}})\geq
C_{2}$ imply that for all $W\in\mathcal{V}$, and $n$ large enough

$\lambda_{\textsf{min}}(\Sigma_{W})\geq\frac{C_{2}}{2}.$

As a result of the last display, we can use a first order Taylor expansion of
the $\log{\sf det}$ to conclude that for all $W\in\mathcal{V}$,

$|\mathcal{R}_{W}^{(1)}|=\left|\frac{1}{2}\left(\log{\sf
det}(\Sigma_{W})-\log{\sf
det}(\Sigma_{\star})\right)\right|\leq\frac{2d^{1/2}}{C_{2}}\left\|\Sigma_{W}-\Sigma_{\star}\right\|_{\textsf{F}}\leq\frac{2d}{C_{2}}\|\Sigma_{W}-\Sigma_{\star}\|_{\textsf{op}}\\\
\leq\frac{4d}{C_{2}}\left(2|\mathcal{H}u_{\star}|_{2}+|\mathcal{H}u_{W}-\mathcal{H}u_{\star}|_{2}\right)|\mathcal{H}u_{W}-\mathcal{H}u_{\star}|_{2}\\\
\leq
4\frac{2^{\kappa}C_{4}d}{C_{2}}\left(2|\mathcal{H}u_{\star}|_{2}+2^{\kappa}C_{4}\epsilon_{0}^{\kappa}\right)\epsilon_{0}^{\kappa}\leq
1,$

again, for all $n$ large enough. For $W\in\mathcal{V}$, let $J_{W}$ be the
$L^{2}$ projector on the linear space spanned by the function
$\\{(\mathcal{H}u_{W})_{i},\;1\leq i\leq d\\}$ in
$L^{2}(\Omega,\mathbb{R},\nu)$. Note that, since
$\lambda_{\textsf{min}}(\Sigma_{W})\geq C_{2}/2$, that sub-space is isomorphic
to $\mathbb{R}^{d}$, and by expressing the calculation in $\mathbb{R}^{d}$, it
easily follows that

$\left|\Phi_{W}^{\textsf{T}}\Sigma_{W}^{-1}\Phi_{W}-|J_{W}\bar{f}_{W}|_{2}^{2}\right|\leq\frac{2}{\lambda
C_{2}},$

where we recall that $\bar{f}_{W}=f-\mathcal{H}_{0}u_{W}$. As a result, for
$W\in\mathcal{V}$,

$\left|\Phi_{W}^{\textsf{T}}\Sigma_{W}^{-1}\Phi_{W}-|\bar{f}_{W}|_{2}^{2}\right|\leq\frac{2}{\lambda
C_{2}}+\left||J_{W}\bar{f}_{W}|_{2}^{2}-|\bar{f}_{W}|_{2}^{2}\right|=\frac{2}{\lambda
C_{2}}+\left|J_{W}\bar{f}_{W}-\bar{f}_{W}\right|_{2}^{2},$

where the equality uses the fact that for all $u\in
L^{2}(\Omega,\mathbb{R},\nu)$,
$|u|_{2}^{2}=|J_{W}u|_{2}^{2}+|(J_{W}-\mathbbm{1})u|_{2}^{2}$. By the
definition of the projector as closest element, and since
$\bar{f}_{W}=f-\mathcal{H}_{0}u_{W}=(\mathcal{H}_{0}u_{\star}-\mathcal{H}_{0}u_{W})+\theta_{\star}^{\textsf{T}}\mathcal{H}u_{\star}$,

$\left|J_{W}\bar{f}_{W}-\bar{f}_{W}\right|_{2}^{2}\leq\left|\theta_{\star}^{\textsf{T}}(\mathcal{H}u_{W})-\bar{f}_{W}\right|_{2}^{2}\leq
2\left|\mathcal{H}_{0}u_{\star}-\mathcal{H}_{0}u_{W}\right|_{2}^{2}+2\left|\theta_{\star}^{\textsf{T}}(\mathcal{H}u_{W})-\theta_{\star}^{\textsf{T}}(\mathcal{H}u_{\star})\right|_{2}^{2}\\\
\leq C\left(1+\|\theta_{\star}\|_{2}^{2}\right)\epsilon_{0}^{2\kappa},$

for some absolute constant $C$ that depends only on $\kappa$ and $C_{4}$. We
conclude that for $W\in\mathcal{V}$,

$|\mathcal{R}_{W}^{(2)}|=\frac{\lambda}{2}\left(|f|_{2}^{2}-\Phi_{W}^{\textsf{T}}\Sigma_{W}^{-1}\Phi_{W}\right)\leq\frac{2}{C_{2}}+C(1+\|\theta_{\star}\|_{2}^{2})\lambda\epsilon_{0}^{2\kappa}\leq\frac{1}{C_{2}}+\frac{Cn\epsilon_{0}^{2}}{2\sigma^{2}},$

where in the ast inequality we have used (2.11). Hence, there exists an
absolute constant $C$, such that for all $n$ large enough,

(A.10) $\int_{\mathcal{V}}e^{\mathcal{R}_{W}}\Pi_{0}(\mathrm{d}W)\geq
Ce^{-\frac{Cn\epsilon_{0}^{2}}{2\sigma^{2}}}\Pi_{0}(\mathcal{V}).$

Lemma A.3 below gives the lower bound

$\Pi_{0}(\mathcal{V})\geq\frac{1}{2}{\sf
exp}\left(-s_{0}(\mathsf{u}+2)\log(q)-\frac{s_{0}}{2}(\|W_{0}\|_{\infty}+1)^{2}-s_{0}\log\left(\frac{Ls_{0}^{1/2}}{\epsilon_{0}}\right)\right),$

which together with assumption (2.8) and the inequality in (A.10) yields

$\int_{\mathcal{V}}e^{\mathcal{R}_{W}}\Pi_{0}(\mathrm{d}W)\geq c\times{\sf
exp}\left(-C\left[\frac{n\epsilon_{0}^{2}}{\sigma^{2}}+s_{0}(\mathsf{u}+2)\log(q)\right]\right).$

From the definition of $s_{0}$ in (2.6), we have

$s_{0}(\mathsf{u}+2)\log(q)\leq 6\mathsf{u}(s_{0}-1)\log(q)\leq
6\frac{n\epsilon_{0}^{2}}{\sigma^{2}},$

which then yields (A.8). The lemma follows from (A.7) and (A.8).

We show here that the prior $\Pi_{0}$ has good contraction properties.

######  Lemma A.3.

  1. 1.

Given $s,r\geq 0$, define

$\mathcal{W}(s,r)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\left\\{W\in\mathbb{R}^{q}:\;\|W\|_{0}\leq
s,\;\mbox{ and }\;\|W\|_{\infty}\leq r\right\\}.$

If $r\geq\sqrt{(1+s)(2+\mathsf{u})\log(q)}$, we have

$\Pi_{0}\left(\mathcal{W}(s,r)\right)\geq 1-\frac{4}{q^{\mathsf{u}(1+s)}}.$

  2. 2.

Assume $q\geq\sqrt{2\pi}$. Fix $W_{0}\in\mathbb{R}^{q}$ with sparsity support
$\Lambda_{0}$. Given $r>0$, let

$\mathcal{V}(W_{0},r)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\left\\{W\in\mathbb{R}^{q}:\;\mathsf{supp}(W)=\Lambda_{0},\;\|W-W_{0}\|_{\infty}\leq
r\right\\}.$ $\Pi_{0}\left(\mathcal{V}(W_{0},r)\right)\geq\frac{1}{2}{\sf
exp}\left(-s_{0}(\mathsf{u}+2)\log(q)-\frac{s_{0}}{2}(\|W_{0}\|_{\infty}+r)^{2}+s_{0}\log(r))\right),$

where $s_{0}=\|W_{0}\|_{0}$.

######  Proof A.4.

  1. 1.

Let $\mathcal{W}^{c}$ be a short for $\mathbb{R}^{q}\setminus W(s,r)$. Let
$\bar{\Pi}_{0}$ denote the joint of $(\Lambda,W)$ in the definition of
$\Pi_{0}$ (see Section 1.5). Then

$\Pi_{0}(\mathcal{W}^{c})=\Pi_{0}(\|\Lambda\|_{0}>s)+\sum_{\Lambda:\;\|\Lambda\|_{0}\leq
s}\bar{\Pi}_{0}(\Lambda)\times\bar{\Pi}_{0}(\|\Lambda\odot
W\|_{\infty}>r|\Lambda).$

If $(\Lambda,W)\sim\bar{\Pi}_{0}$, then $\Lambda$ is an ensemble of iid random
variables drawn from the Bernoulli distribution with success probability
$(1+q^{\mathsf{u}+1})^{-1}$. Hence

$\Pi_{0}(\|\Lambda\|_{0}>s)\leq\sum_{j>s}{q\choose
j}\left(\frac{1}{1+q^{\mathsf{u}+1}}\right)^{j}\left(\frac{q^{\mathsf{u}+1}}{1+q^{\mathsf{u}+1}}\right)^{q-j}\\\
\leq\sum_{j>s}{q\choose j}\left(\frac{1}{q^{\mathsf{u}+1}}\right)^{j}\leq
2\left(\frac{1}{q^{\mathsf{u}}}\right)^{s+1},$

where we use ${q\choose j}\leq q^{j}$, and $q^{\mathsf{u}}\geq 2$. Given
$\Lambda_{d}=1$, $W_{d}\sim\mathbf{N}(0,1)$. Therefore,
$\mathbb{P}(|W_{d}|>t)\leq 2e^{-t^{2}/2}$ for all $t\geq 0$. Hence by union
bound, for $\|\Lambda\|_{0}\leq s$, and since
$r\geq\sqrt{(1+s)(2+\mathsf{u})\log(q)}$, we obtain

$\Pi_{0}\left(\|\Lambda\odot W\|_{\infty}>r\;|\;\Lambda\right)\leq
2e^{-r^{2}/2+\log(s)}\leq\frac{2}{q^{\mathsf{u}(1+s)}}.$

We conclude that

(A.11) $\Pi_{0}(\mathcal{W}^{c})\leq\frac{4}{q^{\mathsf{u}(1+s)}}.$

  2. 2.

We write $\mathcal{V}$ as a short for $\mathcal{V}(W_{0},r)$. We have

$\Pi_{0}(\mathcal{V})=\bar{\Pi}_{0}(\Lambda_{0})\bar{\Pi}_{0}(\|\Lambda\odot
W-W_{0}\|_{\infty}\leq r|\Lambda=\Lambda_{0}).$

Since $\log(1-x)\geq-2x$ for all $0\leq x\leq 1/2$, for $q^{\mathsf{u}}\geq
2/\log(2)$, we have

$\bar{\Pi}_{0}(\Lambda_{0})=\left(\frac{1}{1+q^{\mathsf{u}+1}}\right)^{\|\Lambda_{0}\|_{0}}\left(1-\frac{1}{1+q^{\mathsf{u}+1}}\right)^{q-\|\Lambda_{0}\|_{0}}\\\
=\left(\frac{1}{q^{\mathsf{u}+1}}\right)^{\|\Lambda_{0}\|_{0}}{\sf
exp}\left(q\log\left(1-\frac{1}{1+q^{\mathsf{u}+1}}\right)\right)\\\
\geq\left(\frac{1}{q^{\mathsf{u}+1}}\right)^{\|\Lambda_{0}\|_{0}}{\sf
exp}\left(-\frac{2q}{1+q^{\mathsf{u}+1}}\right)\geq\frac{1}{2}\left(\frac{1}{q^{\mathsf{u}+1}}\right)^{\|\Lambda_{0}\|_{0}}=\frac{1}{2}\left(\frac{1}{q^{\mathsf{u}+1}}\right)^{s_{0}}.$

If $U\sim\textbf{N}(0,1)$, and $t\geq 0$, then for all $a$, and $c\geq|a|$,

$P(|U-a|\leq t)\geq P(c\leq U\leq c+t)=\Phi(c+t)-\Phi(c)\geq
e^{-(c+t)^{2}/2}\frac{t}{\sqrt{2\pi}},$

where $\Phi$ is the cdf of the standard normal distribution. We use this
inequality with $c=\|W_{0}\|_{\infty}$, and we the assumption
$q\geq\sqrt{2\pi}$, we deduce that

$\bar{\Pi}_{0}\left(\|\Lambda\odot W-W_{0}\|_{\infty}\leq
r|\Lambda=\Lambda_{0}\right)\geq\frac{1}{q^{s_{0}}}{\sf
exp}\left(-\frac{s_{0}}{2}(\|W_{0}\|_{\infty}+r)^{2}+s_{0}\log(r)\right).$

Hence

$\bar{\Pi}_{0}(\Lambda_{0})\geq\frac{1}{2}\left(\frac{1}{q^{\mathsf{u}+2}}\right)^{s_{0}}{\sf
exp}\left(-\frac{s_{0}}{2}(\|W_{0}\|_{\infty}+r)^{2}+s_{0}\log(r)\right),$

as claimed.

###  A.2 Ignorability of unsuitable weights

Given integer $s\geq 0$,
$\tau_{s}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\sqrt{(2+\mathsf{u})(1+s)\log(q)}$,
and a constant $C$, we define that

$\mathcal{W}_{0}(s,C)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\left\\{W\in\mathbb{R}^{q}:\;\|W\|_{0}\leq
s,\;\mbox{ and }\;\|W\|_{\infty}\leq\tau_{s},\;\mbox{ and }\;{\cal
J}(u_{W})\leq\frac{Cn\epsilon_{0}^{2}}{\lambda\sigma^{2}}\right\\}.$

Our next result shows that the prior $\Pi_{0}$ puts most of its probability
mass on $\mathcal{W}_{0}(s,C)$.

######  Lemma A.5.

Assume H1-H4. Let $s_{0},\epsilon_{0}$, $W_{0}$ be as in H3. We can find an
absolute constant $C$ such that with

(A.12) $s=Cs_{0},$

it holds,

$\mathbb{E}\left[\Pi\left(\mathcal{W}_{0}(s,C)|\mathcal{D}\right)\right]\geq
1-C_{0}e^{-\frac{n\epsilon_{0}^{2}}{2\sigma^{2}}},$

for some absolute constant $C_{0}$.

######  Proof A.6.

For any measurable set $A\subseteq\mathbb{R}^{q}$, we have

$\Pi(A|\mathcal{D})=\frac{\int_{A}{\sf
exp}\left(-\frac{n}{2\sigma^{2}}\mathbb{P}_{n}(\Delta_{W}-\Delta_{\star})+\mathcal{R}_{W}\right)\Pi_{0}(\mathrm{d}W)}{\int_{\mathbb{R}^{q}}{\sf
exp}\left(-\frac{n}{2\sigma^{2}}\mathbb{P}_{n}(\Delta_{W}-\Delta_{\star})+\mathcal{R}_{W}\right)\Pi_{0}(\mathrm{d}W)}.$

By Lemma A.1, we can find an absolute constant $C_{0}$, such that for all $n$
large enough, and with

$\beta\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\frac{n\epsilon_{0}^{2}}{2\sigma^{2}},$
$\mathbb{E}\left[\Pi\left(A|\mathcal{D}\right)\right]\leq
2e^{-n\epsilon_{0}^{2}/(2\sigma^{2})}+\frac{e^{C0\beta}}{C_{0}}\mathbb{E}\left[\int_{A}{\sf
exp}\left(-\frac{n}{2\sigma^{2}}\mathbb{P}_{n}(\Delta_{W}-\Delta_{\star})+\mathcal{R}_{W}\right)\Pi_{0}(\mathrm{d}W)\right].$

By Fubini’s theorem, the expectation on the right hand side of the last
display is

$\int_{A}e^{\mathcal{R}_{W}}\mathbb{E}\left[{\sf
exp}\left(-\frac{n}{2\sigma^{2}}\mathbb{P}_{n}(\Delta_{W}-\Delta_{\star})\right)\right]\Pi_{0}(\mathrm{d}W).$

By conditioning on ${\bf s}_{1:n}$ we see that the expectation inside the last
integral is equal to $1$ for all $W$. Hence,

(A.13) $\mathbb{E}\left[\Pi\left(A|\mathcal{D}\right)\right]\leq
2e^{-n\epsilon_{0}^{2}/(2\sigma^{2})}+\frac{e^{C_{0}\beta}}{C_{0}}\;\int_{A}e^{\mathcal{R}_{W}}\Pi_{0}(\mathrm{d}W).$

By definition, we have $\mathcal{R}_{W}=-\frac{1}{2}\left(\log{\sf
det}(\Sigma_{W})-\log{\sf det}(\Sigma_{\star})\right)-\lambda{\cal
J}(u_{W})/2$. For all $W\in\mathbb{R}^{q}$, we have

(A.14) $-\frac{1}{2}\left(\log{\sf det}(\Sigma_{W})-\log{\sf
det}(\Sigma_{\star})\right)\leq\frac{d}{2}\log\left(n\|\Sigma_{\star}\|_{\mathsf{op}}\right)\leq
C_{1}\beta,$

where the second inequality uses (2.12). We note also that ${\cal
J}(u_{W})\geq 0$. As a result (A.13) becomes

(A.15) $\mathbb{E}\left[\Pi\left(A|\mathcal{D}\right)\right]\leq
2e^{-n\epsilon_{0}^{2}/(2\sigma^{2})}+\frac{e^{C_{0}\beta}}{C_{0}}\;\int_{A}e^{-\frac{\lambda}{2}{\cal
J}(u_{W})}\Pi_{0}(\mathrm{d}W),$

for some possibly different absolute constant $C_{0}$. We apply this with
$A=\mathcal{W}_{0}(Cs_{0},C)^{c}=\\{W\in\mathbb{R}^{q}:\;{\cal
J}(u_{W})>2C\beta/\lambda\\}\cup[\mathcal{W}_{0}(s)^{c}\cap\\{W\in\mathbb{R}^{q}:\;{\cal
J}(u_{W})\leq 2C\beta/\lambda\\}]$. For $C\geq 1+C_{0}$,

$\mathbb{E}\left[\Pi\left(\mathcal{W}(s)^{c}|\mathcal{D}\right)\right]\leq
2e^{-n\epsilon_{0}^{2}/(2\sigma^{2})}+\frac{1}{C_{0}}e^{-n\epsilon_{0}^{2}/(2\sigma^{2})}+\frac{e^{C_{0}\beta}}{C_{0}}\;\int_{\mathcal{W}(s)^{c}}\Pi_{0}(\mathrm{d}W),$

where $\mathcal{W}(s)=\\{W\in\mathbb{R}^{q}:\;\|W\|_{0}\leq s,\;\mbox{ and
}\;\|W\|_{\infty}\leq\tau_{s}\\}$. Lemma A.3-(1) shows that
$\Pi_{0}(\mathcal{W}(s)^{c})\leq 4{\sf exp}(-\mathsf{u}s\log(q))$. As a
result, by taking $s=Cs_{0}$ with $C\geq C_{0}+1$, we have

$\mathsf{u}s\log(q)=C\mathsf{u}s_{0}\log(q)\geq\frac{Cn\epsilon_{0}^{2}}{2\sigma^{2}}\geq\frac{(C_{0}+1)n\epsilon_{0}^{2}}{2\sigma^{2}},$

and we conclude that

$\mathbb{E}\left[\Pi\left(\mathcal{W}_{0}(s)^{c}|\mathcal{D}\right)\right]\leq
2e^{-n\epsilon_{0}^{2}/(2\sigma^{2})}+\frac{1}{C_{0}}e^{-n\epsilon_{0}^{2}/(2\sigma^{2})}\leq
C\times e^{-n\epsilon_{0}^{2}/(2\sigma^{2})},$

for some absolute constant $C$.

###  A.3 Deviation bounds

The last ingredient of the proof is a concentration inequality for the
empirical process of the log-likelihood ratio that we derive next. The proof
is based on Theorem 5.11 of Geer (2000) that we first present. For any random
variable $X$, its $\kappa$-Bernstein norm $\rho_{\kappa}$ is defined as:

(A.16)
$\rho_{\kappa}(X)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\sqrt{2\kappa^{2}\mathbb{E}\left[e^{\frac{|X|}{\kappa}}-1-\frac{|X|}{\kappa}\right]}\,.$

We recall that we say a random variable satisfies Bernstein condition with
parameter $(K,R)$ if:

$\mathbb{E}[|X|^{m}]\leq\frac{1}{2}m!K^{m-2}R^{2}\ \ \forall\ \ m=2,3,\dots$

It is immediate that if a random variable $X$ satisfies Bernstein condition
with parameters $(K,R)$ then $\rho_{2K}(X)\leq\sqrt{2}R$. Let
$X,X_{1:n}\stackrel{{\scriptstyle i.i.d.}}{{\sim}}P$, $\mathbb{P}_{n}$ their
empirical measure. Let $\mathcal{G}$ be a collection of real-valued functions,
where $g(X)$ satisfies the Bernstein condition with parameters $(K,R)$ for all
$g\in\mathcal{G}$. Let $H_{B,2K}(u,\mathcal{G},P)$ denote the bracketing
entropy of $\mathcal{G}$ with respect to the pseudo-metric $\rho_{2K}$.

######  Theorem A.7 (Theorem 5.11 of Geer (2000)).

Suppose that $g(X)$ satisfies the Bernstein condition with parameters $(K,R)$
for all $g\in\mathcal{G}$. Then there exists a universal constant $C$ such
that for any $a,C_{0},C_{1}$ satisfying:

(A.17) $C_{0}\left(\max\left\\{\int_{0}^{R}\sqrt{H_{B,2K}(u,\mathcal{G},P)}\
du,\ R\right\\}\right)\leq
a\leq\min\left\\{\frac{C_{1}\sqrt{n}R^{2}}{2K},8\sqrt{n}R\right\\}$

and $C^{2}\leq C_{0}^{2}/(C_{1}+1)$:

$\mathbb{P}\left(\sup_{g\in\mathcal{G}}\left|\sqrt{n}(\mathbb{P}_{n}-P)g\right|\geq
a\right)\leq Ce^{-\frac{a^{2}}{C^{2}(C_{1}+1)R^{2}}}\,.$

Given $s\geq 0$, we recall the definition

$\mathcal{W}(s)=\left\\{W\in\mathcal{W}:\;\|W\|_{0}\leq
s,\;\;\|W\|_{\infty}\leq\tau_{s}\right\\},\;\mbox{ where
}\;\tau_{s}=\sqrt{(2+\mathsf{u})(1+s)\log(q)}.$

We also recall that

$\Delta_{W}(y,{\bf
s})\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\left(y-u_{W}({\bf
s})\right)^{2},\;\;\Delta_{\star}(y,{\bf
s})\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\left(y-u_{\star}({\bf
s})\right)^{2},\;\;(y,{\bf s})\in\mathbb{R}\times\Omega.$

We use Theorem A.7 to obtain the following.

######  Lemma A.8.

Assume H1-H2. Fix $s\geq 1$, and let $\mathcal{F}_{s}$, $V_{1}=V_{1}(s)$ and
$V_{2}=V_{2}(s)\geq 6b$ be as in H2. There exists an absolute constant
$M_{0}\geq 1$ such that the following holds. For all $M\geq M_{0}$,
$A_{M}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}Mb\sqrt{V_{1}\log(V_{2}\sqrt{n})}$,
and $\delta>0$ such that

$\frac{2\sqrt{2}}{\sqrt{n}}\leq\frac{A_{M}}{\sqrt{n}}\leq\delta\leq 2b,$

we have

$\mathbb{P}\left[\sup_{W\in\mathcal{W}(s):\;|u_{W}-u_{\star}|_{2}\leq\delta}\;\;\left|\sqrt{n}(\mathbb{P}_{n}-\mathbb{P})(\Delta_{W}-\Delta_{\star})\right|>\frac{A_{M}\delta}{2}\left(1+\frac{\sigma}{b}\right)\right]\leq
Ce^{-\frac{A_{M}^{2}}{Cb^{2}}},$

for some absolute constant $C$.

######  Proof A.9.

We set

$Z_{n}(\delta)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\sup_{W\in\mathcal{W}(s):\;|u_{W}-u_{\star}|_{2}\leq\delta}\;\;\left|\sqrt{n}(\mathbb{P}_{n}-\mathbb{P})(\Delta_{W}-\Delta_{\star})\right|\,.$

Recall that, by definition of $Y$, we have $Y_{i}=u_{\star}({\bf
s}_{i})+\sigma\xi_{i}$ where $\xi_{i}\sim\mathbb{N}(0,1)$, $(\xi_{i},{\bf
s}_{i})$ are independent, and $(\xi_{i},{\bf s}_{i})\sim\mathbb{P}$.
Therefore, we can simplify the difference
$(\Delta_{W}-\Delta_{\star})(\xi,{\bf s})$ as:

$(\Delta_{W}-\Delta_{\star})(\xi,{\bf s})=(u_{\star}({\bf s})-u_{W}({\bf
s}))^{2}+2\sigma\xi(u_{\star}({\bf s})-u_{W}({\bf s}))\,,$

and consequently, we can write:

$\displaystyle Z_{n}(\delta)$
$\displaystyle\leq\sup_{W\in\mathcal{W}(s):\;|u_{W}-u_{\star}|_{2}\leq\delta}\;\left|\sqrt{n}(\mathbb{P}_{n}-\mathbb{P})(u_{\star}({\bf
s})-u_{W}({\bf s}))^{2}\right|$
$\displaystyle\qquad\qquad+2\sigma\sup_{W\in\mathcal{W}(s):\;|u_{W}-u_{\star}|_{2}\leq\delta}\;\left|\sqrt{n}(\mathbb{P}_{n}-\mathbb{P})(\xi
u_{\star}({\bf s})-\xi u_{W}({\bf s}))\right|$
$\displaystyle=Z_{n,1}(\delta)+2\sigma Z_{n,2}(\delta)\,,$

where $Z_{n,1}(\delta)$ and $Z_{n,2}(\delta)$ can also be written as

$Z_{n,1}(\delta)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\sup_{h\in\mathcal{H}_{\delta}}\;\left|\frac{1}{\sqrt{n}}\sum_{i=1}^{n}(h({\bf
s}_{i})-\mathbb{P}(h))\right|,\;\;\;\mbox{ and
}\;\;\;Z_{n,2}(\delta)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\sup_{f\in\bar{\mathcal{F}}_{\delta}}\left|\frac{1}{\sqrt{n}}\sum_{i=1}^{n}f(\xi_{i},X_{i})\right|,$

where
$\mathcal{H}_{\delta}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\\{f^{2},\;f\in\mathcal{F}_{\delta}\\}$,
and where
$\mathcal{F}_{\delta}=\\{u_{W}-u_{\star}:W\in\mathcal{W}(s),\;|u_{W}-u_{\star}|_{2}\leq\delta\\}$,
and $\bar{\mathcal{F}}_{\delta}=\\{(\xi,x)\mapsto\xi
f(x),\;f\in\mathcal{F}_{\delta}\\}$.

#### (A) Deviation bound for $Z_{n,1}(\delta)$

By H2-(1), $|h|_{\infty}\leq b^{2}$ for all $h\in\mathcal{H}_{\delta}$.
Therefore for $f\in\mathcal{F}_{\delta}$, and $m\geq 2$, we have:

$\mathbb{E}[|f^{2}(X)|^{m}]\leq(b^{2})^{m-2}\mathbb{E}[|f^{2}(X)|^{2}]\leq(b^{2})^{m-2}(b\delta)^{2}\leq\frac{1}{2}m!(b^{2})^{m-2}(b\delta)^{2}\,.$

Therefore the function class $\mathcal{H}_{\delta}$ satisfies the Bernstein
condition with parameter $(b^{2},b\delta)$. The bracketing number of
$\mathcal{H}_{\delta}$ satisfies

(A.18)
$\mathcal{H}_{B,2b^{2}}(u\sqrt{2},\mathcal{H}_{\delta},\mathbb{P})\leq\log\mathcal{N}(u/2b,\mathcal{F}_{\delta},L_{\infty})\,,$

for all $0<u<b^{2}$. To see this fix some $u>0$ and let $\\{f_{j}\\}_{1\leq
j\leq N}$ be $(u/2b)$-cover of $\mathcal{F}_{\delta}$ with respect to
$L_{\infty}$ norm. Therefore, for any $f\in\mathcal{F}_{\delta}$ there exists
$1\leq j\leq N$ such that $\|f-f_{j}\|_{\infty}\leq u/2b$. Consequently, we
have $\|f^{2}-f_{j}^{2}\|_{\infty}\leq 2b(u/2b)=u$. Now consider the
collection of brackets $\\{(f_{j}^{2}-u/2,f_{j}^{2}+u/2)\\}_{1\leq j\leq N}$.
For this collection we have,
$\rho_{2b^{2}}^{2}\left(f_{j}^{2}+\frac{u}{2}-f_{j}^{2}+\frac{u}{2}\right)=\rho_{2b^{2}}^{2}(u)\leq(\sqrt{2}u)^{2}$,
for $u\leq b^{2}$. Hence we have established (A.18).

We can then apply Theorem A.7 with $a=A_{M}\delta/2$, $R=b\delta$, $K=b^{2}$,
$C_{1}=1$, $C_{0}^{2}=2C$, where $C$ is the absolute constant in Theorem A.7.
The condition $a\geq C_{0}R$ is satisfied with $M\geq M_{0}\geq 2C_{0}$. For
$\delta\leq 2b$, using H2-(4),

$\displaystyle\int_{0}^{b\delta}\sqrt{\mathcal{H}_{B,2b^{2}}(u,\mathcal{H}_{\delta},\mathbb{P})}\mathrm{d}u$
$\displaystyle=\sqrt{2}\int_{0}^{b\delta/\sqrt{2}}\sqrt{\mathcal{H}_{B,2b^{2}}(u\sqrt{2},\mathcal{H}_{\delta},\mathbb{P})}\
du$
$\displaystyle\leq\sqrt{2}\int_{0}^{b\delta/\sqrt{2}}\sqrt{\log{\mathcal{N}(u/2b,\mathcal{F}_{\delta},L_{\infty})}}\
du$ $\displaystyle\leq
2b\sqrt{2V_{1}}\int_{0}^{\delta/2\sqrt{2}}\sqrt{\log{\left(\frac{V_{2}}{u}\right)}}\
du$ $\displaystyle\leq
2b\delta\sqrt{V_{1}}\sqrt{\log{\left(\frac{2\sqrt{2}V_{2}}{\delta}\right)}}\leq
2b\delta\sqrt{V_{1}\log{\left(V_{2}\sqrt{n}\right)}},$

where the inequality before last follows from the following fact:

$\int_{0}^{c}\sqrt{\log{\left(\frac{C}{u}\right)}}\ du\leq
2c\sqrt{\log{\left(\frac{C}{c}\right)}}\ \ \text{as soon as }\log{(C/c)}\geq
2\,,$

which in our case is implied by $V_{2}\geq 6b$. The last inequality uses
$\delta\geq 2\sqrt{2}/\sqrt{n}$.

Therefore the condition $a\geq
C_{0}\int_{0}^{b\delta}\sqrt{\mathcal{H}_{B,2b^{2}}(u,\mathcal{H}_{\delta},\mathbb{P})}\mathrm{d}u$
is satisfied with $M\geq M_{0}\geq 4C_{0}$. The condition $a\leq
C_{1}\sqrt{n}R^{2}/(2K)$ boils down to $A_{M}/\sqrt{n}\leq\delta$, whereas the
condition $a\leq 8R\sqrt{n}$ boils down to $A_{M}/\sqrt{n}\leq 16b$ which
holds by assumption. Hence by Theorem A.7,

$\mathbb{P}\left(Z_{n,1}(\delta)>\frac{A_{M}\delta}{2}\right)\leq
Ce^{-\frac{A_{M}^{2}}{Cb^{2}}},$

for some absolute constant $C$.

#### (B) Deviation bound for $Z_{n,2}(\delta)$

The argument is similar. Note that, for any
$u_{W}-u_{\star}\in\mathcal{F}_{\delta}$, $|u_{W}-u_{\star}|_{2}\leq\delta$,
and since the function $u_{W}-u_{\star}$ are bounded by $b$ as assumed in H2,
we have:

$\displaystyle\mathbb{E}[|\xi f(X)|^{m}]$
$\displaystyle=\mathbb{E}[|\xi|^{m}]\mathbb{E}[|u_{W}(X)-u_{\star}(X)|^{m}]$
$\displaystyle\leq\mathbb{E}[|\xi|^{m}]b^{m-2}\delta^{2}$
$\displaystyle\leq\frac{1}{2}m!b^{m-2}\delta^{2}\,,$

where the last inequality follows from the centered absolute moment of
Gaussian random variables. Therefore, all $f\in\bar{\mathcal{F}}_{\delta}$
satisfies the Bernstein condition with parameters $(b,\delta)$.

The bracketing entropy satisfies,

(A.19)
$\mathcal{H}_{B,2b}(u\sqrt{2},\bar{\mathcal{F}}_{\delta},\mathbb{P})\leq\log{\mathcal{N}(u/2,\mathcal{F}_{\delta},L_{\infty})}\,\;\;0<u\leq\delta\,.$

To establish (A.19), fix $u\in(0,b]$, and let $\\{f_{1},\dots,f_{N}\\}$ be a
$u/2$ cover of $\mathcal{F}_{\delta}$ with respect to $L_{\infty}$ norm, i.e.

$\sup_{f\in\mathcal{F}_{\delta}}\min_{1\leq j\leq
N}\|f-f_{j}\|_{\infty}\leq\frac{u}{2}\,.$

Now consider the brackets $\\{(\xi,x)\mapsto(\xi f_{j}(x)-u|\xi|/2,\xi
f_{j}(x)+u|\xi|/2),\;1\leq j\leq N\\}$. For any $f$, there exists $f_{j}$ such
that $\|f-f_{j}\|_{\infty}\leq u/2$ due to the property of the covering set,
and it follows that

$\xi f_{j}(x)-\frac{u|\xi|}{2}\leq\xi f(x)\leq\xi f_{j}(x)+\frac{u|\xi|}{2}\ \
\forall\ \ x\,.$

Furthermore, for any $1\leq j\leq N$, since $\mathbb{E}(|u\xi|^{m})\leq
u^{m}m!/2\leq b^{m-2}u^{2}m!/2$,

$\rho_{2b}^{2}\left(\xi f_{j}(x)+u|\xi|/2-\xi
f_{j}(x)+u|\xi|/2\right)=\rho_{2b}^{2}\left(u|\xi|\right)\leq(\sqrt{2}u)^{2},$

which implies (A.19).

Hence we can apply Theorem A.7 with $a=A_{M}\delta/2b$, $R=\delta$, $K=b$,
$C_{1}=1$, $C_{0}^{2}=2C$, where $C$ is the absolute constant in Theorem A.7.
All the conditions of Theorem A.7 can be checked as we did in bounding
$Z_{n,1}(\delta)$. We conclude that

$\mathbb{P}\left(Z_{n,2}(\delta)>\frac{A_{M}\delta}{2b}\right)\leq
Ce^{-\frac{A_{M}^{2}}{Cb^{2}}},$

for some absolute constant $C$. Combining the two bounds, we get

$\mathbb{P}\left(Z_{n}(\delta)>\frac{A_{M}\delta}{2}\left(1+\frac{\sigma}{b}\right)\right)\leq
Ce^{-\frac{A_{M}^{2}}{Cb^{2}}},$

for some absolute constant $C$, which implies the stated result.

###  A.4 Finishing the proof

######  Proof A.10.

Let $s_{0},\epsilon_{0}$, and $W_{0}$ be as in H3. Let $s=Cs_{0}$ as in
(A.12). We set $\tau_{s}=\sqrt{(2+\mathsf{u})(1+s)\log(q)}$,

$\mathcal{W}_{0}(s)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\left\\{W\in\mathbb{R}^{q}:\;\|W\|_{0}\leq
s,\;\|W\|_{\infty}\leq\tau_{s},\;{\cal
J}(u_{W})\leq\frac{Cn\epsilon_{0}^{2}}{\lambda\sigma^{2}}\right\\}.$

We set
$\textsf{r}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}2M_{0}(b+\sigma)\sqrt{V_{1}\log(V_{2}\sqrt{n})/n}$,
where $M_{0}$ is as in Lemma A.8. For $j\geq 1$, we set
$A_{j}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}bM_{0}j\sqrt{V_{1}\log(V_{2}\sqrt{n})}$,
$\textsf{r}_{j}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}j\textsf{r}$. With
$C$ as in Lemma A.1, we set

$\beta\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\frac{C}{\sigma^{2}}n\epsilon_{0}^{2},$

and define

$\mathcal{E}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\left\\{\mathcal{D}:\;\int_{\mathbb{R}^{q}}{\sf
exp}\left(-\frac{n}{2\sigma^{2}}\mathbb{P}_{n}(\Delta_{W}-\Delta_{\star})+\mathcal{R}_{W}\right)\Pi_{0}(\mathrm{d}W)\leq
Ce^{-\beta},\;\;\mbox{ or }\;\right.\\\
\left.\sup_{W\in\mathcal{W}(s):\;\textsf{r}_{j}<|u_{W}-u_{\star}|_{2}\leq\textsf{r}_{j+1}}\;\left|\sqrt{n}(\mathbb{P}_{n}-\mathbb{P})(\Delta_{W}-\Delta_{\star})\right|>\frac{A_{j+1}\textsf{r}_{j+1}}{2}\left(1+\frac{\sigma}{b}\right)\;\;\mbox{
for some }\;j\geq 1\right\\}.$

By Lemma A.1 and Lemma A.8, and for all $n$ large enough

(A.20) $\mathbb{P}(\mathcal{D}\in\mathcal{E})\leq
2e^{-n\epsilon_{0}^{2}/(2\sigma^{2})}+C\sum_{j\geq
1}e^{-A_{j}^{2}/(Cb^{2})}\leq
2e^{-n\epsilon_{0}^{2}/(2\sigma^{2})}+Ce^{-V_{1}\log(V_{2}\sqrt{n})/C},$

for some absolute constant $C$.

We set
$B\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\\{W\in\mathbb{R}^{q}:\;|u_{W}-u_{\star}|_{2}>\textsf{r}\\}$.
We can write $B\cap\mathcal{W}(s)$ as $\cup_{j\geq 1}B_{j}$, where

$B_{j}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\left\\{W\in\mathcal{W}(s):\;Mj\textsf{r}<|u_{W}-u_{\star}|_{2}\leq
M(j+1)\textsf{r}\right\\}.$

We should point out that the union $\cup_{j\geq 1}B_{j}$ is over a finite
number of terms since $B_{j}=\emptyset$ for $Mj\textsf{r}\geq b$. Since
$\Pi(B|\mathcal{D})\leq\Pi(\mathcal{W}_{0}(s)^{c}|\mathcal{D})+\textbf{1}_{\mathcal{E}}(\mathcal{D})+\textbf{1}_{\mathcal{E}^{c}}(\mathcal{D})\Pi(B\cap\mathcal{W}_{0}(s)|\mathcal{D})$,
taking expectation on both sides, and using Lemma A.5, the definition of
$\mathcal{E}$ and (A.20), yields for all $n$ large enough,

$\mathbb{E}\left[\Pi(B|\mathcal{D})\right]\leq
Ce^{-n\epsilon_{0}^{2}/(2\sigma^{2})}+Ce^{-V_{1}\log(V_{2}\sqrt{n})/C}\\\
+\frac{e^{\beta}}{C}\sum_{j\geq
1}\mathbb{E}\left[\textbf{1}_{\mathcal{E}^{c}}(\mathcal{D})\;\int_{B_{j}}{\sf
exp}\left(-\frac{n}{2\sigma^{2}}\mathbb{P}_{n}(\Delta_{W}-\Delta_{\star})+\mathcal{R}_{W}\right)\Pi_{0}(\mathrm{d}W)\right],$

for some absolute constant $C$. We have

$\frac{n}{2\sigma^{2}}\mathbb{P}_{n}(\Delta_{W}-\Delta_{\star})+\mathcal{R}_{W}=-\frac{n}{2\sigma^{2}}\varrho(u_{W},u_{\star})-\frac{n}{2\sigma^{2}}(\mathbb{P}_{n}-\mathbb{P})(\Delta_{W}-\Delta_{\star})+\mathcal{R}_{W}^{\prime},$

where

$\mathcal{R}_{W}^{\prime}=-\frac{1}{2}\left(\log{\sf
det}\Sigma_{\star}-\log{\sf det}\Sigma_{W}\right),$

and satisfies, as we show in the proof of Lemma A.5

$\left|\mathcal{R}_{W}^{\prime}\right|\leq\frac{d}{2}\log\left(\frac{n\|\Sigma_{\star}\|_{\textsf{op}}}{\sigma^{2}}\right)\leq
C_{0}\log(q),$

where the second inequality is our assumption (2.12). Hence, we can find an
absolute constant $C$ such that for all $n$ large enough

(A.21) $\mathbb{E}\left[\Pi(B|\mathcal{D})\right]\leq
Ce^{-n\epsilon_{0}^{2}/(2\sigma^{2})}+Ce^{-V_{1}\log(V_{2}\sqrt{n})/C}\\\
+\frac{e^{\beta}}{C}\sum_{j\geq
1}\mathbb{E}\left[\textbf{1}_{\mathcal{E}^{c}}(\mathcal{D})\;\int_{B_{j}}{\sf
exp}\left(-\frac{n}{2\sigma^{2}}\varrho(u_{W},u_{\star})-\frac{n}{2\sigma^{2}}(\mathbb{P}_{n}-\mathbb{P})(\Delta_{W}-\Delta_{\star})\right)\Pi_{0}(\mathrm{d}W)\right].$

For $W\in B_{j}$, and $\mathcal{D}\notin\mathcal{E}$, the expression inside
the exponential in the last display is bounded from above by

$-\frac{n\textsf{r}_{j}^{2}}{2\sigma^{2}}+\frac{\sqrt{n}A_{j+1}\textsf{r}_{j+1}}{2\sigma^{2}}\left(1+\frac{\sigma}{b}\right)\leq-\frac{n\textsf{r}_{j}^{2}}{4\sigma^{2}},$

with the choice of r. Since

$\sum_{j\geq 1}e^{-\frac{n\textsf{r}_{j}^{2}}{4\sigma^{2}}}\leq
Ce^{-\frac{n\textsf{r}^{2}}{4\sigma^{2}}},$

for some absolute constant $C$, we conclude that for all $n$ large enough,

$\mathbb{E}\left[\Pi(B|\mathcal{D})\right]\leq
C\left(e^{-n\epsilon_{0}^{2}/(2\sigma^{2})}+e^{-V_{1}\log(V_{2}\sqrt{n})/C}+e^{-n\textsf{r}^{2}/(4\sigma^{2})}\right).$

Hence the theorem.

###  A.5 Proof of Theorem 3.1

######  Proof A.11.

First, a simple coupling argument shows that

(A.22)
$\mathsf{W}_{2}^{2}(\Pi^{(\theta)},\Pi^{(\theta)}_{\star})\leq\int_{\mathbb{R}^{q}}\mathsf{W}_{2}^{2}\left(\textbf{N}(\hat{\theta}_{W},\frac{1}{\lambda}\Sigma_{W}^{-1}),\;\textbf{N}(\theta_{\star},\frac{1}{\lambda}\Sigma_{\star}^{-1})\right)\Pi^{(W)}(\mathrm{d}W|\mathcal{D}).$

We recall (see e.g. Bhatia et al. (2019)) that for symmetric and positive
definite matrices $\Lambda_{1},\Lambda_{2}$,

$\mathsf{W}_{2}^{2}\left(\textbf{N}(\mu_{1},\Lambda_{1}),\textbf{N}(\mu_{2},\Lambda_{2})\right)=\|\mu_{2}-\mu_{1}\|_{2}^{2}+\textsf{Tr}\left(\Lambda_{1}+\Lambda_{2}-2\left(\Lambda_{1}^{1/2}\Lambda_{2}\Lambda_{1}^{1/2}\right)^{1/2}\right)\\\
\leq\|\mu_{2}-\mu_{1}\|_{2}^{2}+\|\Lambda_{1}^{1/2}-\Lambda_{2}^{1/2}\|_{\textsf{F}}^{2},$

where the inequality follows from Theorem 1 of (Bhatia et al. (2019)). By the
Hemmen-Ando inequality (Ando and van Hemmen (1980) Proposition 2.1),

$\|\Lambda_{1}^{1/2}-\Lambda_{2}^{1/2}\|_{\textsf{F}}\leq\frac{\|\Lambda_{1}-\Lambda_{2}\|_{\textsf{F}}}{\sqrt{\lambda_{\textsf{min}}(\Lambda_{1})}+\sqrt{\lambda_{\textsf{min}}(\Lambda_{2})}}.$

Hence

(A.23)
$\mathsf{W}_{2}^{2}\left(\textbf{N}(\mu_{1},\Lambda_{1}),\textbf{N}(\mu_{2},\Lambda_{2})\right)\leq\|\mu_{2}-\mu_{1}\|_{2}^{2}+\frac{\|\Lambda_{1}-\Lambda_{2}\|_{\textsf{F}}^{2}}{\left(\sqrt{\lambda_{\textsf{min}}(\Lambda_{1})}+\sqrt{\lambda_{\textsf{min}}(\Lambda_{2})}\right)^{2}}.$

We apply this bound, to conclude that there exists a constant $c$ (that we can
take as $d/\lambda_{\textsf{max}}(\Sigma_{\star})$ such that for all
$W\in\mathbb{R}^{q}$,

$\mathsf{W}_{2}^{2}\left(\textbf{N}(\hat{\theta}_{W},\frac{1}{\lambda}\Sigma_{W}^{-1}),\;\textbf{N}(\theta_{\star},\frac{1}{\lambda}\Sigma_{\star}^{-1})\right)\leq\|\hat{\theta}_{W}-\theta_{\star}\|_{2}^{2}+\frac{c}{\lambda}\|\Sigma_{W}^{-1}-\Sigma_{\star}^{-1}\|^{2}_{\textsf{op}}.$

We consider the integrand of (A.22) under two scenarios.

#### Case 1: $W$ is such that $|u_{W}-u_{\star}|_{2}\leq M\textsf{r}$

By (A.9), and for $W$ such that $|u_{W}-u_{\star}|_{2}\leq M\textsf{r}$, for
some absolute onstant $C$ we have

$\|\Sigma_{W}-\Sigma_{\star}\|_{\textsf{op}}\leq
C\left(|\mathcal{H}u_{\star}|_{2}+(M\textsf{r})^{\kappa}\right)(M\textsf{r})^{\kappa},$

for all $n$ large enough. So, since
$\lambda_{\textsf{min}}(\Sigma_{\star})>C_{2}$, we conclude that for all such
$W$ under consideration, and for all $n$ large enough,
$\lambda_{\textsf{min}}(\Sigma_{W})>C_{2}/2$.

For all $W$,

$(\hat{\theta}_{W}-\theta_{\star})^{\textsf{T}}\mathcal{H}_{1}u_{W}=\hat{\theta}_{W}^{\textsf{T}}\mathcal{H}_{1}u_{W}-\theta_{\star}^{\textsf{T}}\mathcal{H}_{1}u_{\star}-\theta_{\star}^{\textsf{T}}(\mathcal{H}_{1}u_{W}-\mathcal{H}_{1}u_{\star})\\\
=\hat{\theta}_{W}^{\textsf{T}}\mathcal{H}_{1}u_{W}-\mathcal{H}_{0}u_{W}-f+(\mathcal{H}_{0}u_{W}-\mathcal{H}_{0}u_{\star})-\theta_{\star}^{\textsf{T}}(\mathcal{H}_{1}u_{W}-\mathcal{H}_{1}u_{\star}).$

Therefore

$|(\hat{\theta}_{W}-\theta_{\star})^{\textsf{T}}\mathcal{H}_{1}u_{W}|_{2}^{2}\leq
2|f-\mathcal{H}_{0}u_{W}-\hat{\theta}_{W}^{\textsf{T}}\mathcal{H}_{1}u_{W}|_{2}^{2}+c\times|\mathcal{H}u_{W}-\mathcal{H}u_{\star}|_{2}^{2}.$

Since $\hat{\theta}_{W}$ minimizes
$\theta\mapsto|f-\mathcal{H}_{0}u_{W}-\theta^{\textsf{T}}\mathcal{H}_{1}u_{W}|_{2}+\|\theta\|_{2}^{2}/\lambda$,
we have

$|f-\mathcal{H}_{0}u_{W}-\hat{\theta}_{W}^{\textsf{T}}\mathcal{H}_{1}u_{W}|_{2}^{2}\leq|f-\mathcal{H}_{0}u_{W}-\theta_{\star}^{\textsf{T}}\mathcal{H}_{1}u_{W}|_{2}^{2}+\frac{\|\theta_{\star}\|_{2}^{2}}{\lambda}\\\
\leq
c\times|\mathcal{H}u_{W}-\mathcal{H}u_{\star}|_{2}^{2}+\frac{\|\theta_{\star}\|_{2}^{2}}{\lambda}.$

Using the above, and H4, we conclude that we can find a constant $c$ such that
for $W$ such that $|u_{W}-u_{\star}|_{2}\leq M\textsf{r}$, it holds

$|(\hat{\theta}_{W}-\theta_{\star})^{\textsf{T}}\mathcal{H}_{1}u_{W}|_{2}^{2}\leq
c\textsf{r}^{2\kappa}.$

However,
$|(\hat{\theta}_{W}-\theta_{\star})^{\textsf{T}}\mathcal{H}_{1}u_{W}|_{2}^{2}=(\hat{\theta}_{W}-\theta_{\star})^{\textsf{T}}(\Sigma_{W}-(1/\lambda)\textbf{1}_{d})(\hat{\theta}_{W}-\theta_{\star})$.
And since $\lambda_{\textsf{min}}(\Sigma_{W})>C/2$ as seen above, and since
$\lambda\to\infty$, as $n\to\infty$, for all $n$ large enough we have

$|(\hat{\theta}_{W}-\theta_{\star})^{\textsf{T}}\mathcal{H}_{1}u_{W}|_{2}^{2}\geq
C\|\hat{\theta}_{W}-\theta_{\star}\|_{2}^{2},$

for some absolute constant $C$. In conclusion, for $W$ such that
$|u_{W}-u_{\star}|_{2}\leq M\textsf{r}$, and for all $n$ large enough,

$\|\hat{\theta}_{W}-\theta_{\star}\|_{2}^{2}\leq c\textsf{r}^{2\kappa},$

and

$\mathsf{W}_{2}^{2}\left(\textbf{N}(\hat{\theta}_{W},\frac{1}{\lambda}\Sigma_{W}^{-1}),\;\textbf{N}(\theta_{\star},\frac{1}{\lambda}\Sigma_{\star}^{-1})\right)\leq
c\textsf{r}^{2\kappa}.$

#### Case 2: $W$ is such that $|u_{W}-u_{\star}|_{2}>M\textsf{r}_{n}$

Since the smallest eigenvalue of $\Sigma_{W}$ is at least $1/\lambda$, there
exists a constant $c$ such that for all $W\in\mathbb{R}^{q}$, we have

$\|\Sigma_{W}^{-1}-\Sigma_{\star}^{-1}\|_{\textsf{op}}\leq c\lambda.$

Similarly,

$\|\hat{\theta}_{W}-\theta_{\star}\|_{2}\leq\|\hat{\theta}_{W}\|_{2}+\|\theta_{\star}\|_{2},$

and

$\|\hat{\theta}_{W}\|_{2}^{2}=\hat{\theta}_{W}^{\textsf{T}}\Sigma_{W}^{1/2}\Sigma_{W}^{-1}\Sigma_{W}^{1/2}\hat{\theta}_{W}\leq\lambda\hat{\theta}_{W}^{\textsf{T}}\Sigma_{W}\hat{\theta}_{W}=\lambda\Phi_{W}^{\textsf{T}}\Sigma_{W}^{-1}\Phi_{W}\leq\lambda|\bar{f}_{W}|_{2}^{2}\leq
2\lambda(|f|_{2}^{2}+|\mathcal{H}_{0}u_{W}|_{2}^{2}).$

We deduce that for all $W$ such that $|u_{W}-u_{\star}|_{2}>M\textsf{r}$, we
have

$\mathsf{W}_{2}^{2}\left(\textbf{N}(\hat{\theta}_{W},\frac{\sigma^{2}}{\lambda}\Sigma_{W}^{-1}),\;\textbf{N}(\theta_{\star},\frac{\sigma^{2}}{\lambda}\Sigma_{\star}^{-1})\right)\leq
c\lambda\left(1+|\mathcal{H}_{0}u_{W}|_{2}^{2}\right),$

for some constant $c$. Given the behaviors of the $\mathsf{W}_{2}$ distance
between
$\textbf{N}(\hat{\theta}_{W},\frac{\mathbbm{1}}{\lambda}\Sigma_{W}^{-1})$ and
$\textbf{N}(\theta_{\star},\frac{1}{\lambda}\Sigma_{\star}^{-1})$ obtained in
the cases above, we return to (A.22) to write

$\mathsf{W}_{2}^{2}(\Pi^{(\theta)},\Pi^{(\theta)}_{\star})\leq
c\textsf{r}^{2\kappa}+c\lambda\int_{\mathbb{R}^{q}}\left(1+|\mathcal{H}_{0}u_{W}|_{2}^{2}\right)\textbf{1}_{\\{|u_{W}-u_{\star}|_{2}>M\textsf{r}_{n}\\}}\Pi^{(W)}(\mathrm{d}W|\mathcal{D}).$

Taking the expectation on both sides and with a similar argument as in the
proof of Theorem 2.2 yields

$\mathbb{E}\left[\mathsf{W}_{2}^{2}(\Pi^{(\theta)},\Pi^{(\theta)}_{\star})\right]\leq
c\left(\textsf{r}^{2\kappa}+\lambda
e^{-n\epsilon_{0}^{2}/(2\sigma^{2})}+\lambda
e^{-V_{1}\log(V_{2}\sqrt{n})/C}\right),$

which is the stated bound.

###  A.6 Proof of Theorem 2.4

######  Proof A.12.

The result is based on a generalization of Markov’s polynomial inequality due
to Harris (2008) Theorem 1. The result shows that differential operators are
bounded operators when restricted to polynomials, and an upper on their norms
is given.

Let $X,Y$ be a real Banach space with norms $\|\cdot\|_{X}$, and
$\|\cdot\|_{Y}$ respectively. A function $P:X\to Y$ is a homogeneous
polynomial of degree $d\geq 0$ if $P(x)=L(x,\cdots,x)$, where $L$ is a linear
symmetric map from $\underbrace{X\times\cdots\times X}_{d\mbox{ times }}$ to
$Y$ (for $d=0$, these are constant functions on $X$). A function $P:X\to Y$ is
a polynomial of degree $d$ is $P=\sum_{j=0}^{d}P_{j}$, where $P_{j}$ is a
homogeneous polynomial of degree $j$. We let $\wp_{d}(X,Y)$ be the set of all
polynomials $P:\;X\to Y$ with degree at most $d$. When $Y=\mathbb{R}$, we
write $\wp_{d}(X)$. If $f:X\to Y$ has Frechet derivatives to order $k$ we
write $\nabla^{(k)}f(x)$ to denote its $k$-th order Frechet derivative, and
for $z\in X$, $\nabla^{(k)}f(x)\cdot
z^{k}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\nabla^{(k)}f(x)(z,\ldots,z)$.
Furthermore, we set

$\|\nabla^{(k)}f(x)\|\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\sup_{z\in
X:\;\|z\|_{X}\leq 1}\|\nabla^{(k)}f(x)(z,\ldots,z)\|_{Y}.$

Let $T_{d}(t)=\cos(d\arccos(t))$, $t\in[-1,1]$ denote the Chebyshev polynomial
of degree $d$, and $T^{(j)}_{d}(t)$ its $j$-th order derivative. The following
lemma is due to Harris (2008) Theorem 1.

######  Lemma A.13.

For $P\in\wp_{d}(\mathbb{R}^{m})$, and $k\geq 1$, let $P^{(k)}$ denote the
$k$-th order derivative of $P$. We have

$\sup_{x:\;\|x\|_{2}\leq 1}\|P^{(k)}(x)\|\leq
T_{d}^{(k)}(1)\;\times\sup_{x:\;\|x\|_{2}\leq 1}|P(x)|,$

where
$\|P^{(k)}(x)\|\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\sup_{y:\|y\|_{2}\leq
1}|P^{(k)}(x)(y,\ldots,y)|$.

We turn to the proof of Theorem 2.4. Let $\alpha=\bar{\alpha}$ as defined in
the statement of the theorem, and fix $x_{0}\in\mathsf{int}_{\alpha}(\Omega)$.
The function $u$ has derivatives to the order $\lfloor\beta\rfloor$ at
$x_{0}$. By Taylor approximation we have, for all $z\in\mathbb{R}^{m}$, with
$\|z\|_{2}\leq 1$, setting
$r\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\lfloor\beta\rfloor-1$,

$u(x_{0}+\alpha z)=\\\ u(x_{0})+\alpha\nabla u(x_{0})\cdot
z+\cdots+\frac{\alpha^{r}}{r!}\nabla^{(r)}u(x_{0})\cdot
z^{r}+\frac{\alpha^{\lfloor\beta\rfloor}}{r!}\int_{0}^{1}(1-t)^{r}\nabla^{(\lfloor\beta\rfloor)}u\left(x_{0}+\alpha
tz\right)\cdot z^{\lfloor\beta\rfloor}\mathrm{d}t\\\ =u(x_{0})+\alpha\nabla
u(x_{0})\cdot
z+\cdots+\frac{\alpha^{\lfloor\beta\rfloor}}{\lfloor\beta\rfloor!}\nabla^{(\lfloor\beta\rfloor)}u(x_{0})\cdot
z^{\lfloor\beta\rfloor}+R_{u}(z),$

where, using the fact that $u\in\mathcal{C}^{\beta}(\Omega,M)$, the remainder
satisfies

$|R_{u}(z)|=\frac{\alpha^{\lfloor\beta\rfloor}}{r!}\left|\int_{0}^{1}(1-t)^{r}\left(\nabla^{(\lfloor\beta\rfloor)}u\left(x_{0}+\alpha
tz\right)-\nabla^{(\lfloor\beta\rfloor)}u(x_{0})\right)\cdot
z^{\lfloor\beta\rfloor}\mathrm{d}t\right|\\\
\leq\frac{\alpha^{\lfloor\beta\rfloor}}{r!}\int_{0}^{1}(1-t)^{r}M(\alpha
t)^{\beta-\lfloor\beta\rfloor}\|z\|_{2}^{\beta-\lfloor\beta\rfloor}\|z\|_{\infty}^{\lfloor\beta\rfloor}\mathrm{d}t\leq\frac{M\alpha^{\beta}}{\lfloor\beta\rfloor!}.$

A similar expansion holds for $\tilde{u}$. And since $x_{0}+\alpha
z\in\Omega$, we have $|u(x_{0}+\alpha z)-\tilde{u}(x_{0}+\alpha
z)|\leq\varepsilon$. Therefore, setting

$P(z)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}u(x_{0})-\tilde{u}(x_{0})+\alpha(\nabla
u(x_{0})-\nabla\tilde{u}(x_{0}))\cdot z\\\
+\cdots+\frac{\alpha^{\lfloor\beta\rfloor}}{\lfloor\beta\rfloor!}\left(\nabla^{(\lfloor\beta\rfloor)}u(x_{0})-\nabla^{(\lfloor\beta\rfloor)}\tilde{u}(x_{0})\right)\cdot
z^{\lfloor\beta\rfloor},$

we conclude that

(A.24) $\sup_{z\in\mathbb{R}^{m}:\;\|z\|_{2}\leq
1}|P(z)|\leq\varepsilon+\frac{2M\alpha^{\beta}}{\lfloor\beta\rfloor!}$

We note that the function $z\mapsto P(z)$ is infinitely differentiable
polynomial of degree $\lfloor\beta\rfloor$ on $\mathbb{R}^{m}$, and for
$0\leq\tau\leq\lfloor\beta\rfloor$, and $y\in\mathbb{R}^{m}$,

(A.25) $\nabla^{(\tau)}P(z)\cdot
y^{\tau}=\sum_{j=\tau}^{\lfloor\beta\rfloor}\frac{\alpha^{j}}{(j-\tau)!}\left(\nabla
u^{(j)}(x_{0})-\nabla\tilde{u}^{(j)}(x_{0})\right)\cdot\left(z^{j-\tau},y^{\tau}\right).$

On the other hand, the function
$x\mapsto\nabla^{(\tau)}u(x):\;\Omega\mapsto\mathcal{L}(\mathbb{R}^{m}\times\cdots\times\mathbb{R}^{m};\mathbb{R})$
has derivatives to order $\lfloor\beta\rfloor-\tau$ and
$(\nabla^{(j)}\nabla^{(\tau)}u(x_{0})\cdot y^{j})\cdot
x^{\tau}=\nabla^{(\tau+j)}u(x_{0})\cdot(y^{j},x^{\tau})$. Setting
$k=\lfloor\beta\rfloor-\tau$, its Taylor expansion at $x_{0}$ yields

$\nabla^{(\tau)}u(x_{0}+\alpha
z)=\underbrace{\nabla^{(\tau)}u(x_{0})+\alpha\nabla^{(\tau+1)}u(x_{0})\cdot
z+\cdots+\frac{\alpha^{k}}{k!}\nabla^{(\lfloor\beta\rfloor)}u(x_{0})\cdot
z^{k}}_{Q(z)}+R^{\prime}(z),$

and

$\nabla^{(\tau)}\tilde{u}(x_{0}+\alpha
z)=\underbrace{\nabla^{(\tau)}\tilde{u}(x_{0})+\alpha\nabla^{(s+1)}\tilde{u}(x_{0})\cdot
z+\cdots+\frac{\alpha^{k}}{k!}\nabla^{(\lfloor\beta\rfloor)}\tilde{u}(x_{0})\cdot
z^{k}}_{\tilde{Q}(z)}+\tilde{R}^{\prime}(z).$

Using the same calculations as above, we check that the remainders
$R^{\prime}$ and $\tilde{R}^{\prime}$ satisfy

$\sup_{z\in\mathbb{R}^{d}:\;\|z\|_{2}\leq
1}\;\|R^{\prime}(z)\|+\|\tilde{R}^{\prime}(z)\|\leq\frac{2M\alpha^{\beta-\tau}}{(\lfloor\beta\rfloor-\tau)!}.$

Using the last display, and noting that
$\alpha^{\tau}(Q(z)-\tilde{Q}(z))=\nabla^{(\tau)}P(z)$ given in (A.25), it
follows that for all $z$ in the unit ball,

$\|\nabla^{(\tau)}u(x_{0}+\alpha z)-\nabla^{(\tau)}\tilde{u}(x_{0}+\alpha
z)\|\leq\frac{1}{\alpha^{\tau}}\|\nabla^{(\tau)}P(z)\|+\frac{2M\alpha^{\beta-\tau}}{(\lfloor\beta\rfloor-\tau)!}.$

Therefore, by Markov’s polynomial inequality (Lemma A.13), for all $z$ in the
unit ball

$\|\nabla^{(\tau)}u(x_{0}+\alpha z)-\nabla^{(\tau)}\tilde{u}(x_{0}+\alpha
z)\|\leq\frac{T^{(\tau)}_{\lfloor\beta\rfloor}(1)}{\alpha^{\tau}}\left(\varepsilon+\frac{2M\alpha^{\beta}}{\lfloor\beta\rfloor!}\right)+\frac{2M\alpha^{\beta-\tau}}{(\lfloor\beta\rfloor-\tau)!}\\\
=\frac{A}{\alpha^{\tau}}+B\alpha^{\beta-\tau},$

where

$A=\varepsilon T^{(\tau)}_{\lfloor\beta\rfloor}(1),\;\;B=\frac{2M\times
T^{(\tau)}_{\lfloor\beta\rfloor}(1)}{\lfloor\beta\rfloor!}+\frac{2M}{(\lfloor\beta\rfloor-\tau)!}.$

With the choice

$\alpha=\left(\frac{\tau
A}{(\beta-\tau)B}\right)^{1/\beta}\leq\left(\frac{\tau\varepsilon\lfloor\beta\rfloor!}{2M(\beta-\tau)}\right)^{1/\beta},$

we get

$\|\nabla^{(\tau)}u(x_{0})-\nabla^{(\tau)}\tilde{u}(x_{0})\|\leq
CM^{\frac{\tau}{\beta}}\varepsilon^{\frac{\beta-\tau}{\beta}},$

for some constant $C$ that depends only on $\tau$ and $\beta$. In fact, $C$
can be taken as

$C=\frac{\beta
T^{(\tau)}_{\lfloor\beta\rfloor}(1)}{\beta-\tau}\left(\frac{\beta-\tau}{\tau}\right)^{\tau/\beta}\leq\frac{\beta}{\tau}T^{(\tau)}_{\lfloor\beta\rfloor}(1).$

###  A.7 Proof of Theorem 2.3

We start with a lemma that establishes that the covering number of
$\Gamma(\Theta)$ grows polynomially near $0$, which consequently implies that
$\Gamma(\Theta)$ has a finite VC dimension. This effectively proves that
$\Gamma(\Theta)$ is a finite-dimensional subset of
$\mathcal{C}^{\beta}(\Omega)$, serving as a cornerstone of our main theorem.

######  Lemma A.14.

Given any ${\epsilon}>0$, we have:

$\mathcal{N}\left({\epsilon},\Gamma(\Theta),L_{\infty}\right)\leq
C\left(\frac{c}{{\epsilon}}\right)^{d}\,,$

for some constants $C,c>0$.

######  Proof A.15.

Recall that $\Theta$ is a compact subset of ${\mathbb{R}}^{d}$, which implies
${\sf vol}(\Theta)<\infty$. Therefore, from the standard covering number
calculation, we have:

$\mathcal{N}({\epsilon},\Theta,L_{2})\leq\frac{{\sf vol}(\Theta)}{{\sf
vol}(B)}\left(\frac{3}{{\epsilon}}\right)^{d}\triangleq
C\left(\frac{3}{{\epsilon}}\right)^{d}\,,$

where $B$ is the unit ball (centered at the origin) in ${\mathbb{R}}^{d}$
(e.g., see Theorem 14.2 of Yihong Wu’s lecture notes, ECE598, Spring 2016). We
next claim that:

$\mathcal{N}\left({\epsilon},\Gamma(\Theta),L_{\infty}\right)\leq\mathcal{N}({\epsilon}/L,\Theta,L_{2})\,.$

which will complete the proof. To show this, suppose
$\\{\theta_{1},\dots,\theta_{n}\\}$ are the center of $({\epsilon}/L)$-covers
of $\Theta$, i.e.

$\sup_{\theta\in\Theta}\min_{1\leq j\leq
n}\|\theta_{j}-\theta\|_{2}\leq{\epsilon}/L\,.$

where $n=\mathcal{N}({\epsilon}/L,\Theta,L_{2})$. Now consider the images
corresponding to these $\theta$’s, meaning,
$\\{\gamma_{1},\dots,\gamma_{n}\\}$ where $\gamma_{i}=\Gamma(\theta_{i})$.
Then we have:

$\sup_{\theta\in\Theta}\min_{1\leq j\leq
n}\|\Gamma(\theta_{j})-\Gamma(\theta)\|_{\infty}\leq
L\sup_{\theta\in\Theta}\min_{1\leq j\leq n}\|\theta_{j}-\theta\|_{2}\leq
L({\epsilon}/L)={\epsilon}\,.$

Hence, the proof of the lemma is complete with $c=3L$.

######  Proof A.16 (Proof of Theorem 2.3).

We use the standard rate theorem (e.g., see 3.2.5 of Vaart and Wellner (1997))
to establish the convergence rate of $|\hat{u}-u_{\star}|_{2}^{2}$. Define,
for simplicity, $M(u)=\mathbb{E}[(Y-u({\bf s}))^{2}]$ and
$M_{n}(u)=(1/n)\sum_{i}(Y_{i}-u({\bf s}_{i}))^{2}$. Therefore, from (2.14) we
have $\hat{u}$ is the minimizer of $M_{n}(u)$ over $\Gamma(\Theta)$.
Furthermore, we have:

$M(u)-M(u_{\star})=|u-u_{\star}|_{2}^{2}\,.$

To establish the rate of convergence, we need to find the modulus of
continuity function $\phi_{n}(\delta)$ that satisfies:

$\mathbb{E}\left[\sup_{|u-u_{\star}|_{2}\leq\delta}\left|\left(M_{n}(u)-M(u)\right)-\left(M_{n}(u_{\star})-M(u_{\star})\right)\right|\right]\leq\frac{\phi_{n}(\delta)}{\sqrt{n}}\,.$

Let $\xi_{i}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}(Y_{i}-u_{\star}({\bf
s}_{i}))/\sigma$. Some simple algebra yields:

$\displaystyle\left|\left(M_{n}(u)-M(u)\right)-\left(M_{n}(u_{\star})-M(u_{\star})\right)\right|$
$\displaystyle=\left|(\mathbb{P}_{n}-\mathbb{P})(u({\bf s})-u_{\star}({\bf
s}))^{2}+\sigma(\mathbb{P}_{n}-\mathbb{P})\xi(u({\bf s})-u_{\star}({\bf
s}))\right|$ $\displaystyle\leq\left|(\mathbb{P}_{n}-\mathbb{P})(u({\bf
s})-u_{\star}({\bf
s}))^{2}\right|+\sigma\left|(\mathbb{P}_{n}-\mathbb{P})\xi(u({\bf
s})-u_{\star}({\bf s}))\right|$

Therefore, by symmetrization we have (using $(\eta_{1},\dots,\eta_{n})$ i.i.d
Rademacher random variables):

$\displaystyle\mathbb{E}\left[\sup_{|u-u_{\star}|_{2}\leq\delta}\left|\left(M_{n}(u)-M(u)\right)-\left(M_{n}(u_{\star})-M(u_{\star})\right)\right|\right]$
$\displaystyle\leq\mathbb{E}\left[\sup_{|u-u_{\star}|_{2}\leq\delta}\left|(\mathbb{P}_{n}-\mathbb{P})(u({\bf
s})-u_{\star}({\bf
s}))^{2}\right|\right]+\sigma\mathbb{E}\left[\sup_{|u-u_{\star}|_{2}\leq\delta}\left|(\mathbb{P}_{n}-\mathbb{P})\xi(u({\bf
s})-u_{\star}({\bf s}))\right|\right]$ $\displaystyle\leq
2\mathbb{E}\left[\sup_{|u-u_{\star}|_{2}\leq\delta}\left|\frac{1}{n}\sum_{i}\eta_{i}(u({\bf
s}_{i})-u_{\star}({\bf
s}_{i}))^{2}\right|\right]+\sigma\mathbb{E}\left[\sup_{|u-u_{\star}|_{2}\leq\delta}\left|\frac{1}{n}\sum_{i}\xi_{i}(u({\bf
s}_{i})-u_{\star}({\bf s}_{i}))\right|\right]$ $\displaystyle\leq
2B\mathbb{E}\left[\sup_{|u-u_{\star}|_{2}\leq\delta}\left|\frac{1}{n}\sum_{i}\eta_{i}(u({\bf
s}_{i})-u_{\star}({\bf
s}_{i}))\right|\right]+\sigma\mathbb{E}\left[\sup_{|u-u_{\star}|_{2}\leq\delta}\left|\frac{1}{n}\sum_{i}\xi_{i}(u({\bf
s}_{i})-u_{\star}({\bf s}_{i}))\right|\right]$ $\displaystyle\leq
C_{1}\left(\delta\sqrt{\frac{d}{n}\log{\left(\frac{C_{2}}{\delta}\right)}}+\frac{d}{n}\log{\left(\frac{C_{2}}{\delta}\right)}\right)$

Here, in the second last inequality, we use Leduox-Talagrand contraction
inequality. The last line follows Dudley’s chaining bound and some standard
calculations (e.g., proof of Theorem 8.7 of Sen (2018)). Therefore, a valid
choice of $\phi_{n}(\delta)$ is:

$\phi_{n}(\delta)=C_{1}\left(\delta\sqrt{d\log{\left(\frac{C_{2}}{\delta}\right)}}+\frac{d}{\sqrt{n}}\log{\left(\frac{C_{2}}{\delta}\right)}\right)\,.$

Therefore, we conclude that $\|\hat{u}-u_{\star}\|_{2}=O_{p}(\delta_{n})$
where $\delta_{n}$ satisfies:

$\sqrt{n}\delta_{n}^{2}\geq\phi_{n}(\delta_{n})\,.$

Some simple algebra yields that a valid choice for $\delta_{n}$ is
$\sqrt{(d/n)\log{(n/d)}}$, which completes the proof.

## References

  * Ando and van Hemmen (1980) Ando, T. and van Hemmen, J. L. (1980).  An inequality for trace ideals.  Communications in Mathematical Physics 76 143 – 148\.
  * Atchade and Bhattacharyya (2019) Atchade, Y. and Bhattacharyya, A. (2019).  An approach to large-scale quasi-bayesian inference with spike-and-slab priors.
  * Atchade and Wang (2023) Atchade, Y. and Wang, L. (2023).  A fast asynchronous Markov chain Monte Carlo sampler for sparse Bayesian inference.  Journal of the Royal Statistical Society Series B: Statistical Methodology 85 1492–1516.
  * Belomestny et al. (2022) Belomestny, D., Naumov, A., Puchkin, N. and Samsonov, S. (2022).  Simultaneous approximation of a smooth function and its derivatives by deep neural networks with piecewise-polynomial activations.
  * Bhatia et al. (2019) Bhatia, R., Jain, T. and Lim, Y. (2019).  On the bures-wasserstein distance between positive definite matrices.  Expositiones Mathematicae 37 165–191.
  * Both et al. (2021) Both, G.-J., Choudhury, S., Sens, P. and Kusters, R. (2021).  Deepmod: Deep learning for model discovery in noisy data.  Journal of Computational Physics 428 109985\.
  * Chen et al. (2020) Chen, Z., Liu, Y. and Sun, H. (2020).  Physics-informed learning of governing equations from scarce data.  Nature Communications 12.
  * Cuomo et al. (2022) Cuomo, S., Cola, V. S. D., Giampaolo, F., Rozza, G., Raissi, M. and Piccialli, F. (2022).  Scientific machine learning through physics-informed neural networks: Where we are and what’s next.  Journal of Scientific Computing 92.
  * de Hoop et al. (2023) de Hoop, M. V., Kovachki, N. B., Nelsen, N. H. and Stuart, A. M. (2023).  Convergence rates for learning linear operators from noisy data.  SIAM/ASA Journal on Uncertainty Quantification 11 480–513.
  * Devroye et al. (2023) Devroye, L., Mehrabian, A. and Reddad, T. (2023).  The total variation distance between high-dimensional gaussians with the same mean.
  * Flamary et al. (2021) Flamary, R., Courty, N., Gramfort, A., Alaya, M. Z., Boisbunon, A., Chambon, S., Chapel, L., Corenflos, A., Fatras, K., Fournier, N., Gautheron, L., Gayraud, N. T., Janati, H., Rakotomamonjy, A., Redko, I., Rolet, A., Schutz, A., Seguy, V., Sutherland, D. J., Tavenard, R., Tong, A. and Vayer, T. (2021).  Pot: Python optimal transport.  Journal of Machine Learning Research 22 1–8.
  * Geer (2000) Geer, S. A. (2000).  Empirical Processes in M-estimation, vol. 6.  Cambridge university press.
  * Giordano and Kekkonen (2020) Giordano, M. and Kekkonen, H. (2020).  Bernstein–von mises theorems and uncertainty quantification for linear inverse problems.  SIAM/ASA Journal on Uncertainty Quantification 8 342–373.
  * Harris (2008) Harris, L. (2008).  Multivariate markov polynomial inequalities and chebyshev nodes.  Journal of Mathematical Analysis and Applications 338 350–357.
  * Hwang (1980) Hwang, C.-R. (1980).  Laplace’s Method Revisited: Weak Convergence of Probability Measures.  The Annals of Probability 8 1177 – 1182.
  * Lu et al. (2021) Lu, J., Shen, Z., Yang, H. and Zhang, S. (2021).  Deep network approximation for smooth functions.  SIAM Journal on Mathematical Analysis 53 5465–5506.
  * Lu et al. (2022) Lu, Y., Chen, H., Lu, J., Ying, L. and Blanchet, J. (2022).  Machine learning for elliptic PDEs: Fast rate generalization bound, neural scaling law and minimax optimality.  In International Conference on Learning Representations.
  * McGoff et al. (2015) McGoff, K., Mukherjee, S. and Pillai, N. (2015).  Statistical inference for dynamical systems: A review.  Statistics Surveys 9 209 – 252.
  * Mishra and Molinaro (2021) Mishra, S. and Molinaro, R. (2021).  Estimates on the generalization error of physics-informed neural networks for approximating a class of inverse problems for PDEs.  IMA Journal of Numerical Analysis 42 981–1022.
  * Mishra and Molinaro (2022) Mishra, S. and Molinaro, R. (2022).  Estimates on the generalization error of physics-informed neural networks for approximating PDEs.  IMA Journal of Numerical Analysis 43 1–43.
  * Monard et al. (2021) Monard, F., Nickl, R. and Paternain, G. P. (2021).  Statistical guarantees for Bayesian uncertainty quantification in nonlinear inverse problems with Gaussian process priors.  The Annals of Statistics 49 3255 – 3298.
  * Nickl et al. (2020) Nickl, R., van de Geer, S. and Wang, S. (2020).  Convergence rates for penalized least squares estimators in pde constrained regression problems.  SIAM/ASA Journal on Uncertainty Quantification 8 374–413.
  * Pardo (2006) Pardo, L. (2006).  Statistical Inference Based on Divergence Measures, vol. 185 of Statistics: Textbooks and Monographs.  Chapman & Hall/CRC, Boca Raton, FL.
  * Raissi et al. (2019) Raissi, M., Perdikaris, P. and Karniadakis, G. (2019).  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  Journal of Computational Physics 378 686–707.
  * Ramsay et al. (2007) Ramsay, J. O., Hooker, G., Campbell, D. and Cao, J. (2007).  Parameter estimation for differential equations: A generalized smoothing approach.  Journal of the Royal Statistical Society. Series B (Statistical Methodology) 69 741–796.
  * Rasmussen and Williams (2005) Rasmussen, C. E. and Williams, C. K. I. (2005).  Gaussian Processes for Machine Learning.  Adaptive Computation and Machine Learning, The MIT Press, Cambridge, MA.
  * Schmidt-Hieber (2020) Schmidt-Hieber, J. (2020).  Nonparametric regression using deep neural networks with relu activation function.  Annals of Statistics 48.
  * Sen (2018) Sen, B. (2018).  A gentle introduction to empirical process theory and applications.  Lecture Notes, Columbia University 11 28–29.
  * Shin et al. (2023) Shin, Y., Zhang, Z. and Karniadakis, G. E. (2023).  Error estimates of residual minimization using neural networks for linear pdes.  Journal of Machine Learning for Modeling and Computing 4 73–101.
  * Sirignano and Spiliopoulos (2018) Sirignano, J. and Spiliopoulos, K. (2018).  Dgm: A deep learning algorithm for solving partial differential equations.  Journal of computational physics 375 1339–1364.
  * Spokoiny (2020) Spokoiny, V. (2020).  Bayesian inference for nonlinear inverse problems.
  * Stone (1982) Stone, C. J. (1982).  Optimal Global Rates of Convergence for Nonparametric Regression.  The Annals of Statistics 10 1040 – 1053.
  * Sun et al. (2020) Sun, L., Gao, H., Pan, S. and Wang, J.-X. (2020).  Surrogate modeling for fluid flows based on physics-constrained deep learning without simulation data.  Computer Methods in Applied Mechanics and Engineering 361 112732\.
  * Taheri et al. (2021) Taheri, M., Xie, F. and Lederer, J. (2021).  Statistical guarantees for regularized neural networks.  Neural Networks 142 148–161.
  * Vaart and Wellner (1997) Vaart, A. v. d. and Wellner, J. A. (1997).  Weak convergence and empirical processes with applications to statistics.  Journal of the Royal Statistical Society-Series A Statistics in Society 160 596–608.
  * Xun et al. (2013) Xun, X., Cao, J., Mallick, B., Maity, A. and Carroll, R. J. (2013).  Parameter estimation of partial differential equation models.  Journal of the American Statistical Association 108 1009–1020.
  * Yang et al. (2021) Yang, L., Meng, X. and Karniadakis, G. E. (2021).  B-pinns: Bayesian physics-informed neural networks for forward and inverse pde problems with noisy data.  Journal of Computational Physics 425 109913\.
  * Yarotsky (2017) Yarotsky, D. (2017).  Error bounds for approximations with deep relu networks.  Neural Networks 94 103–114.
  * Yarotsky and Zhevnerchuk (2020) Yarotsky, D. and Zhevnerchuk, A. (2020).  The phase diagram of approximation rates for deep neural networks.  In Proceedings of the 34th International Conference on Neural Information Processing Systems. NIPS’20, Curran Associates Inc.
  * Yu et al. (2018) Yu, B. et al. (2018).  The deep ritz method: a deep learning-based numerical algorithm for solving variational problems.  Communications in Mathematics and Statistics 6 1–12.