will be enough to exhibit a limit object that explains the patterns visible on
Figure 1. This is done in propositions 10 and 11, which are statements on the
regularized log marginal likelihood
$\mathcal{L}_{\mathrm{reg}}(x_{0},\lambda)$ defined in equation (4.62) rather
than on $\mathcal{L}(x_{0},\lambda)$. Note however that proving a complete
mathematical statement linking the behaviours of $\mathcal{L}(x_{0},\lambda)$
and $\mathcal{L}_{\mathrm{reg}}(x_{0},\lambda)$ is lacking for the moment.

####  4.5.2. Regularizing the point source and corresponding approximations

First, we regularize $F_{t}$ thanks to a mollifier $\varphi(x)$ on
$\mathbb{R}^{3}$ which we choose to be radially symmetric as in [26], section
4.2.1. Define $\varphi_{R}(y)=\varphi(y/R)/R^{3}$, then a
$\mathcal{C}^{\infty}_{c}$ regularization of $F_{t}$ is obtained by

$\displaystyle f_{t}^{\mathrm{R}}(x):=(F_{t}*\varphi_{R})(x)\ \ \ \forall
x\in\mathbb{R}^{3}$

As $F_{t}$, $f_{t}^{\mathrm{R}}$ exhibits radial symmetry. We will next use
the following approximated objects :

  * •

Note
$k_{x_{0}}^{\mathrm{reg}}((x,t),(x^{\prime},t^{\prime})):=f_{t}^{\mathrm{R}}(x-x_{0})f_{t^{\prime}}^{\mathrm{R}}(x^{\prime}-x_{0})$
a regularized version of the limit object in Proposition 9. Proposition 9
states that in some sense, when $R$ approaches $0$,
$k_{\mathrm{nor,x_{0}}}^{\mathrm{R}}$ is close to $k_{x_{0}}^{\mathrm{reg}}$.

  * •

We also assume that $w(x_{i},t_{j})$ may be approximated by
$\tilde{w}(x_{i},t_{j})=f_{t_{j}}^{\mathrm{R}}(x_{i}-x_{0}^{*})$ as in the
point source limit, $v_{0}=\delta_{x_{0}^{*}}$ and in that case we would have
$w(x_{i},t_{j})=(F_{t_{j}}*v_{0})(x_{i})=F_{t_{j}}(x_{i}-x_{0}^{*})$
(forgetting for a second that $F_{t}$ is not defined pointwise). We thus
introduce the column vector of “approximated observations”
$W=\big{(}\tilde{w}(x_{i},t_{j})\big{)}_{i,j}$ and we suppose that $W$ is
ordered in the following way :

$\displaystyle W=\begin{pmatrix}W_{1}\\\ \cline{1-1}\cr\vdots\\\
\cline{1-1}\cr W_{q}\end{pmatrix}$

where $W_{i}$ corresponds to the sensor $i$ :
$W_{i}=(\tilde{w}(x_{i},t_{1}),...,\tilde{w}(x_{i},t_{N}))^{T}$.

These approximated objects enable us to introduce the “regularized” log
marginal likelihood built by replacing $k$ with $k_{x_{0}}^{\mathrm{reg}}$ and
$w_{\mathrm{obs}}$ by $W$ :

(4.62)
$\displaystyle\mathcal{L}_{\mathrm{reg}}(x_{0},\lambda):=W^{T}(K_{x_{0}}^{\mathrm{reg}}+\lambda
I_{n})^{-1}W+\log\det(K_{x_{0}}^{\mathrm{reg}}+\lambda I_{n})$

with $K_{x_{0}}^{\mathrm{reg}}=k_{x_{0}}^{\mathrm{reg}}(X,X)$. Instead of
studying $\mathcal{L}(x_{0},\lambda)$, we now rather analyse
$\mathcal{L}_{\mathrm{reg}}(x_{0},\lambda)$; we expect
$\mathcal{L}_{\mathrm{reg}}(x_{0},\lambda)$ to be close to
$\mathcal{L}(x_{0},\lambda)$ in some sense, although here we do not state any
approximation theorem w.r.t. the approximation of $\mathcal{L}(x_{0},\lambda)$
with $\mathcal{L}_{\mathrm{reg}}(x_{0},\lambda)$.

Because of the tensor product shape of $k_{x_{0}}^{\mathrm{reg}}$,
$K_{x_{0}}^{\mathrm{reg}}$ has rank one :
$K_{x_{0}}^{\mathrm{reg}}=F_{x_{0}}F_{x_{0}}^{T}$ with

$\displaystyle F_{x_{0}}=\begin{pmatrix}F_{x_{0}}^{1}\\\
\cline{1-1}\cr\vdots\\\ \cline{1-1}\cr F_{x_{0}}^{q}\end{pmatrix}$

and
$F_{x_{0}}^{i}=(f_{t_{1}}^{\mathrm{R}}(x_{i}-x_{0}),...,f_{t_{N}}^{\mathrm{R}}(x_{i}-x_{0}))^{T}$.
Below, we write a proposition which describes the asymptotic behaviour of
$\mathcal{L}_{\mathrm{reg}}(x_{0},\lambda)$ in the limit of
$\lambda\rightarrow 0$. This limit corresponds to noiseless observations; we
will see that the limit object in Proposition 10 provides an explanation of
the patterns of Figure 1. We note $||\cdot||$ and $\langle\cdot,\cdot\rangle$
the usual Euclidean norm and dot product in $\mathbb{R}^{n}$.

######  Proposition 10 (Asymptotic behaviour of
$\mathcal{L}_{\mathrm{reg}}(x_{0},\lambda)$ when $\lambda\rightarrow 0$).

Let $\varepsilon>0$ and
$E_{\varepsilon}:=\\{x_{0}\in\mathbb{R}^{3}:\rho(K_{x_{0}}^{\mathrm{reg}})>\varepsilon\\}$
where $\rho(A)$ denotes the spectral radius of $A$. Define the correlation
coefficient between $F_{x_{0}}$ and $W$ by
$r(x_{0})=\text{Corr}(F_{x_{0}},W)=\langle
F_{x_{0}},W\rangle/(||W||\times||F_{x_{0}}||)\in[-1,1]$. We set $r(x_{0})=0$
if $F_{x_{0}}=0$.

Then we have the pointwise convergence

$\displaystyle\forall x_{0}\in\mathbb{R}^{3},\ \ \
\big{|}\lambda\mathcal{L}_{\mathrm{reg}}(x_{0},\lambda)-||W||^{2}\big{(}1-r(x_{0})^{2}\big{)}\big{|}=O_{\lambda\rightarrow
0}(\lambda\log\lambda)$

and the uniform convergence on $E_{\varepsilon}$

$\displaystyle\sup_{x_{0}\in
E_{\varepsilon}}\big{|}\lambda\mathcal{L}_{\mathrm{reg}}(x_{0},\lambda)-||W||^{2}\big{(}1-r(x_{0})^{2}\big{)}\big{|}=O_{\lambda\rightarrow
0}(\lambda\log\lambda)$

###### Proof.

Suppose first that $||F_{x_{0}}||^{2}=0$. Then obviously, $r(x_{0})=0$ and
$\mathcal{L}_{\mathrm{reg}}(x_{0},\lambda)=||W||^{2}/\lambda+n\log\lambda$
which shows that

(4.63)
$\displaystyle\big{|}\lambda\mathcal{L}_{\mathrm{reg}}(x_{0},\lambda)-||W||^{2}\big{|}=O_{\lambda\rightarrow
0}(\lambda\log\lambda)$

Now, suppose that $||F_{x_{0}}||^{2}>0$. Let $\varepsilon>0$ such that
$||F_{x_{0}}||^{2}\geq\varepsilon$.

We first deal with the first term in (4.62). Using the Sherman–Morrison
formula ([55], section 2.7.1), we may invert
$(K_{x_{0}}^{\mathrm{reg}}+\lambda I_{n})$ explicitly :

$\displaystyle(K_{x_{0}}^{\mathrm{reg}}+\lambda
I_{n})^{-1}=\frac{1}{\lambda}I_{n}-\frac{1}{\lambda^{2}}\frac{F_{x_{0}}F_{x_{0}}^{T}}{1+\frac{1}{\lambda}F_{x_{0}}^{T}F_{x_{0}}}=\frac{1}{\lambda}\Big{(}I_{n}-\frac{F_{x_{0}}F_{x_{0}}^{T}}{\lambda+||F_{x_{0}}||^{2}}\Big{)}$

One also easily derives the determinant term as $F_{x_{0}}F_{x_{0}}^{T}$ only
has one non zero eigenvalue equal to $||F_{x_{0}}||^{2}$, since
$(F_{x_{0}}F_{x_{0}}^{T})F_{x_{0}}=F_{x_{0}}(F_{x_{0}}^{T}F_{x_{0}})=||F_{x_{0}}||^{2}F_{x_{0}}$
:

(4.64) $\displaystyle\log\det(K_{x_{0}}^{\mathrm{reg}}+\lambda
I_{n})=(n-1)\log\lambda+\log\lambda+||F_{x_{0}}||^{2})$

(The same argument shows that
$\rho(K_{x_{0}}^{\mathrm{reg}})=||F_{x_{0}}||^{2}$.) Thus,

$\displaystyle\mathcal{L}_{\mathrm{reg}}(x_{0},\lambda)$
$\displaystyle=W^{T}(K_{x_{0}}^{\mathrm{reg}}+\lambda
I_{n})^{-1}W+\log\det(K_{x_{0}}^{\mathrm{reg}}+\lambda I_{n})$
$\displaystyle=\frac{1}{\lambda}\bigg{(}||W||^{2}-\frac{\langle
F_{x_{0}},W\rangle^{2}}{\lambda+||F_{x_{0}}||^{2}}\bigg{)}+(n-1)\log\lambda+\log(\lambda+||F_{x_{0}}||^{2})$
(4.65) $\displaystyle=\frac{||W||^{2}}{\lambda}\bigg{(}1-\frac{\langle
F_{x_{0}},W\rangle^{2}}{||W||^{2}(\lambda+||F_{x_{0}}||^{2})}\bigg{)}+(n-1)\log\lambda+\log(\lambda+||F_{x_{0}}||^{2})$

and

$\displaystyle\lambda\mathcal{L}_{\mathrm{reg}}(x_{0},\lambda)$
$\displaystyle-||W||^{2}(1-r(x_{0})^{2})$ (4.66)
$\displaystyle=||W||^{2}\bigg{(}\frac{\langle
F_{x_{0}},W\rangle^{2}}{||W||^{2}||F_{x_{0}}||^{2}}-\frac{\langle
F_{x_{0}},W\rangle^{2}}{||W||^{2}(\lambda+||F_{x_{0}}||^{2})}\bigg{)}$
$\displaystyle\hskip
15.0pt+(n-1)\lambda\log\lambda+\lambda\log(\lambda+||F_{x_{0}}||^{2})$

But for the term in (4.5.2) that is multiplied by $||W||^{2}$,

$\displaystyle\frac{\langle
F_{x_{0}},W\rangle^{2}}{||W||^{2}||F_{x_{0}}||^{2}}-\frac{\langle
F_{x_{0}},W\rangle^{2}}{||W||^{2}(\lambda+||F_{x_{0}}||^{2})}$
$\displaystyle=\frac{\langle
F_{x_{0}},W\rangle^{2}}{||W||^{2}}\bigg{(}\frac{1}{||F_{x_{0}}||^{2}}-\frac{1}{\lambda+||F_{x_{0}}||^{2}}\bigg{)}$
$\displaystyle=\frac{\langle
F_{x_{0}},W\rangle^{2}}{||W||^{2}}\frac{\lambda}{||F_{x_{0}}||^{2}(\lambda+||F_{x_{0}}||^{2})}$
(4.67) $\displaystyle\leq
r(x_{0})^{2}\frac{\lambda}{\lambda+||F_{x_{0}}||^{2}}\leq\frac{\lambda}{||F_{x_{0}}||^{2}}\leq\frac{\lambda}{\varepsilon}$

and obviously, since $\lambda>0$,

(4.68) $\displaystyle\frac{\langle
F_{x_{0}},W\rangle^{2}}{||W||^{2}||F_{x_{0}}||^{2}}-\frac{\langle
F_{x_{0}},W\rangle^{2}}{||W||^{2}(\lambda+||F_{x_{0}}||^{2})}\geq 0$

Also, one sees that $F_{x_{0}}=0$ as soon as $\sup_{i}|x_{0}-x_{i}|>cT+R$, ie
$x_{0}$ is too far from the receivers for them to capture non zero signal
during the time interval $[0,T]$. Thus the function
$x_{0}\longmapsto||F_{x_{0}}||^{2}$ is zero outside of a compact set. It is
obviously continuous on $\mathbb{R}^{3}$ and is thus bounded on
$\mathbb{R}^{3}$ by some constant $M>0$. Using this together with (4.5.2) and
(4.68) inside (4.5.2) and assuming that $\lambda\leq 1$ yields

(4.69)
$\displaystyle\big{|}\lambda\mathcal{L}_{\mathrm{reg}}(x_{0},\lambda)-||W||^{2}(1-r(x_{0})^{2})\big{|}\leq\frac{\lambda}{\varepsilon}||W||^{2}+(n-1)|\lambda\log\lambda|+\lambda\log(M+1)$

which shows the uniform convergence statement as well as the pointwise one
(together with (4.63)). ∎

The set $E_{\varepsilon}$ is a set on which all the matrices
$K_{x_{0}}^{\mathrm{reg}}$ are all uniformly large enough in the spectral
radius sense.

######  Remark 4.7.

In the proof of Proposition 10, the second term in (4.62) has no influence in
the limit object and only pollutes the rate of convergence. Discarding it
leads to a $O_{\lambda\rightarrow 0}(\lambda)$ rate of convergence.

It also makes sense to inspect the limiting case $N\rightarrow\infty$, which
corresponds to having the sampling frequency of the captors go to infinity. In
this case, the discrete objects in Proposition 10 behave as Riemann sums if
the time steps $t_{k}$ are equally spaced; as a result, we obtain integrals in
the limit $N\rightarrow\infty$. We state a weak asymptotic behaviour result of
$\mathcal{L}_{\mathrm{reg}}(x_{0},\lambda)$ in that limit, with a very similar
limit object. Highlight the dependence in $N$ in
$\mathcal{L}_{\mathrm{reg}}(x_{0},\lambda)$ by noting it instead
$\mathcal{L}_{\mathrm{reg}}^{N}(x_{0},\lambda)$.

######  Proposition 11 (Asymptotic behaviour of
$\mathcal{L}_{\mathrm{reg}}^{N}(x_{0},\lambda)$ when $N\rightarrow\infty$).

Define the following vector valued functions $L^{2}([0,T],\mathbb{R}^{q})$:

$\displaystyle\forall t\in[0,T],\hskip 11.0ptI_{\mathrm{u}}(t)$
$\displaystyle:=\big{(}\tilde{u}(x_{1},t),...,\tilde{u}(x_{q},t)\big{)}^{T}$
$\displaystyle\forall t\in[0,T],\hskip 7.0ptI_{x_{0}}(t)$
$\displaystyle:=\big{(}f_{t}^{\mathrm{R}}(x_{1}-x_{0}),...,f_{t}^{\mathrm{R}}(x_{q}-x_{0})\big{)}^{T}$

and the scalar valued function

(4.70) $\displaystyle r_{\infty}(x_{0})$ $\displaystyle:=\frac{\langle
I_{\mathrm{u}},I_{x_{0}}\rangle}{||I_{\mathrm{u}}||\ ||I_{x_{0}}||}$

whenever $||I_{x_{0}}||\neq 0$; in (4.70), the norms and the dot product are
those of the usual euclidean structure of $L^{2}([0,T],\mathbb{R}^{q})$.
Suppose that for all $k\in\\{1,...,N\\}$, $t_{k}=T(k-1)/(N-1)$, i.e. the
$t_{k}$ are equally spaced in $[0,T]$.

Then whenever $x_{0}$ is such that $||I_{x_{0}}||\neq 0$, we have the
following pointwise convergence at $x_{0}$

(4.71)
$\displaystyle\frac{\lambda}{N}\mathcal{L}_{\mathrm{reg}}^{N}(x_{0},\lambda)\xrightarrow[N\rightarrow\infty]{}||I_{\mathrm{u}}||_{2}^{2}\big{(}1-r_{\infty}(x_{0})^{2}\big{)}+q\lambda\log\lambda$

###### Proof.

In all concerned mathematical objects, we highlight the $N$ dependency with an
exponent, i.e. $W^{N}$, $F_{x_{0}}^{N}$, etc. We use the exact same tool as in
the previous proof, namely that we have

$\displaystyle\mathcal{L}_{\mathrm{reg}}^{N}(x_{0},\lambda)$
$\displaystyle=\frac{||W^{N}||^{2}}{\lambda}\bigg{(}1-\frac{\langle
F_{x_{0}}^{N},W^{N}\rangle^{2}}{||W^{N}||^{2}\big{(}\lambda+||F_{x_{0}}^{N}||^{2}\big{)}}\bigg{)}+(n-1)\log\lambda+\log(\lambda+||F_{x_{0}}^{N}||^{2})$

But

  * •

$||W^{N}||^{2}=\sum_{i=1}^{q}\sum_{k=1}^{N}\tilde{u}(x_{i},t_{k})^{2}$

  * •

$||F_{x_{0}}^{N}||^{2}=\sum_{i=1}^{q}\sum_{k=1}^{N}f_{t_{k}}^{\mathrm{R}}(x_{i}-x_{0})^{2}$

  * •

$\langle
F_{x_{0}}^{N},W^{N}\rangle=\sum_{i=1}^{q}\sum_{k=1}^{N}f_{t_{k}}^{\mathrm{R}}(x_{i}-x_{0})\tilde{u}(x_{i},t_{k})$

Since the time steps are equally spaced, we can study the limit
$N\rightarrow\infty$ of the above objects thanks to Riemann sums. When
$N\rightarrow\infty$,

(4.72) $\displaystyle\frac{1}{N}||W^{N}||^{2}$
$\displaystyle\longrightarrow\sum_{i=1}^{q}\int_{0}^{T}\tilde{u}(x_{i},t)^{2}dt=||I_{\mathrm{u}}||^{2}$
(4.73) $\displaystyle\frac{1}{N}||F_{x_{0}}^{N}||^{2}$
$\displaystyle\longrightarrow\sum_{i=1}^{q}\int_{0}^{T}f_{t}(x_{i}-x_{0})^{2}dt=||I_{x_{0}}||^{2}$
(4.74) $\displaystyle\frac{1}{N}\langle W^{N},F_{x_{0}}^{N}\rangle$
$\displaystyle\longrightarrow\sum_{i=1}^{q}\int_{0}^{T}\tilde{u}(x_{i},t)f_{t}(x_{i}-x_{0})dt=\langle
I_{\mathrm{u}},I_{x_{0}}\rangle$

Suppose that $x_{0}$ is such that $||I_{x_{0}}||\neq 0$, then because of
(4.73), the quantity $||F_{x_{0}}^{N}||$ is bounded from below by
$||I_{x_{0}}||/2$ for sufficiently large $N$. We then have the following
convergence

(4.75) $\displaystyle\frac{\langle
F_{x_{0}}^{N},W^{N}\rangle^{2}}{||W^{N}||^{2}(\lambda+||F_{x_{0}}^{N}||^{2})}=\frac{(\frac{1}{N}\langle
F_{x_{0}}^{N},W^{N}\rangle)^{2}}{\frac{1}{N}||W^{N}||^{2}(\frac{\lambda}{N}+\frac{1}{N}||F_{x_{0}}^{N}||^{2})}\xrightarrow[N\rightarrow\infty]{}r_{\infty}(x_{0})$

Likewise, since $n=qN$, when $N\rightarrow\infty$ :

$\displaystyle\frac{(n-1)\log\lambda}{N}$
$\displaystyle+\frac{1}{N}\log(\lambda+||F_{x_{0}}||^{2})$ (4.76)
$\displaystyle=\frac{(Nq-1)\log\lambda}{N}+\frac{\log
N}{N}+\frac{1}{N}\log\Big{(}\frac{\lambda}{N}+\frac{1}{N}||F_{x_{0}}||^{2}\Big{)}\xrightarrow[N\rightarrow\infty]{}q\log\lambda$

which shows the announced result. ∎

######  Remark 4.8.

One could surely obtain better rates of convergence in Proposition 11 using
basic convergence results for Riemann sums and regularity assumptions, but we
restrict ourselves to a $o_{N\rightarrow\infty}(1)$ convergence to avoid
further technicalities.

####  4.5.3. Discussion : why do we observe that the actual source position is
located at the intersection of spheres centered on receivers in Figure 1 ?

The answer can be found by analysing the limit term in proposition 10 (or the
one in Proposition 11). Note $L(x_{0})$ the said limit object from Proposition
10 :

(4.77) $\displaystyle
L(x_{0})=||W||^{2}\big{(}1-r(x_{0})^{2}\big{)}=||W||^{2}\bigg{(}1-\frac{{}_{B}ig(\sum_{i=1}^{q}\langle
F_{x_{0}}^{i},W_{i}\rangle\Big{)}^{2}}{||W||^{2}||F_{x_{0}}||^{2}}\bigg{)}$

Note $T_{i}$ the time of arrival of the point source wave at sensor $i$ :
$|x_{i}-x_{0}^{*}|=c^{*}T_{i}$. Define also the sphere $S_{i}$ centered on
$x_{i}$ and the thin spherical shell $A_{i}$ of thickness $2R$ that surrounds
it

(4.78) $\displaystyle S_{i}$ $\displaystyle:=S(x_{i},cT_{i})$ (4.79)
$\displaystyle A_{i}$
$\displaystyle:=A(x_{i},cT_{i}-R,cT_{i}+R)=\overline{B(x_{0},cT_{i}+R)\setminus
B(x_{0},cT_{i}-R)}$

(i) $L(x_{0})$ reaches a local minima over the whole sphere $S_{i}$. When
$x_{0}$ is located inside $A_{i}$, the subvectors $W_{i}$ and $F_{x_{0}}^{i}$
of $W$ and $F_{x_{0}}$ respectively become almost colinear because
$f_{t}^{\mathrm{R}}$ is radially symmetric. They become exactly colinear when
$x_{0}\in S_{i}$. This maximizes the term $\langle F_{i},W_{i}\rangle$ in
virtue of the Cauchy-Schwarz inequality. When $x_{0}$ lies on one and only one
of those spherical shells $A_{i}$, the other terms $\langle
F_{j},W_{j}\rangle$ are all zero.

(ii) The local minima of $L(x_{0})$ located at the intersection of $2$ or more
spheres $S_{i}$ are smaller. More generally, when $I$ is a subset of
$\\{1,...,n\\}$ and when $x_{0}\in\bigcap_{i\in
I}A_{i}\setminus\bigcap_{j\notin I}A_{j}$, the term $\sum_{i\in I}\langle
F_{i},W_{i}\rangle$ is (almost) maximized while $\sum_{j\notin I}\langle
F_{j},W_{j}\rangle=0$, which explains why the intersection of spheres are
darker coloured than the other parts of the spheres in Figure 1.

(iii) The spheres $S_{i}$ intersect at a single point, which is exactly
$x_{0}^{*}$ as well as the global minima of $L(x_{0})$. The quantity
$r(x_{0})$ reaches a global maximum when all subvectors $W_{i}$ and
$F_{x_{0}}^{i}$ are colinear, which is the case only when
$x_{0}\in\bigcap_{i}S_{i}$. When there are at least 4 sensors, the
intersection of all the spheres $\bigcap_{i}S_{i}$ is reduced to at most one
point. Recall that we have supposed that $c=c^{*}$ : this implies that
$x_{0}^{*}\in\bigcap_{i}S_{i}$, and thus the minimum of $L(x_{0})$ is located
at $x_{0}=x_{0}^{*}$.

Note that if the speed $c$ in $k_{x_{0}}^{\mathrm{R}}$ does not correspond to
the real speed $c^{*}$, the intersection $\bigcap_{i}S_{i}$ will be empty.
Additionally, from an optimization point of view, numerically solving
$\inf_{x_{0}}\mathcal{L}(x_{0},\lambda)$ is obviously highly non convex and
none of our numerical experiments lead to the correct solution.

###  4.6. Computational Speedups

In the Kriging formulas (2.3) and (2.4) as well as in the log marginal
likelihood equation (2.9), we are interested in inverting matrices of the form
$K+\lambda I$ where $K$ in positive semidefinite. Using basic algebra, we show
an exact computational shortcut to computing $(K+\lambda I)^{-1}$ when $K$ has
null lines, which is typically the case when using kernels (4.56) because of
finite speed propagation. It may provide huge computational speedups in the
limit of small source size as is the case in the previous subsection.

####  4.6.1. Inversion shortcut for Tichonov-regularized PSD matrices with
null columns

Let $K$ be a symmetric semi-positive definite matrix of size $n$. Suppose it
has $p$ non zero columns $i_{1},...,i_{p}$ and $q=n-p$ zero columns
$i_{p+1},...,i_{n}$. Since $K$ is symmetric, its lines $i_{1},...,i_{p}$ are
also zero. Finally, a simple necessary and sufficient condition for its column
$j$ to be zero is $K_{jj}=0$ because of the Cauchy-Schwarz inequality for PSD
matrices : $|a_{ij}|\leq\sqrt{a_{ii}}\sqrt{a_{jj}}$. Define the matrix
$\tilde{K}$ of the extracted non-zero columns of $K$, of size $p\times p$ by
$\tilde{K}_{k,l}=K_{i_{k},i_{l}}\forall k,l\in\\{1,...,p\\}$. Note
$\mathcal{E}=(e_{1},...,e_{n})$ the canonical basis of $\mathbb{R}^{n}$ and
$\mathcal{F}=(e_{i_{1}},...,e_{i_{n}})$ a permutation according to the indexes
$i_{1},...,i_{p}$ and $i_{p+1},...,i_{n}$ defined above (that permutation is
unique up to additional permutations of the null columns of $K$). Then for the
orthogonal permutation matrix $P$ such that

(4.80) $\displaystyle\forall k\in\\{1,...,n\\},\ \ \ Pe_{k}=e_{i_{k}}$

we have

(4.81) $\displaystyle
K=P^{T}\times\begin{pmatrix}\begin{matrix}\tilde{K}\end{matrix}&\vline&\mbox{\large
0}_{p,q}\\\ \hline\cr\mbox{\large 0}_{q,p}&\vline&\begin{matrix}\mbox{\large
0}_{q,q}\end{matrix}\end{pmatrix}\times P$

Then we have blockwise,

(4.82) $\displaystyle(K+\lambda
I_{n})^{-1}=P^{T}\times\begin{pmatrix}\begin{matrix}(\tilde{K}+\lambda
I_{p})^{-1}\end{matrix}&\vline&\mbox{\large 0}_{p,q}\\\ \hline\cr\mbox{\large
0}_{q,p}&\vline&\begin{matrix}\frac{1}{\lambda}I_{q}\end{matrix}\end{pmatrix}\times
P$

Note that matrix $P$ is easy to find as one only has to compute the diagonal
elements of $K$ in order to find the null lines/columns of $K$, since $K$ is
PSD. That is, we only need to invert a $p\times p$ matrix. For $a$ in
$\mathbb{R}^{n}$, define $a_{in}\in\mathbb{R}^{p}$ the first $p$ components of
$Pa$ and $a_{out}\in\mathbb{R}^{n-p}$ the last $n-p$ components of $Pa$, ie

(4.83) $\displaystyle Pa=\begin{pmatrix}(Pa)_{1}\\\ \vdots\\\ (Pa)_{p}\\\
\cline{1-1}\cr(Pa)_{p+1}\\\ \vdots\\\
(Pa)_{n}\end{pmatrix}=\begin{pmatrix}{a}_{in}\\\ \cline{1-1}\cr
a_{out}\end{pmatrix}$

Using (4.82) and (4.83), we have immediately that

(4.84) $\displaystyle a^{T}(K+\lambda
I_{n})^{-1}a=a_{in}^{T}(\tilde{K}+\lambda
I_{p})^{-1}a_{in}+\frac{1}{\lambda}||a_{out}||^{2}$

In the term $a^{T}(K+\lambda I_{n})^{-1}a$, $a_{in}$ corresponds to the
coordinates of $a$ that are ”kept” by the non zero elements of $K$ and
$a_{out}$ corresponds to those that are left out by the non zero elements of
$K$.

####  4.6.2. Computational shortcuts for Kriging

We may now use the previous section to speed up GPR computations when using a
covariance kernel $k_{\theta}$. This is summarized in Proposition 12 in which
we use the following notations.

  * •

Note $X=(x_{1},...,x_{n})^{T}\in\mathbb{R}^{n}$ the observation locations,
$Y\in\mathbb{R}^{n}$ the vector of observations ($Y_{i}=u(x_{i})$ for
noiseless observations).

  * •

Note $K_{\theta}:=k_{\theta}(X,X)$ the associated covariance matrix,
$P_{\theta}$ and $\tilde{K}_{\theta}$ the matrices $P$ and $\tilde{K}$ defined
by (4.80) for $K=K_{\theta}$, $p_{\theta}$ the number of non zero columns of
$K_{\theta}$ and $q_{\theta}=n-p_{\theta}$.

  * •

As in (4.83), note $X_{in}^{\theta}$ the first $p_{\theta}$ coordinates of
$P_{\theta}X$ and $X_{out}^{\theta}$ the $q_{\theta}$ last, and likewise for
$Y_{in}^{\theta}$ and $Y_{out}^{\theta}$.

Plugging equation (4.82) in the log-marginal likelihood and Kriging formulas,
we obtain Proposition 12.

######  Proposition 12.

If using a covariance kernel $k_{\theta}$ with hyperparameter $\theta$ for
GPR, we have the following exact computational shortcuts.
$(i)$ For the log-marginal likelihood,

$\displaystyle\mathcal{L}(\theta,\lambda)$
$\displaystyle=(Y_{in}^{\theta})^{T}(\tilde{K}_{\theta}+\lambda
I_{p_{\theta}})^{-1}Y_{in}^{\theta}$ (4.85)
$\displaystyle+\frac{1}{\lambda}||Y_{out}^{\theta}||^{2}+\log\det(\tilde{K}_{\theta}+\lambda
I_{p_{\theta}})+q_{\theta}\log\lambda$

$(ii)$ For the Kriging formulas,

(4.86) $\displaystyle\tilde{m}(x)$ =
$m(x)+k_{\theta}(X_{in}^{\theta},x)^{T}(\tilde{K}_{\theta}+\lambda
I_{p_{\theta}})^{-1}(Y_{in}^{\theta}-m(X_{in}^{\theta}))$ (4.87)
$\displaystyle\tilde{k}(x,x^{\prime})$ =
$k_{\theta}(x,x^{\prime})-k_{\theta}(X_{in}^{\theta},x)^{T}(\tilde{K}_{\theta}+\lambda
I_{p_{\theta}})^{-1}k_{\theta}(X_{in}^{\theta},x^{\prime})$

$(iii)$ if $k_{\theta}(x,x)=0$ then $\tilde{m}(x)=m(x)$ and
$\tilde{k}(x,x^{\prime})=0$.

###### Proof.

As in equation (4.84), we have immediately that

(4.88) $\displaystyle Y^{T}(K_{\theta}+\lambda
I_{n})^{-1}Y=(Y_{in}^{\theta})^{T}(\tilde{K}_{\theta}+\lambda
I_{p_{\theta}})^{-1}Y_{in}^{\theta}+\frac{1}{\lambda}||Y_{out}^{\theta}||^{2}$

It is also obvious from (4.81) that $\log\det(K_{\theta}+\lambda
I_{n})=\log\det(\tilde{K}_{\theta}+\lambda
I_{p_{\theta}})+q_{\theta}\log\lambda$, which proves (12).

For the Kriging equations, suppose without loss of generality that $m\equiv
0$. The Kriging mean is, as in (4.84),

$\displaystyle\tilde{m}(x)$ $\displaystyle=k(X,x)^{T}(K+\lambda I_{n})^{-1}Y$
(4.89) $\displaystyle=k(X_{in}^{\theta},x)^{T}(\tilde{K}_{\theta}+\lambda
I_{p_{\theta}})^{-1}Y_{in}^{\theta}+\frac{1}{\lambda}k(X_{out}^{\theta},x)^{T}Y_{out}^{\theta}$

Now, suppose that the observation location $x_{i}$ is such that
$k_{\theta}(x_{i},x_{i})=(K_{\theta})_{ii}=0$, which is exactly the case of
all the coordinates of $X_{out}^{\theta}$ (by definition). Then from the
Cauchy-Schwarz inequality (2.2), the function $x\longmapsto k(x,x_{i})$ is
identically zero. Thus for all $x$, we have
$k(X_{out}^{\theta},x)^{T}Y_{out}^{\theta}=0$ which, when plugged in (4.6.2),
proves (4.86). The exact same reasoning leads to (4.87). Finally, for the same
reason, if $k_{\theta}(x,x)=0$ then the vector $k_{\theta}(X_{in},x)$ is null,
which implies that $\tilde{m}(x)=0$. ∎

We remind that formulas (12), (4.86) and (4.87) are exact. There are actually
three shortcuts for the Kriging equations (4.86) and (4.87). First, the fast
matrix inversion, which is especially useful during the log marginal
likelihood optimization step (hyperparameter calibration). Second, the
selection of the pertinent data points in $k_{\theta}(X_{in},x)$ by using
$X_{in}$ instead of $X$. Third, the selection of the pertinent evaluation
points thanks to the point $(iii)$ of Proposition 12. These last two are
useful when $\tilde{m}$ and/or $\tilde{k}$ have to be evaluated many times.
Note also that all the mathematical objects on the right-hand side of
equations (12), (4.86) and (4.87) depend on $\theta$ and have to be recomputed
for each new value of $\theta$. However this is fast since checking which
columns of $K_{\theta}$ are zero only requires to compute the diagonal of
$K_{\theta}$ : this requires $n$ calls to $k_{\theta}$.

######  Remark 4.9.

In equation (12), the term $\frac{1}{\lambda}||Y_{out}^{\theta}||^{2}$ can be
understood as a very heavy penalization term that sanctions non zero
observations that do not lie in the support of $k_{\theta}$. When applying
this to kernel (4.10), the term $\frac{1}{\lambda}||Y_{out}^{\theta}||^{2}$
penalizes the non zero observations that do not lie in the space time domain
(aka light cone) allowed by Huygens’ principle.

####  4.6.3. Expected numerical gain when using the truncated wave kernel
(4.10)

In this case, to examine if a point $(x,t)$ is such that
$k_{\theta}((x,t),(x,t))=0$, one only needs to check that $(x,t)$ lies in the
light cone that corresponds to $\theta$ (among others, $\theta$ contains
$(x_{0},R,c)$). Explicitly, the condition to check is whether or not $x\in
B(x_{0},c|t|+R)\setminus B(x_{0},(c|t|-R)_{+})$ or equivalently,
$c|t|-R\leq|x-x_{0}|\leq c|t|+R$. This is actually much cheaper than
evaluating $k_{\theta}((x,t),(x,t))$.

For the log-marginal likelihood : using Cholesky decomposition to invert
$\tilde{K}_{\theta}+\lambda I_{p_{\theta}}$, only $p_{\theta}^{3}/3$ flops are
required versus $n^{3}/3$ previously (in both cases, the computation of the
determinant term in (12) is immediate from the Cholesky decomposition).

For the Kriging formulas : here we suppose that we are performing GPR using
the truncated wave kernel (4.10). Suppose that the hyperparameter vector
$\theta$ has been previously estimated : in particular, it provides estimated
wave speed $c$, radius $R$ and source position $x_{0}$. Now, equation (4.86)
states that

(4.90)
$\displaystyle\tilde{m}(x,t)=\sum_{i=1}^{p_{\theta}}a_{i}k_{\theta}\big{(}(x,t),(X_{in}^{\theta})_{i}\big{)}$

for some fixed real numbers $a_{1},...,a_{p_{\theta}}$. Suppose that we want
to compute $\tilde{m}(x,t)$ over a cubic spatial grid $[0,L]^{3}$ of side
length $L$ at a fixed time $t$. Suppose that the space step $\Delta x$ of this
grid verifies $\Delta x\ll L$ : then the number of grid points inside a volume
$V$ is approximately $V/\Delta x^{3}$ (each grid point occupies a cubic volume
of $\Delta x^{3}$). In compliance with the strong Huygens principle, computing
$\tilde{m}(x,t)$ over the whole grid using Proposition 12 results in
$N_{\theta}$ calls to $k$, with

(4.91) $\displaystyle N_{\theta}\simeq
p_{\theta}\times\frac{\text{Vol}\Big{(}[B(x_{0},ct+R)\setminus
B(x_{0},(ct-R)_{+})]\cap[0,L]^{3}\Big{)}}{\Delta x^{3}}$

The set $B(x_{0},ct+R)\setminus B(x_{0},(ct-R)_{+})$ corresponds to the
admissible domain allowed by the strong Huygens principle at time $t$. This is
to be compared to $N_{\mathrm{full}}=n\times L^{3}/\Delta x^{3}$, the number
of calls to $k$ needed without using Proposition 12; obviously,
$N_{\theta}<N_{\mathrm{full}}$ and actually, $N_{\theta}\ll N_{\mathrm{full}}$
when $R$ is small. More precise estimations of the number of grid points
inside the volume in equation (4.91) can be derived, though this is a
difficult research topic in itself : see for instance [7] or the Gauss circle
problem for its 2D counterpart. Also, the volume appearing in (4.91) can be
explicitly derived but we consider the additional insights this provides to be
negligible.

##  5\. Numerical experiments

In this section, we showcase WIGPR on functions $w$ that are solutions of
problem (4.1), using kernels (4.47) and (4.48) separately as well as together
as in (4.10). The goal is twofold : reconstructing the target function $w$,
which amounts to reconstructing its initial conditions, and estimating the
physical parameters attached. Note that with badly estimated physical
parameters, the reconstruction step is more or less bound to fail, especially
with inaccurate wave speed $c$ and/or source centers $x_{0}^{u}$ and
$x_{0}^{v}$. Throughout this section, we only consider compactly supported and
radially symmetric initial conditions, for which we can use the formulas
(4.56) and (4.57).

At the end of this section, we also briefly show that perfect transparent
boundary conditions are automatically enforced on the numerical
reconstructions of $w$ provided by WIGPR. This is true simply because the
Kriging means provided by WIGPR are built to be exact solutions of the free
space problem (4.1). Transparent boundary conditions are a well known research
topic in the numerical analysis community ([11],[24], [35],[39]), which arises
from the need to numerically simulate free space (unbounded) PDE problems on
numerical simulation domains that are inherently finite.

###  5.1. General test protocol for WIGPR

####  5.1.1. Test cases

Three general test cases for WIGPR are considered here. They correspond to
different initial conditions $(u_{0},v_{0})$ :

  1. (1)

A test case for $k_{\mathrm{u}}^{\mathrm{wave}}$ described in Subsection 5.2,
for which $v_{0}=0$.

  2. (2)

A test case for $k_{\mathrm{v}}^{\mathrm{wave}}$ described in Subsection 5.3,
for which $u_{0}=0$.

  3. (3)

A test case for
$k_{\mathrm{u}}^{\mathrm{wave}}+k_{\mathrm{v}}^{\mathrm{wave}}$ described in
Subsection 5.4 using both of the non zero initial conditions of the two
previous test cases.

For each test case, the complete procedure described below is performed.

####  5.1.2. Numerical simulation

Given initial conditions $u_{0}$ and $v_{0}$, we numerically simulate the
solution $w$ over a given time period. We use a basic two step explicit finite
difference time domain (FDTD) numerical scheme for the wave equation as
described in [10], equation A.24, over the cube $[0,1]^{3}$. The space unit is
the meter (each side is $L=1$ meter), the time unit is the second. We also use
second order Engquist-Majda transparent boundary conditions [24]. The
simulation parameters are as follows. We use a sample rate $SR=200\ Hz$ (time
step $\Delta t=1/200$ s), a space step $\Delta x=43\ mm$, a wave speed $c=0.5\
m/s$. The simulation duration is $T=1.5\ s=0.75\ L/c$.

####  5.1.3. Database generation

30 sensors are scattered in the cube $[0.2,0.8]^{3}$ using a Latin hypercube
repartition and a minimax space filling algorithm. Signal outputs correspond
to time series for each sensor, with a sample rate of $50\ Hz$, so $75$ data
points altogether equally spanned over the interval $[0,1.5]\ s$ for each
sensor, which leads to $30\times 75=2250$ observations. Each signal is then
artificially polluted by a centered Gaussian white noise with standard
deviation $\sigma_{noise}=0.09$. This value corresponds to around $10\%$ of
the maximal amplitude of the signals, see Figures 2, 4 and 6.

####  5.1.4. WIGPR

We perform WIGPR on (portions of) the dataset obtained above, using the kergp
package [21] from R [57]. For that we use kernels (4.56) and/or (4.57) which
are “fast” to evaluate, with $K_{v}$ and $k_{u}^{0}$ both 1D $5/2-$Matérn
kernels. This Matérn kernel is stationary and writes, in term of the increment
$h=x-x^{\prime}$,

(5.1) $\displaystyle
k_{5/2}(h)=\sigma^{2}\bigg{(}1+\frac{|h|}{\rho}+\frac{h^{2}}{3\rho^{2}}\bigg{)}\exp\bigg{(}-\frac{|h|}{\rho}\bigg{)}$

It has two additional hyperparameters : $\rho$ and $\sigma^{2}$. $\rho$ is the
length scale of the kernel (5.1) and should correspond to the typical
variation length scale of the function approximated with GPR. $\sigma^{2}$ is
the variance of the kernel. We tackle two different questions related to
WIGPR.

  1. ($Q_{1}$)

We study how well the physical parameters (propagation speed $c$, source
position $x_{0}$ and size $R$) can be estimated with WIGPR. For this, we first
select $N_{s}$ time series corresponding to the first $N_{s}$ sensors with
$N_{s}\in\\{1,2,3,4,5,6,8,10,15,20,25,30\\}$. The corresponding database
encompasses $75\times N_{s}$ data points. For this database, we perform
negative log marginal likelihood minimization to estimate the corresponding
hyperparameters, which are either

$\displaystyle\theta=\begin{cases}(x_{0}^{a},R_{a},\theta_{k_{a}^{0}},c,\lambda)\in\mathbb{R}^{8},a\in\\{u,v\\},&\text{
if }v_{0}=0\hskip 1.0pt\text{ or }\hskip 1.0ptu_{0}=0\\\
(x_{0}^{u},R_{u},\theta_{k_{u}^{0}},x_{0}^{v},R_{v},\theta_{k_{v}^{0}},c,\lambda)\in\mathbb{R}^{14}&\text{
if }v_{0}\neq 0\text{ and }u_{0}\neq 0\end{cases}$

$\lambda$ corresponds to $\sigma^{2}$ in (2.9), and is viewed as an additional
hyperparameter in the log marginal likelihood. We use a COBYLA optimization
algorithm to solve (2.5) and a multistart procedure with $n_{mult}=100$
different starting points. That is, $100$ different values of $\theta_{0}$ are
scattered over an hypercube $H\subset\mathbb{R}^{8}$ or
$H\subset\mathbb{R}^{14}$, and the COBYLA log marginal likelihood optimization
procedure is run using each value of $\theta_{0}$ as a starting point. Among
all the obtained results, the hyperparameter value that provides the maximal
log marginal likelihood is selected. The multistart procedure mitigates the
risk of getting stuck in local maxima. COBYLA is a gradient-free optimization
method used in kergp and available in the nloptr package from R. We then
reconstruct the initial conditions using WIGPR.

  2. ($Q_{2}$)

We study the sensibility of the reconstruction step w.r.t. the sensor
locations. Consider $40$ different Latin hypercube layouts of the $30$
sensors, each obtained with a minimax space filling algorithm. For each
layout, we provide the correct hyperparameter values to the model; these
valued are described in each test case. We then reconstruct the initial
conditions using GPR and $N_{s}$ sensors, with
$N_{s}\in\\{1,5,10,15,20,25,30\\}$. $L^{p}$ relative errors (see equation
(5.2)) are computed between the reconstructed initial condition and the real
initial condition, and for each number of sensors $N_{s}$, statistics over the
40 different datasets for these $L^{p}$ errors are summarized in boxplots.

In both cases, the approximated initial position $\tilde{u}_{0}$ is recovered
by evaluating the WIGPR Kriging mean at $t=0$ over a 3D grid and the initial
speed $\tilde{v}_{0}$ is recovered by evaluating the Kriging mean at $t=0$ and
$t=dt=10^{-7}$ over the same 3D grid:
$\tilde{v}_{0}\simeq(\tilde{m}(\cdot,dt)-\tilde{m}(\cdot,0))/dt$ with
$dt=10^{-7}$.

####  5.1.5. Result visualization

For each test case are provided an example of captured data (observation
points), for one particular design of experiment (out of 40) for 6 sensors
(out of 30). This displayed data takes the form of a sequence of 6 time series
(Figures 2, 4 and 6). For each test case, an example of initial condition
reconstruction achieved thanks to WIGPR is also showcased. These examples take
the form of a slice of the corresponding 3D functions evaluated on some plane
$z=Cst$. Figures are obtained with MATLAB [47].

####  5.1.6. Numerical indicators

Precise numerical results are displayed in appendix A; they are as follows.

For ($Q_{1}$) are plotted the distances between the correct physical
parameters and the estimated ones, depending on the number of sensors used.
Additionally, for every $p\in\\{1,2,\infty\\}$ are plotted relative $L^{p}$
errors $e_{p,rel}$ defined below depending on the number of sensors used.

(5.2) $\displaystyle
e_{p,rel}^{u}=\frac{||u_{0}-\tilde{u}_{0}||_{p}}{||u_{0}||_{p}}\ \ \text{ and
}\ \ e_{p,rel}^{v}=\frac{||v_{0}-\tilde{v}_{0}||_{p}}{||v_{0}||_{p}}$

A relative error of $50\%$ means that
$||u_{0}-\tilde{u}_{0}||_{p}=0.5||u_{0}||_{p}$. A relative error of over
$100\%$ means that $||u_{0}-\tilde{u}_{0}||_{p}\geq||u_{0}||_{p}$, in which
case the null estimator $\hat{u}_{0}=0$ performs better than the estimator
$\tilde{u}_{0}$ provided by WIGPR. Note that obtaining a small $L^{p}$
relative error is harder when the target function has a small $L^{p}$ norm.

For ($Q_{2}$) are plotted boxplots of the relative $L^{p}$ errors over the 40
different sensor layouts, depending on the number of sensors used.

Integrals for the $L^{p}$ error plots are approximated using Riemann sums over
a 3D grid with space step $dx=0.01$.

###  5.2. Test case for $k_{\mathrm{u}}^{\mathrm{wave}}$

In this test case, $v_{0}$ is assumed null and thus we set $k_{\mathrm{v}}=0$,
which yields $k_{\mathrm{v}}^{\mathrm{wave}}=0$. We thus use
$k_{\mathrm{u}}^{\mathrm{wave}}$ defined in (4.57) for GPR. We use the 1D
Matérn kernel (5.1) for $k_{u}^{0}$ in equation (4.57). The initial condition
$u_{0}$ is a radial raised cosine described as follows. Set
$x_{0}=(0.65,0.3,0.5)^{T}$, $R=0.25$ and $A=5$, the corresponding initial
conditions (IC) are given by :

$\displaystyle\begin{cases}u_{0}(x)&=A\mathbbm{1}_{[0,R]}(|x-x_{0}|)\Bigg{(}1+\cos\bigg{(}\frac{\pi|x-x_{0}|}{R}\bigg{)}\Bigg{)}\\\
v_{0}(x)&=0\end{cases}$

See Figure 3, left column, for a graphical representation. For problem
($Q_{1}$), the optimization domain is chosen to be the following hypercube

$\displaystyle\theta$ $\displaystyle=(x_{0},R,\rho,\sigma^{2},c,\lambda)$
(5.3)
$\displaystyle\in[0,1]^{3}\times[0.03,0.5]\times[0.02,2]\times[0.1,5]\times[0.2,0.8]\times[10^{-8},10^{-2}]$

For problem ($Q_{2}$), the hyperparameter $\theta_{0}$ provided to the model
is

(5.4)
$\displaystyle\theta_{0}=(x_{0},R,(\rho,\sigma^{2}),c,\lambda)=((0.65,0.3,0.5),0.3,(0.2,3),0.5,\sigma_{noise}^{2})$

with $\sigma_{noise}^{2}=0.0081$. The value of $0.2$ provided for $\rho$ is a
visual estimation of the length scale of $u_{0}$ based on Figure 3.
Additionally, the couple $(0.2,3)$ provided for $(\rho,\sigma^{2})$ is more or
less confirmed by the estimations of $(\rho,\sigma^{2})$ performed for problem
($Q_{1}$), test case#1, visible in Table 1.

Figure 2. Example of captured signals for the test case #1 Figure 3. Test case
#1 : top and lateral view of the reconstruction of $u_{0}$ provided by WIGPR
in comparison with $u_{0}$. Left : $u_{0}$. Right : reconstruction of $u_{0}$
using WIGPR. 20 sensors were used. The images correspond to the 3D solutions
evaluated at $z=0.5$.

####  5.2.1. Discussion on the numerical results

For problem ($Q_{1}$), Figures 10(a) and 10(b) show that the physical
parameters $x_{0}$ and $c$ are well estimated. Figure 10(c) shows that the
source size parameter $R$ is overestimated, as could be expected from Remark
4.5. The relative errors in Figures 10(d), 10(e) and 10(f) show that the
overall function reconstruction is satisfying, with around $10\%$ relative
error. For problem ($Q_{2}$) (figures 14(a), 14(c) and 14(e)), the relative
errors stagnate at around $2\%$. The data summarized by the boxplots are not
very scattered, as is visible from the interquartile range (IQR) of the
boxplots. The interquartile range is the difference between the $3^{rd}$ and
the $1^{st}$ quartiles, and remains between $0.5\%$ and $2\%$ in the three
figures aforementioned. This means that for this test case, the reconstruction
step is not very sensitive to the sensors layout when they are scattered as a
Latin hypercube. For problems ($Q_{1}$) and ($Q_{2}$), the results are very
good overall. Note though that this is a simple test case : the function to be
reconstructed is differentiable and exhibits rather small variations compared
to the next test case.

###  5.3. Test case for $k_{\mathrm{v}}^{\mathrm{wave}}$

We use the kernel (4.56) and the 1D Matérn kernel (5.1) for $K_{v}$ in
equation (4.56). The initial condition $v_{0}$ is a radial ring cosine
described as follows. Set $x_{0}=(0.3,0.6,0.7)^{T}$, $R_{1}=0.05$,
$R_{2}=0.15$, $A=50$, the corresponding IC are given by :

$\displaystyle\begin{cases}u_{0}(x)&=0\\\
v_{0}(x)&=A\mathbbm{1}_{[R_{1},R_{2}]}(|x-x_{0}|)\Bigg{(}1+\cos\bigg{(}\frac{2\pi(|x-x_{0}|-\frac{R_{1}+R_{2}}{2})}{R_{2}-R_{1}}\bigg{)}\Bigg{)}\end{cases}$

See Figure 5, left column, for a graphical representation. For problem
($Q_{1}$), the optimization domain is chosen to be the one from the previous
test case defined in (5.2). For problem ($Q_{2}$), the hyperparameter
$\theta_{0}$ provided to the model is

(5.5)
$\displaystyle\theta_{0}=(x_{0},R,(\rho,\sigma^{2}),c,\lambda)=((0.3,0.6,0.7),0.15,(0.03,3),0.5,\sigma_{noise}^{2})$

with $\sigma_{noise}^{2}=0.0081$. As in the test case #1, the value of $0.03$
provided for $\rho$ is a visual estimation of the length scale of $v_{0}$
based on Figure 5. Additionally, the couple $(0.03,3)$ provided for
$(\rho,\sigma^{2})$ is more or less confirmed by the estimations of
$(\rho,\sigma^{2})$ performed for problem ($Q_{1}$), test case #2, visible in
Table 1.

Figure 4. Example of captured signals for the test case #2 Figure 5. Test case
#2 : top and lateral view of the reconstruction of $v_{0}$ provided by WIGPR
in comparison with $v_{0}$. Left : $v_{0}$. Right : reconstruction of $v_{0}$
using WIGPR. 25 sensors were used. The images correspond to the 3D solutions
evaluated at $z=0.7$.

####  5.3.1. Discussion of the numerical results

For problem ($Q_{1}$), Figures 11(a) and 11(b) show that the physical
parameters $x_{0}$ and $c$ are well estimated. As in the first test case,
Figure 11(c) shows that the source size parameter $R$ is overestimated, see
Remark 4.5. However, the relative error plots are not as good as in the
previous example. The poor $L^{1}$ relative error between $60\%$ and $100\%$
in Figure 11(f) is an indicator that the reconstructed function (WIGPR Kriging
mean) is supported on a set significantly larger than the original function,
as the $L^{1}$ norm is minimized for sparse functions. This is indeed visible
on Figure 5, and is partly due to the fact that $R$ was overestimated. The
$L^{2}$ and $L^{\infty}$ plots in Figures 11(d) and 11(e) are not too bad,
with $20\%$ to $30\%$ relative error. Although these errors are large when
compared to the test case #1, WIGPR still managed to capture the structure of
the initial condition, as is visible on Figure 5. For problem ($Q_{2}$), the
relative error plots are better than for ($Q_{1}$). The $L^{2}$ and $L^{1}$
relative errors stagnate at around $6\%$ with an IQR of around $2\%$ (figures
14(b) and 14(f)). The $L^{\infty}$ relative error stagnates at around $10\%$
with an IQR of around $5\%$ (figure 14(d)). Note finally that this test case
is more difficult than the previous one : the initial condition to reconstruct
has a smaller support and exhibits much larger variations.

###  5.4. Test case for
$k_{\mathrm{u}}^{\mathrm{wave}}+k_{\mathrm{v}}^{\mathrm{wave}}$

Here, we combine the initial conditions from the two previous test cases. Set
$x_{0}^{u}=(0.65,0.3,0.5)^{T}$, $R_{u}=0.25$, $A_{u}=2.5$,
$x_{0}^{v}=(0.3,0.6,0.7)^{T}$, $R_{1}^{v}=0.05$, $R_{2}^{v}=0.15$ and
$A_{v}=30$. The corresponding IC are given by :

$\displaystyle\begin{cases}u_{0}(x)&=A_{u}\mathbbm{1}_{[0,R_{u}]}(|x-x_{0}^{u}|)\Bigg{(}1+\cos\bigg{(}\frac{\pi|x-x_{0}^{u}|}{R_{u}}\bigg{)}\Bigg{)}\\\
v_{0}(x)&=A_{v}\mathbbm{1}_{[R_{1}^{v},R_{2}^{v}]}(|x-x_{0}^{v}|)\Bigg{(}1+\cos\bigg{(}\frac{2\pi(|x-x_{0}^{v}|-\frac{R_{1}^{v}+R_{2}^{v}}{2})}{R_{2}^{v}-R_{1}^{v}}\bigg{)}\Bigg{)}\end{cases}$

See Figures 3 and 3, left columns, for graphical representations. For problem
($Q_{1}$), the optimization domain is chosen to be the following hypercube

$\displaystyle\theta=$
$\displaystyle(x_{0}^{u},R_{u},(\rho_{u},\sigma_{u}^{2}),x_{0}^{v},R_{v},(\rho_{v},\sigma_{v}^{2}),c,\lambda)$
$\displaystyle\in$
$\displaystyle[0,1]^{3}\times[0.05,0.4]\times[0.02,2]\times[0.1,5]$ (5.6)
$\displaystyle\times$
$\displaystyle[0,1]^{3}\times[0.05,0.4]\times[0.02,2]\times[0.1,5]\times[0.2,0.8]\times[10^{-8},10^{-2}]$

For problem ($Q_{2}$), the hyperparameter value $\theta_{0}$ provided to the
model is

(5.7)
$\displaystyle\theta_{0}=((0.65,0.3,0.5),0.3,(0.2,0.3),(0.3,0.6,0.7),0.15,(0.03,3),0.5,\sigma_{noise}^{2})$

with $\sigma_{noise}^{2}=0.0081$. The provided values for
$(\rho_{u},\sigma_{u}^{2})$ and $(\rho_{v},\sigma_{v}^{2})$ are those of the
two previous test cases; see Subsections 5.2 and 5.3.

Figure 6. Example of captured signals for the test case #3 Figure 7. Test case
#3 : top and lateral view of the reconstruction of $u_{0}$ provided by WIGPR
in comparison with $u_{0}$. Left : $u_{0}$. Right : reconstruction of $u_{0}$
using WIGPR. 30 sensors were used. The images correspond to the 3D solutions
evaluated at $z=0.5$. Figure 8. Test case #3 : top and lateral view of the
reconstruction of $v_{0}$ provided by WIGPR in comparison with $v_{0}$. Left :
$v_{0}$. Right : reconstruction of $v_{0}$ using WIGPR. 30 sensors were used.
The images correspond to the 3D solutions evaluated at $z=0.7$.

####  5.4.1. Discussion of the numerical results

Figures 12(a), 12(b) and 12(e) show that the physical parameters $x_{0}^{u}$,
$x_{0}^{v}$ and $c$ are well estimated. Figures 12(c) and 12(d) show that the
source sizes $R_{u}$ and $R_{v}$ are overestimated, see Remark 4.5. For the
initial position $u_{0}$, Figures 13(a) and 13(c) show that with $N_{s}=30$
sensors, WIGPR yielded good results with around $10\%$ $L^{2}$ and
$L^{\infty}$ of relative errors. However, for the other values of $N_{s}$,
WIGPR did not perform as well with around $50\%$ of $L^{2}$ and $L^{\infty}$
relative error. As in the previous test case, the $L^{1}$ relative error plot
in Figure 13(e) shows that the reconstructed function $\tilde{u}_{0}$ is
usually supported on larger sets than the true function $u_{0}$. With $30$
sensors, WIGPR yielded $30\%$ relative error and between $60\%$ and $90\%$ of
relative $L^{1}$ error for the other values of $N_{s}$. The reconstruction of
$u_{0}$ in the case $N_{s}=30$ is displayed in Figure 7, showing a
satisfactory result. For the reconstruction of the initial speed $v_{0}$, all
the numerical attempts failed. The best result is attained for $N_{s}=30$
sensors with $65\%$ of $L^{2}$ and $L^{\infty}$ relative errors (13(b) and
13(d)). But for this case, Figure 8 shows that WIGPR failed to capture the
structure of $v_{0}$. One explanation of this is that for the test case #3,
the value of the length scale hyperparameter $\rho_{v}$ was overestimated by a
factor of at least 10 (see the $4^{th}$ row of Table 1 versus the $3^{rd}$
row). Therefore the WIGPR model was never able to capture the fast variations
of $v_{0}$. For problem ($Q_{2}$), the results are better. For the initial
position $u_{0}$, the results are very similar to those of the test case #1.
Figures 15(a), 15(c) and 15(e) show that relative errors stagnate at around
$2\%$ to $4\%$ for the three relative errors. The corresponding IQR are of the
same order of magnitude : $1\%$ to $2\%$. For the initial speed $v_{0}$
(figures 15(b), 15(d) and 15(f)), the results are very similar to those of the
test case #2 but multiplied by a factor between 3 and 6. The $L^{2}$ and
$L^{1}$ relative errors now stagnate below $30\%$ and the $L^{\infty}$
relative stagnates below $45\%$. The corresponding IQR are of around $10\%$
for the $L^{2}$ and $L^{1}$ errors and $15\%$ for the $L^{\infty}$ error. Note
finally that this test case combines the difficulties of the first and second
test cases.

Table 1. Estimated Matérn hyperparams. and noise var. for ($Q_{1}$) | Test case | $N_{s}=10$ | $N_{s}=15$ | $N_{s}=20$ | $N_{s}=25$ | $N_{s}=30$
---|---|---|---|---|---|---
$\rho_{u}$ | 1 | 0.117 | 0.11 | 0.043 | 0.184 | 0.144
3 | 1.145 | 0.357 | 0.891 | 0.529 | 0.088
$\rho_{v}$ | 2 | 0.105 | 0.149 | 0.023 | 0.068 | 0.046
3 | 1.453 | 1.703 | 0.595 | 0.482 | 0.353
$\sigma_{u}^{2}$ | 1 | 4.97 | 3.88 | 4.23 | 3.09 | 4.89
3 | 3.9 | 4.97 | 3.84 | 4.64 | 4.89
$\sigma_{v}^{2}$ | 2 | 4.25 | 5 | 1.61 | 4.17 | 4.89
3 | 3.54 | 3.05 | 3.5 | 3.27 | 3.26
$\lambda$ | 1 | 0.01 | 0.0099 | 0.0088 | 0.0086 | 0.0096
2 | 0.0097 | 0.0096 | 0.0036 | 0.0094 | 0.0095
3 | 0.0098 | 0.01 | 0.0099 | 0.0083 | 0.0094

####  5.4.2. Noise level estimation

The estimated noise variance is displayed in the last three rows of Table 1.
In all the test cases, the noise variance $\lambda$ was rather well estimated
for $N_{s}=25$ or $N_{s}=30$, relatively to the target value of
$\sigma_{noise}^{2}=0.0081$.

###  5.5. Transparent boundary conditions

The search for perfect or approximated transparent boundary conditions for
numerical schemes is an active research topic. They became a popular research
topic in 1970’ [24] and several different strategies have since been set up to
deal with them : see e.g. [11],[24], [35], or [39] for an overview of these
different methods. The need for them rises from the fact that numerical
simulation domains are necessarily finite and thus bounded; or at the very
least, discretized thanks to a finite grid. To numerically solve a free space
PDE such as (4.1), one first uses an adapted numerical scheme (finite
differences, finite elements, etc) for the desired PDE applied in the interior
domain, i.e. at the interior points of the discretization grid. At the
boundaries of the numerical domain, one is then forced to encode boundary
conditions : common examples are the Dirichlet and Von Neumann boundary
conditions, which cause reflections on the boundaries of the domain. For the
free space problem, the question is then the following : can we numerically
implement boundary conditions that produce no reflections?

Figure 9. Numerical example showing that the outgoing waves produced by WIGPR
do not suffer from numerical reflection issues.

Figure 9 shows that WIGPR provides numerical solutions for problem (4.1) that
do not suffer from numerical reflection issues. The reason is simple : these
numerical reconstructions are evaluations on a grid of functions that are
exact solutions of the free space problem (4.1). While being a promising start
for new types of transparent boundary conditions, the discussion above is not
enough to provide viable transparent boundary conditions in the sense of
numerical schemes. Given the covariance kernels from (4.9) and (4.8), two
questions should now be investigated. 1) How can WIGPR be associated to
classical numerical schemes such as finite difference schemes for the wave
equation? 2) Is the corresponding numerical transparent boundary condition
numerically stable?

##  6\. Conclusion and perspectives

In Section 3, we have presented a new result that provides a simple
characterization of the second order stochastic processes whose trajectories
verify linear differential constraints within the framework of generalized
functions. This result is one of the starting points of Section 4, in which a
Gaussian process model of the wave equation was described and thoroughly
studied. This model was derived by putting a GP prior on the wave equation’s
initial conditions, following a Bayesian approach when considering these
initial conditions to be unknown. We then showed that WIGPR naturally provides
an approximation of the initial conditions of the wave equation, as expected
with a Bayesian approach w.r.t. the initial conditions. The radial symmetry
WIGPR formulas from Section 4 were then showcased in Section 5, where two
practical questions were tackled. First, WIGPR can correctly estimate certain
physical parameters attached to the corresponding wave equation, such as wave
speed or source position. When these parameters are well estimated, the
initial condition reconstruction step is rather satisfactory in terms of
$L^{1}$, $L^{2}$ and $L^{\infty}$ relative errors. Second, we observed that
the WIGPR reconstruction step is not very sensitive to the layout of sensors
when they are spread according to a Latin hypercube, assuming that the correct
set of hyperparameters is provided to the model.

Future possible investigations may concern the practical use of the more
general formula (4.10) without any symmetry assumptions. To compute the
convolution efficiently, one may then resort to Fast Fourier transforms (FFT).
The case of the two dimensional wave equation is also of practical interest,
notably in oceanography [43], and presents many different properties than its
3D counterpart ([27], p.80). It would thus deserve a theoretical and practical
study in its own right when coupled with GPR. Moreover, proposition 2
constitutes a first step towards understanding PDE constrained stochastic
processes in an integrated sense; different functional analysis frameworks can
now be considered, for example by assuming Sobolev regularity of the
coefficients of $T$ instead of $\mathcal{C}^{k}$ regularity or by trying to
work in smaller function spaces than $\mathscr{D}^{\prime}(\mathcal{D})$ by
transferring only a part of the derivatives of the PDE to the test function. A
well known particular case is when the differential operator can be “split in
two” by this procedure : it is then alternatively viewed as a bilinear form
over some Hilbert space. This leads to the variational formulation of the PDE
and is the canonical approach to studying elliptic PDEs.

#### Acknowledgements

We would like to thank the Service Hydrographique et Océanographique de la
Marine (SHOM) for funding this work, and Rémy Baraille in particular for his
personal involvement in the project. We thank the GMM laboratory of INSA
Toulouse for providing us access to their servers on which we ran our
numerical computations.

## Appendix A : Numerical Results

##  7\. Hyperparameter estimation and related error plots

###  7.1. Test case #1 : $k_{\mathrm{u}}^{\mathrm{wave}}$

Here,
$\theta=(x_{0}^{u},R_{u},(\rho_{u},\sigma_{u}^{2}),c,\lambda)\in\mathbb{R}^{8}$.
The target value is $\theta^{*}=\big{(}(0.65,0.3,0.5),0.25,(\sim
0.15,?),0.5,0.0081\big{)}$. The values $(\rho_{u},\sigma_{u}^{2})$ are the
hyperparameters of the underlying Matérn kernel $k_{u}^{0}$. The a priori
target value of $0.15$ for $\rho_{u}$ is a visual estimation of the
characteristic length scale of $u_{0}$ from the figure 3. The corresponding
variance $\sigma_{u}^{2}$ is difficult to estimate a priori and its target
value is left at “?”.

(a) 3D distance $|x_{0}-\hat{x}_{0}|$

(b) Distance $|c-\hat{c}|$

(c) Estimated source radius

(d) $L^{2}$ rel. error

(e) $L^{\infty}$ rel. error

(f) $L^{1}$ rel. error

Figure 10. Hyperparameter estimation and error plots for the test case #1

###  7.2. Test case #2 : $k_{\mathrm{v}}^{\mathrm{wave}}$

Here,
$\theta=(x_{0}^{v},R_{v},(\rho_{v},\sigma_{v}^{2}),c,\lambda)\in\mathbb{R}^{8}$.
The target value is $\theta^{*}=\big{(}(0.3,0.6,0.7),0.15,(\sim
0.05,?),0.5,0.0081\big{)}$. The values $(\rho_{v},\sigma_{v}^{2})$ are the
hyperparameters of the underlying Matérn kernel $k_{v}^{0}$. The a priori
target value of $0.05$ for $\rho_{v}$ is a visual estimation of the
characteristic length scale of $v_{0}$ from the figure 5. The corresponding
variance $\sigma_{v}^{2}$ is difficult to estimate a priori and its target
value is left at “?”.

(a) 3D distance $|x_{0}-\hat{x}_{0}|$

(b) Distance $|c-\hat{c}|$

(c) Estimated source radius

(d) $L^{2}$ rel. error

(e) $L^{\infty}$ rel. error

(f) $L^{1}$ rel. error

Figure 11. Hyperparameter estimation and error plots for the test case #2

###  7.3. Test case #3 :
$k_{\mathrm{u}}^{\mathrm{wave}}+k_{\mathrm{v}}^{\mathrm{wave}}$

Here,
$\theta=(x_{0}^{u},R_{u},(\rho_{u},\sigma_{u}^{2}),x_{0}^{v},R_{v},(\rho_{v},\sigma_{v}^{2}),c,\lambda)\in\mathbb{R}^{14}$.
The target value is $\theta^{*}=\big{(}(0.65,0.3,0.5),0.3,(\sim
0.15,?),(0.3,0.6,0.7),0.15,(\sim 0.05,?),0.5,\lambda\big{)}$. The values
$(\rho_{u},\sigma_{u}^{2})$ and $(\rho_{v},\sigma_{v}^{2})$ are the
hyperparameters of the underlying Matérn kernels $k_{u}^{0}$ and $k_{v}^{0}$.
The a priori target values of $0.15$ for $\rho_{u}$ and $0.05$ for $\rho_{v}$
are visual estimations of the characteristic length scale of $u_{0}$ and
$v_{0}$ from the figures 3 and 5. The corresponding variances $\sigma_{u}^{2}$
and $\sigma_{v}^{2}$ are difficult to estimate a priori and the corresponding
target values are left at “?”.

(a) 3D distance $|x_{0}^{u}-\hat{x}_{0}^{u}|$

(b) 3D distance $|x_{0}^{v}-\hat{x}_{0}^{v}|$

(c) Estimated source radius for $u_{0}$

(d) Estimated source radius for $v_{0}$

(e) Distance $|c-\hat{c}|$

Figure 12. Hyperparameter estimation for the test case #3

(a) $L^{2}$ rel. error for $u_{0}$

(b) $L^{2}$ rel. error for $v_{0}$

(c) $L^{\infty}$ rel. error for $u_{0}$

(d) $L^{\infty}$ rel. error for $v_{0}$

(e) $L^{1}$ rel. error for $u_{0}$

(f) $L^{1}$ rel. error for $v_{0}$

Figure 13. Error plots for the test case #3

##  8\. Experience plan sensibility results

In all the presented figures are represented the following objects. Each box
plot shows the median, the first and the third quartiles of a dataset
corresponding to results obtained on the 40 different receiver dispositions.
The “targets” (dot inside a circle) correspond to the median of each boxplot.
The dashed line links the mean of each box plot. The black crosses are the
mean of each box plot, through which each dashed curve pass.

(a) $L^{2}$ rel. error for $u_{0}$ (test case #1)

(b) $L^{2}$ rel. error for $v_{0}$ (test case #2)

(c) $L^{\infty}$ rel. error for $u_{0}$ (test case #1)

(d) $L^{\infty}$ rel. error for $v_{0}$ (test case #2)

(e) $L^{1}$ rel. error for $u_{0}$ (test case #1)

(f) $L^{1}$ rel. error for $v_{0}$ (test case #2)

Figure 14. Error plots for the test cases #1 and #2

(a) $L^{2}$ rel. error for $u_{0}$ (test case #3)

(b) $L^{2}$ rel. error for $v_{0}$ (test case #3)

(c) $L^{\infty}$ rel. error for $u_{0}$ (test case #3)

(d) $L^{\infty}$ rel. error for $v_{0}$ (test case #3)

(e) $L^{1}$ rel. error for $u_{0}$ (test case #3)

(f) $L^{1}$ rel. error for $v_{0}$ (test case #3)

Figure 15. Error plots for the test case #3

## References

  * [1] R. J. Adler and J. E. Taylor, _Random fields and geometry_ , Springer-Verlag New York, NY, 2007.
  * [2] C. G. Albert and K. Rath, _Gaussian process regression for data fulfilling linear differential equations with localized sources_ , Entropy 22 (2020), no. 2.
  * [3] P. A. Alvarado, M. A. Alvarez, G. Daza-Santacoloma, A. Orozco, and G. Castellanos-Dominguez, _A latent force model for describing electric propagation in deep brain stimulation: A simulation study_ , 2014 36th Annual International Conference of the IEEE Engineering in Medicine and Biology Society, 2014, pp. 2617–2620.
  * [4] M. Álvarez, D. Luengo, and N. Lawrence, _Linear latent force models using Gaussian processes_ , IEEE Transactions on Pattern Analysis and Machine Intelligence 35 (2013), 2693–2705.
  * [5] W. F. Ames, _Numerical methods for partial differential equations (second edition)_ , Academic Press, second edition ed., 1977.
  * [6] Habib Ammari (ed.), _Mathematical modeling in biomedical imaging II. optical, ultrasound, and opto-acoustic tomographies_ , Lecture Notes in Mathematics, vol. 2035, Springer, Berlin, Heidelberg, 2012.
  * [7] L. Arkhipova, _Number of lattice points in a sphere_ , Moscow University Mathematics Bulletin 63 (2008), 214–215.
  * [8] Jean-Marc Azaïs and Mario Wschebor, _Level sets and extrema of random processes and fields_ , Wiley & Sons, 2009.
  * [9] A. Berlinet and C. Thomas-Agnan, _Reproducing kernel Hilbert spaces in probability and statistics_ , Springer US, 2004.
  * [10] S. Bilbao, _Wave and scattering methods for numerical simulations_ , John Wiley & Sons, Ltd, 2004.
  * [11] Jean-Pierre Bérenger, _Perfectly matched layer (PML) for computational electromagnetics_ , Synthesis Lectures on Computational Electromagnetics 2 (2007), no. 1, 1–117.
  * [12] Martin D. Buhmann, _Radial basis functions: Theory and implementations_ , Cambridge Monographs on Applied and Computational Mathematics, Cambridge University Press, 2003.
  * [13] Yifan Chen, Bamdad Hosseini, Houman Owhadi, and Andrew M. Stuart, _Solving and learning nonlinear PDEs with Gaussian processes_ , Journal of Computational Physics 447 (2021), 110668.
  * [14] Igor Cialenco, Gregory E Fasshauer, and Qi Ye, _Approximation of stochastic partial differential equations by a kernel-based collocation method_ , International Journal of Computer Mathematics 89 (2012), no. 18, 2543–2561.
  * [15] J. Cockayne, C. Oates, T. Sullivan, and M. Girolami, _Probabilistic numerical methods for partial differential equations and bayesian inverse problems_ , ArXiv (2017).
  * [16] Jon Cockayne, Chris Oates, Tim Sullivan, and Mark Girolami, _Probabilistic numerical methods for pde-constrained bayesian inverse problems_ , AIP Conference Proceedings 1853 (2017), no. 1, 060001.
  * [17] Patrick R Conrad, Mark Girolami, Simo Särkkä, Andrew Stuart, and Konstantinos Zygalakis, _Statistical analysis of differential equations: introducing probability measures on numerical solutions_ , Statistics and Computing 27 (2017), no. 4, 1065–1082.
  * [18] S. L. Cotter, M. Dashti, and A. M. Stuart, _Approximation of Bayesian inverse problems for PDEs_ , SIAM Journal on Numerical Analysis 48 (2010), no. 1, 322–345.
  * [19] Robert C Dalang and Marta Sanz-Solé Sol, _Hölder-Sobolev regularity of the solution to the stochastic wave equation in dimension three_ , vol. 199, Memoirs of the American Mathematical Society, no. 931, American Mathematical Society, 2009.
  * [20] Masoumeh Dashti and Andrew M. Stuart, _The Bayesian approach to inverse problems_ , Handbook of Uncertainty Quantification (Cham) (Roger Ghanem, David Higdon, and Houman Owhadi, eds.), Springer International Publishing, 2017, pp. 311–428.
  * [21] Y. Deville, D. Ginsbourger, and O. Roustant. Contributors: Nicolas Durrande., _kergp: Gaussian process laboratory_ , 2021, R package version 0.5.5.
  * [22] R. M. Dudley, _Real analysis and probability_ , 2 ed., Cambridge Studies in Advanced Mathematics, Cambridge University Press, 2002.
  * [23] Dean G Duffy, _Green’s functions with applications_ , second edition ed., Chapman and Hall/CRC, 2015.
  * [24] Bjorn Engquist and Andrew Majda, _Absorbing boundary conditions for the numerical simulation of waves_ , Mathematics of Computation 31 (1977), no. 139, 629–651.
  * [25] S. Esmailzadeh, Z. Medina-Cetina, J.W. Kang, and L.F. Kallivokas, _Varying dimensional Bayesian acoustic waveform inversion for 1d semi-infinite heterogeneous media_ , Probabilistic Engineering Mechanics 39 (2015), 56–68.
  * [26] Lawrence C Evans and Ronald F Garzepy, _Measure theory and fine properties of functions, revised edition (1st ed.)_ , Chapman and Hall/CRC, 2015.
  * [27] L.C. Evans, _Partial differential equations_ , Graduate studies in mathematics, American Mathematical Society, 1998.
  * [28] G. Fasshauer, _Meshfree approximation methods with MATLAB_ , Interdisciplinary Mathematical Sciences, 2007.
  * [29] Gregory E Fasshauer and Qi Ye, _Reproducing kernels of generalized Sobolev spaces via a Green function approach with distributional operators_ , Numerische Mathematik 119 (2011), no. 3, 585–611.
  * [30] by same author, _Reproducing kernels of Sobolev spaces via a Green kernel approach with differential operators and boundary operators_ , Advances in Computational Mathematics 38 (2013), no. 4, 891–921.
  * [31] D. Ginsbourger, O. Roustant, and N. Durrande, _On degeneracy and invariances of random fields paths with applications in Gaussian process modelling_ , Journal of Statistical Planning and Inference (2016), 170 :117 – 128.
  * [32] T. Graepel, _Solving Noisy Linear Operator Equations by Gaussian Processes: Application to Ordinary and Partial Differential Equations_ , Proceedings of the 20th International Conference on Machine Learning, AAAI Press, 2003, pp. 234–241.
  * [33] M. Gulian, A. Frankel, and L. Swiler, _Gaussian process regression constrained by boundary value problems_ , Computer Methods in Applied Mechanics and Engineering 388 (2022), 114117.
  * [34] Lars Hörmander, _The analysis of linear partial differential operators I: Distribution theory and Fourier analysis_ , Springer, 2015.
  * [35] Fang Q. Hu, _Absorbing boundary conditions_ , International Journal of Computational Fluid Dynamics 18 (2004), no. 6, 513–522.
  * [36] Svante Janson, _Gaussian Hilbert spaces_ , Cambridge Tracts in Mathematics, Cambridge University Press, 1997.
  * [37] C. Jidling, J. Hendriks, N. Wahlstrom, A. Gregg, T. Schon, C. Wensrich, and A. Wills, _Probabilistic modelling and reconstruction of strain_ , Nuclear Instruments & Methods in Physics Research Section B-beam Interactions With Materials and Atoms 436 (2018), 141–155.
  * [38] C. Jidling, N. Wahlström, A. Wills, and T. B. Schön, _Linearly constrained Gaussian processes_ , Advances in Neural Information Processing Systems (I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, eds.), vol. 30, Curran Associates, Inc., 2017\.
  * [39] Maria Kazakova and Pascal Noble, _Discrete transparent boundary conditions for the linearized Green–Naghdi system of equations_ , SIAM Journal on Numerical Analysis 58 (2020), no. 1, 657–683.
  * [40] P. Kuchment and L. Kunyansky, _Mathematics of photoacoustic and thermoacoustic tomography_ , Handbook of Mathematical Methods in Imaging v. 2, Ch. 19 (2010), 817–866.
  * [41] M. Lange-Hegermann, _Algorithmic linearly constrained Gaussian processes_ , Advances in Neural Information Processing Systems (S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, eds.), vol. 31, Curran Associates, Inc., 2018.
  * [42] by same author, _Linearly constrained Gaussian processes with boundary conditions_ , Proceedings of The 24th International Conference on Artificial Intelligence and Statistics (A. Banerjee and K. Fukumizu, eds.), Proceedings of Machine Learning Research, vol. 130, PMLR, 13–15 Apr 2021, pp. 1090–1098.
  * [43] David Lannes and Philippe Bonneton, _Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation_ , Physics of Fluids 21 (2009), no. 1, 016601.
  * [44] Jean-François Le Gall, _Mouvement brownien, martingales et calcul stochastique_ , Springer, 2013.
  * [45] A. F. López-Lopera, N. Durrande, and M. Álvarez, _Physically-inspired Gaussian process models for post-transcriptional regulation in drosophila_ , IEEE/ACM Transactions on Computational Biology and Bioinformatics 18 (2021), 656–666.
  * [46] M. Álvarez, D. Luengo, and N. D. Lawrence, _Latent force models_ , Proceedings of the Twelth International Conference on Artificial Intelligence and Statistics (Hilton Clearwater Beach Resort, Clearwater Beach, Florida USA) (D. van Dyk and M. Welling, eds.), Proceedings of Machine Learning Research, vol. 5, PMLR, 16–18 Apr 2009, pp. 9–16.
  * [47] The Mathworks, Inc., Natick, Massachusetts, _MATLAB version 9.8.0.1721703 (R2020a) Update 7_ , 2020.
  * [48] Fábio Macêdo Mendes and Edson Alves da Costa Júnior, _Bayesian inference in the numerical solution of Laplace’s equation_ , AIP Conference Proceedings 1443 (2012), no. 1, 72–79.
  * [49] F. J. Narcowich and J. Ward, _Generalized Hermite interpolation via matrix-valued conditionally positive definite functions_ , Mathematics of Computation 63 (1994), 661–687.
  * [50] N.C. Nguyen and J. Peraire, _Gaussian functional regression for linear partial differential equations_ , Computer Methods in Applied Mechanics and Engineering 287 (2015), 69–89.
  * [51] A. P. Nikitin, ST John, A. Solin, and S. Kaski, _Non-separable spatio-temporal graph kernels via SPDEs_ , ArXiv (2021).
  * [52] Houman Owhadi, _Bayesian numerical homogenization_ , Multiscale Modeling & Simulation 13 (2015), no. 3, 812–828.
  * [53] M. Álvarez P. A. Alvarado and A. Orozco, _A three spatial dimension wave latent force model for describing excitation sources and electric potentials produced by deep brain stimulation_ , arXiv (2016).
  * [54] A. Pazy, _Semigroups of linear operators and applications to partial differential equations_ , Applied Mathematical Sciences, 1983.
  * [55] William H Press, Saul A Teukolsky, William T Vetterling, and Brian P Flannery, _Numerical recipes 3rd edition: The art of scientific computing_ , Cambridge university press, 2007.
  * [56] Zenith Purisha, Carl Jidling, Niklas Wahlström, Thomas B Schön, and Simo Särkkä, _Probabilistic approach to limited-data computed tomography reconstruction_ , Inverse Problems 35 (2019), no. 10, 105004.
  * [57] R Core Team, _R: A language and environment for statistical computing_ , R Foundation for Statistical Computing, Vienna, Austria, 2020.
  * [58] M. Raissi, P. Perdikaris, and G. E. Karniadakis, _Machine learning of linear differential equations using Gaussian processes_ , Journal of Computational Physics 348 (2017), 683–693.
  * [59] by same author, _Numerical Gaussian processes for time-dependent and nonlinear partial differential equations_ , SIAM Journal on Scientific Computing 40 (2018), no. 1, A172–A198.
  * [60] C. E. Rasmussen and C.K.I. Williams, _Gaussian processes for machine learning_ , the MIT Press, 2006.
  * [61] K. Ritter, _Average-case analysis of numerical problems_ , Lecture Notes in Mathematics, Springer Berlin Heidelberg, 2007.
  * [62] W. Rudin, _Functional analysis_ , McGraw-Hill, New York, 1991.
  * [63] S. Särkkä, M. Álvarez, and N. Lawrence, _Gaussian process latent force models for learning and stochastic control of physical systems_ , IEEE Transactions on Automatic Control 64 (2019), 2953–2960.
  * [64] Simo Särkkä, _Linear operators and stochastic partial differential equations in Gaussian process regression_ , Artificial Neural Networks and Machine Learning – ICANN 2011 (Berlin, Heidelberg) (Timo Honkela, Włodzisław Duch, Mark Girolami, and Samuel Kaski, eds.), Springer Berlin Heidelberg, 2011, pp. 151–158.
  * [65] R. Schaback, _Solving the Laplace equation by meshless collocation using harmonic kernels_ , Advances in Computational Mathematics 31 (2009), 457–470.
  * [66] H.H. Schaefer, _Topological vector spaces_ , Graduate Texts in Mathematics, Springer-Verlag New York, 1999.
  * [67] M. Scheuerer and M. Schlather, _Covariance models for divergence-free and curl-free random vector fields_ , Stochastic Models 28 (2012), 433 – 451\.
  * [68] Michael Scheuerer, _Regularity of the sample paths of a general second order random field_ , Stochastic Processes and their Applications 120 (2010), no. 10, 1879–1897.
  * [69] D. Serre, _Systems of conservation laws_ , Cambridge University Press, 1999\.
  * [70] S.Lang, _Real and functional analysis_ , Graduate Texts in Mathematics, Springer-Verlag New York, 1993.
  * [71] A. Solin and Simo Särkkä, _Hilbert space methods for reduced-rank Gaussian process regression_ , Statistics and Computing 30 (2020), 419–446.
  * [72] Arno Solin, _Stochastic Differential Equation Methods for Spatio-Temporal Gaussian Process Regression_ , Doctoral thesis, School of Science, 2016, pp. 68 + app. 72.
  * [73] A. M. Stuart, _Inverse problems: A Bayesian perspective_ , Acta Numerica 19 (2010), 451–559.
  * [74] Laura P. Swiler, Mamikon Gulian, Ari L. Frankel, Cosmin Safta, and John D. Jakeman, _A survey of constrained Gaussian process regression: Approaches and implementation challenges_ , Journal of Machine Learning for Modeling and Computing 1 (2020), no. 2, 119–156.
  * [75] A. Tarantola, _Inverse problem theory and methods for model parameter estimation_ , SIAM, 2005.
  * [76] F. Treves, _Topological vector spaces, distributions and kernels_ , Dover books on mathematics, Dover Publications, 2006.
  * [77] Ricardo Carrizo Vergara, Denis Allard, and Nicolas Desassis, _A general framework for SPDE-based stationary random fields_ , Bernoulli 28 (2022), no. 1, 1–32.
  * [78] G. Vidossich, _Characterization of separability for LF-spaces_ , Annales de l’institut Fourier tome 18 (1968), no. 2, p. 87–90.
  * [79] N. Wahlstrom, M. Kok, T. B. Schön, and F. Gustafsson, _Modeling magnetic fields using Gaussian processes_ , 2013 IEEE International Conference on Acoustics, Speech and Signal Processing (2013), 3522–3526.
  * [80] Holger Wendland, _Scattered data approximation_ , vol. 17, Cambridge university press, 2004.


This work has been funded by the German Research Foundation (DFG) within the Collaborative Research Center (CRC) 1053 MAKI.

[1]
authorEdoardo Amaldi,
authorStefano Coniglio,
authorArie M.C.A. Koster &
authorMartin Tieves
(year2016): titleOn the computational complexity
of the virtual network embedding problem.
journalElectronic Notes in Discrete
Mathematics, pp. pages213–220,

[2]
authorAnthony Anjorin,
authorNils Weidmann,
authorRobin Oppermann,
authorLars Fritsche &
authorAndy Schürr
(year2020): titleAutomating Test Schedule
Generation with Domain-Specific Languages: A Configurable, Model-Driven
In: booktitleProc. of the Int. Conf. on Model
Driven Engineering Languages and Systems, seriesMODELS '20,
publisherACM, p. pages320–331,

[3]
authorStephen P. Bradley,
authorArnoldo C. Hax &
authorThomas L. Magnanti
(year1977): titleApplied Mathematical

[4]
authorMichael R. Bussieck &
authorAlex Meeraus
(year2004): titleGeneral Algebraic Modeling
System (GAMS).
publisherSpringer, 10.1007/978-1-4613-0215-5_8.

[5]
authorMartin Fleck,
authorJavier Troya &
authorManuel Wimmer
(year2016): titleSearch-based model
journalJournal of Software: Evolution and
Process, pp. pages1081–1117, 10.1002/smr.1804.

[6]
authorCharles L. Forgy
(year1982): titleRete: A Fast Algorithm for the
Many Pattern/Many Object Pattern Match Problem.
journalArtificial Intelligence, p.
pages17–37, 10.1016/0004-3702(82)90020-0.

[7]
authorRobert Fourer,
authorDavid M. Gay &
authorBrian W. Kernighan
(year2002): titleAMPL: A Modeling Language for
Mathematical Programming.
publisherCengage Learning, 10.1287/mnsc.36.5.519.

[8]
authorSebastian Götz,
authorJohannes Mey,
authorRené Schöne &
authorUwe Aßmann
(year2018): titleA JastAdd- and ILP-based
Solution to the Software-Selection and Hardware-Mapping-Problem at the TTC
In: booktitleProc. of Transformation Tool
Contest, seriesTTC@STAF '18, publisherCEUR-WS.org,
pp. pages31–36.

[9]
authorSebastian Götz,
authorJohannes Mey,
authorRené Schöne &
authorUwe Aßmann
(year2018): titleQuality-based
Software-Selection and Hardware-Mapping as Model Transformation Problem.
In: booktitleProc. of Transformation Tool
Contest, seriesTTC@STAF '18, publisherCEUR-WS.org,
pp. pages3–11.

[10]
authorCarl Hewitt, authorPeter
Bishop & authorRichard Steiger (year1973): titleA Universal
Modular ACTOR Formalism for Artificial Intelligence.
In: booktitleProc. of the Int. joint Conf. on
Artificial Intelligence, seriesIJCAI ’73,
publisherMorgan Kaufmann Publishers Inc., p.
pages235–245, 10.5555/1624775.1624804.

[11]
authorThorsten Koch
(year2006): titleRapid Mathematical Programming
or How to Solve Sudoku Puzzles in a Few Seconds.
In: booktitleOperations Research Proceedings,
seriesGOR '05, publisherSpringer, pp.
pages21–26, 10.1007/3-540-32539-5_4.

[12]
authorErhan Leblebici,
authorAnthony Anjorin &
authorAndy Schürr
(year2017): titleInter-model Consistency
Checking Using Triple Graph Grammars and Linear Optimization Techniques.
In: booktitleProc. of the Int. Conf. on
Fundamental Approaches to Software Engineering, seriesFASE '17,
publisherSpringer, pp. pages191–207,

[13]
authorDavid G. Luenberger &
authorYinyu Ye (year1984):
titleLinear and Nonlinear Programming.
publisherSpringer, 10.1007/978-3-319-18842-3.

[14]
authorUwe Pohlmann &
authorMarcus Hüwe
(year2019): titleModel-driven allocation
engineering: specifying and solving constraints based on the example of
automotive systems.
journalAutomated Software Engineering, pp.
pages315–378, 10.1007/s10515-018-0248-3.

[15]
authorRamesh Raman &
authorIgnacio E. Grossmann
(year1994): titleModelling and computational
techniques for logic based integer programming.
journalComputers & Chemical Engineering, pp.
pages563–578, 10.1016/0098-1354(93)E0010-7.

[16]
authorSiqi Shen, authorVincent
Van Beek & authorAlexandru
Iosup (year2015):
titleStatistical Characterization of Business-Critical
Workloads Hosted in Cloud Datacenters.
In: booktitleProc. of the Int. Symposium on
Cluster Computing and the Grid, seriesCCGrid '15,
publisherACM, pp. pages465–474,

[17]
authorStefan Tomaszek,
authorErhan Leblebici,
authorLin Wang & authorAndy
Schürr (year2018):
titleVirtual Network Embedding: Reducing the Search Space by
Model Transformation Techniques.
In: booktitleProc. of the Int. Conf. on Theory
and Practice of Model Transformation, seriesICMT '18,
publisherSpringer, pp. pages59–75,

[18]
authorStefan Tomaszek,
authorErhan Leblebici,
authorLin Wang & authorAndy
Schürr (year2018):
titleModel-driven Development of Virtual Network Embedding
Algorithms with Model Transformation and Linear Optimization Techniques.
In: booktitleModellierung 2018,
publisherGesellschaft für Informatik e.V., pp.

[19]
authorStefan Tomaszek,
authorRoland Speith &
authorAndy Schürr
(year2021): titleVirtual network embedding:
ensuring correctness and optimality by construction using model
transformation and integer linear programming techniques.
journalSoftware and Systems Modeling, pp.
pages1299–1332, 10.1007/s10270-020-00852-z.

[20]
authorDániel Varró,
authorGábor Bergmann,
authorÁbel Hegedüs,
authorÁkos Horváth,
authorIstván Ráth &
authorZoltán Ujhelyi
(year2016): titleRoad to a Reactive and
Incremental Model Transformation Platform: Three Generations of the VIATRA
journalSoftware and Systems Modeling, p.
pages609–629, 10.1007/s10270-016-0530-4.

[21]
authorGergely Varró &
authorFrederik Deckwerth
(year2013): titleA Rete Network Construction
Algorithm for Incremental Pattern Matching.
In: booktitleProc. of the Int. Conf. on Theory
and Practice of Model Transformations, seriesICMT ’13,
publisherSpringer, pp. pages125–140,

[22]
authorMarkus Weckesser,
authorMalte Lochau,
authorMichael Ries &
authorAndy Schürr
(year2018): titleMathematical Programming for
Anomaly Analysis of Clafer Models.
In: booktitleProc. of the Int. Conf. on Model
Driven Engineering Languages and Systems, seriesMODELS '18,
publisherACM, p. pages34–44,


Local-basis Difference Potentials Method for elliptic PDEs in complex geometry

Qing Xia

Department of Mathematics, KTH Royal Institute of Technology, Stockholm, Sweden

We develop efficient and high-order accurate finite difference methods for elliptic partial differential equations in complex geometry in the Difference Potentials framework. The main novelty of the developed schemes is the use of local basis functions defined at near-boundary grid points. The use of local basis functions allow unified numerical treatment of (i) explicitly and implicitly defined geometry; (ii) geometry of more complicated shapes, such as those with corners, multi-connected domain, etc; and (iii) different types of boundary conditions. This geometrically flexible approach is complementary to the classical difference potentials method using global basis functions, especially in the case where a large number of global basis functions are needed to resolve the boundary, or where the optimal global basis functions are difficult to obtain. Fast Poisson solvers based on FFT are employed for standard centered finite difference stencils regardless of the designed order of accuracy. Proofs of convergence of difference potentials in maximum norm are outlined both theoretically and numerically.

difference potentials, fast Poisson solver, geometric flexibility, high-order finite difference, local basis functions

§ ACKNOWLEDGEMENT

The author would like to thank COMSOL for funding. The author is grateful to Jeffrey Banks, Elias Jarlebrin, Jennifer Ryan, Anna-Karin Tornberg and the group, Siyang Wang and Sara Zahedi for many helpful discussions.