# Bayesian Modeling with Spatial Curvature Processes

Aritra Haldera, Sudipto Banerjeeb and Dipak K. Deyc
aDepartment of Biostatistics,
Drexel University, Philadelphia, PA, USA.
bDepartment of Biostatistics,
University of California, Los Angeles, CA, USA.
cDepartment of Statistics,
University of Connecticut, Storrs, CT, USA

###### Abstract

Spatial process models are widely used for modeling point-referenced variables
arising from diverse scientific domains. Analyzing the resulting random
surface provides deeper insights into the nature of latent dependence within
the studied response. We develop Bayesian modeling and inference for rapid
changes on the response surface to assess directional curvature along a given
trajectory. Such trajectories or curves of rapid change, often referred to as
_wombling_ boundaries, occur in geographic space in the form of rivers in a
flood plain, roads, mountains or plateaus or other topographic features
leading to high gradients on the response surface. We demonstrate fully model
based Bayesian inference on directional curvature processes to analyze
differential behavior in responses along wombling boundaries. We illustrate
our methodology with a number of simulated experiments followed by multiple
applications featuring the Boston Housing data; Meuse river data; and
temperature data from the Northeastern United States.

Keywords— Bayesian modeling, Directional Curvature, Gaussian Processes,
Wombling.

##  1 Introduction

Spatial data science manifests in a variety of domains including environmental
and geographical information systems (GIS) (Webster & Oliver, 2007; Burrough
et al., 2015; Schabenberger & Gotway, 2017; Plant, 2018), digital cartography
and terrain modeling (Law et al., 2000; Santner et al., 2003; Jones, 2014;
Vaughan, 2018), imaging (Winkler, 2003; Chiu et al., 2013; Dryden & Mardia,
2016), spatial econometrics and land use (LeSage & Pace, 2009), public health
and epidemiology (Elliot et al., 2000; Waller & Gotway, 2004; Lawson, 2013)
and public policy (Haining, 1993; Wise & Craglia, 2007). Spatial data analysis
seeks to estimate an underlying spatial surface representing the process
generating the data. Specific inferential interest resides with local features
of the surface including rates of change of the process at points and along
“spatial boundaries” to understand the behavior of the underlying process and
identify lurking explanatory variables or risk factors. This exercise is often
referred to as “wombling”, named after a seminal paper by Womble (1951); (also
see Gleyze et al., 2001). For regionally aggregated data, it identifies
boundaries delineating neighboring regions and has been used to study health
disparities (Lu & Carlin, 2005; Li et al., 2015; Gao et al., 2022) and
ecological boundaries (Fitzpatrick et al., 2010). For point-referenced data,
where variables are mapped at locations within an Euclidean coordinate frame
with a sufficiently smooth spatial surface, it refers to estimating spatial
gradients and identifying boundaries representing large gradients (Banerjee et
al., 2003; Banerjee & Gelfand, 2006; Qu et al., 2021).

Our current contribution develops Bayesian inference for spatial curvature
along curves on Euclidean domains. Modeling curvature will require smoothness
considerations of the process (Adler, 1981; Kent, 1989; Stein, 1999; Banerjee
& Gelfand, 2003). Observations over a finite set of locations from these
processes cannot visually inform about smoothness. Therefore, smoothness of
the process is specified from mechanistic considerations which can be
introduced through prior specifications as needed. While Bayesian inference
for first order derivatives and directional gradients have received
considerable attention (see, e.g., Morris et al., 1993; Banerjee et al., 2003;
Majumdar et al., 2006; Liang et al., 2009; Heaton, 2014; Terres & Gelfand,
2015; Wang & Berger, 2016; Terres & Gelfand, 2016; Wang et al., 2018; Qu et
al., 2021, for inferential developments involving spatial gradients from
diverse modeling and application perspectives) such processes inform about
directional change, but do not enable inference on curvature (departure from
flatness) of the spatial surface.

Analyzing surface roughness from sampling considerations can be traced at
least as far back as Greenwood (1984). We offer full inference with
uncertainty quantification about spatial curvature at a point and average
curvature along a curve from observed data after accounting for explanatory
variables. Considering second-order finite differences we establish a valid
spatial curvature process as a limit of such finite difference processes. When
formulating directional curvature, we favor the normal direction corresponding
to a chosen curve and devise a “wombling” measure to track curvature of the
surface along the curve. We derive and exploit analytical expressions of
higher order processes to avoid numerical finite differences. The Bayesian
inferential framework delivers exact posterior inference for the above
constructs on the response as well as latent (or residual) processes.

Section 2 develops the directional curvature processes through a differential
operator. Section 3 develops the vector analytic framework for curvilinear
wombling using curvature processes. Section 4 builds a hierarchical model to
exploit the preceding distribution theory and conduct curvature analysis on
the response and the latent process. Section 5 presents detailed simulation
experiments for assessing directional gradients and curvatures. Section 6
considers applications to three different data sets: Boston housing data,
Meuse river data, and Northeastern US Temperatures (the third data is
presented in the Supplement).

##  2 Spatial Curvature Processes

Let $\\{Y(\mbox{$\bf s$}):\mbox{$\bf s$}\in{\cal
S}\subset\mbox{$\mathbb{R}^{d}$}\\}$ be a univariate weakly stationary random
field with zero mean, finite second moment and a positive definite covariance
function $K(\mbox{$\bf s$},\mbox{$\bf
s$}^{\prime})=\mathrm{Cov}\left(Y(\mbox{$\bf s$}),Y(\mbox{$\bf
s$}^{\prime})\right)$ for locations $\mbox{$\bf s$},\mbox{$\bf
s$}^{\prime}\in\mbox{$\mathbb{R}^{d}$}$. In particular, under isotropy we
assume $K(\mbox{$\bf s$},\mbox{$\bf
s$}^{\prime})=\widetilde{K}\left(||\mbox{$\bf s$}-\mbox{$\bf
s$}^{\prime}||\right)$, where $||\mbox{$\bf s$}-\mbox{$\bf s$}^{\prime}||$ is
the Euclidean distance between the locations $\mbox{$\bf s$},\mbox{$\bf
s$}^{\prime}$ (Matérn, 2013). Building upon notions of mean square smoothness
(see, e.g., Stein, 1999) at an arbitrary location $\mbox{$\bf s$}_{0}$ in
$\mathbb{R}^{d}$, we focus upon second order differentiability, $\displaystyle
Y(\mbox{$\bf s$}_{0}+h\mbox{$\bf u$})=Y(\mbox{$\bf s$}_{0})+h\mbox{$\bf
u$}^{\top}\nabla Y(\mbox{$\bf s$}_{0})+h^{2}\mbox{$\bf
u$}^{\top}\nabla^{2}Y(\mbox{$\bf s$}_{0})\mbox{$\bf u$}/2+r_{2}(\mbox{$\bf
s$}_{0},h^{2}||\mbox{$\bf u$}||)$, where $r_{2}(\mbox{$\bf
s$}_{0},h^{2}||\mbox{$\bf u$}||)/h^{2}\to 0$ as $h\to 0$ in the $L_{2}$ sense
and $\nabla$ and $\nabla^{2}$ are the gradient and Hessian operators,
respectively.

For the scalar $h$ and unit vectors $\bf u$, $\bf v$, we define $\displaystyle
Y^{(2)}_{\mbox{$\bf u$},\mbox{$\bf v$},h}(\mbox{$\bf s$}_{0})=(Y(\mbox{$\bf
s$}_{0}+h(\mbox{$\bf u$}+\mbox{$\bf v$}))-Y(\mbox{$\bf s$}_{0}+h\mbox{$\bf
u$})-Y(\mbox{$\bf s$}_{0}+h\mbox{$\bf v$})+Y(\mbox{$\bf s$}_{0}))/h^{2}$ to be
the second order finite difference processes in the directions $\bf u$, $\bf
v$ at scale $h$. Being a linear function of stationary processes it is well-
defined. Passing to limits, $D^{(2)}_{\mbox{$\bf u$},\mbox{$\bf
v$}}Y(\mbox{$\bf s$}_{0})=\lim_{h\to 0}Y^{(2)}_{\mbox{$\bf u$},\mbox{$\bf
v$},h}(\mbox{$\bf s$}_{0})$. Provided the limit exists, $D^{(2)}_{\mbox{$\bf
u$},\mbox{$\bf u$}}Y(\mbox{$\bf s$}_{0})$ is defined as the directional
curvature process. If $Y(\mbox{$\bf s$})$ is a mean square second order
differentiable process in $\mathbb{R}^{d}$ for every $\mbox{$\bf
s$}\in\mbox{$\mathbb{R}^{d}$}$ then $D^{(2)}_{\mbox{$\bf u$},\mbox{$\bf
v$}}Y(\mbox{$\bf s$})=\mbox{$\bf u$}^{\top}\nabla^{2}Y(\mbox{$\bf
s$})\mbox{$\bf v$}$ is well-defined with $D^{(2)}_{\mbox{$\bf u$},\mbox{$\bf
v$}}Y(\mbox{$\bf s$})=\lim_{h\to 0}\left(h^{2}\mbox{$\bf
u$}^{\top}\nabla^{2}Y(\mbox{$\bf s$})\mbox{$\bf
v$}+\widetilde{r_{2}}\right)/h^{2}=\mbox{$\bf u$}^{\top}\nabla^{2}Y(\mbox{$\bf
s$})\mbox{$\bf v$}$, where $\widetilde{r_{2}}=r_{2}(\mbox{$\bf
s$},h^{2}||\mbox{$\bf u$}+\mbox{$\bf v$}||)-r_{2}(\mbox{$\bf
s$},h^{2}||\mbox{$\bf u$}||)-r_{2}(\mbox{$\bf s$},h^{2}||\mbox{$\bf v$}||)$.
In practice, we need only work with computing these derivatives for an
orthonormal basis of $\mathbb{R}^{d}$, say the Euclidean canonical unit
vectors along each axis $\\{\mbox{$\bf e$}_{1},\ldots,\mbox{$\bf e$}_{d}\\}$.
If $\mbox{$\bf u$}=\sum_{i=1}^{d}u_{i}\mbox{$\bf e$}_{i}$, and $\mbox{$\bf
v$}=\sum_{i=1}^{d}v_{i}\mbox{$\bf e$}_{i}$ are arbitrary unit vectors, we can
compute $D^{(2)}_{\mbox{$\bf u$},\mbox{$\bf v$}}Y(\mbox{$\bf
s$})=\sum_{i=1}^{d}\sum_{j=1}^{d}u_{i}D^{(2)}_{\mbox{$\bf e$}_{i},\mbox{$\bf
e$}_{j}}Y(\mbox{$\bf s$})v_{j}$. The directional curvature process is linear
in the sense that $D^{(2)}_{-\mbox{$\bf u$},-\mbox{$\bf v$}}Y(\mbox{$\bf
s$})=D^{(2)}_{\mbox{$\bf u$},\mbox{$\bf v$}}Y(\mbox{$\bf s$})$,
$D^{(2)}_{\mbox{$\bf u$},-\mbox{$\bf v$}}Y(\mbox{$\bf
s$})=D^{(2)}_{-\mbox{$\bf u$},\mbox{$\bf v$}}Y(\mbox{$\bf
s$})=-D^{(2)}_{\mbox{$\bf u$},\mbox{$\bf v$}}Y(\mbox{$\bf s$})$. Since
$D^{(2)}_{\mbox{$\bf w$},\mbox{$\bf w$}}Y(\mbox{$\bf s$})=||\mbox{$\bf
w$}||^{2}D^{(2)}_{\mbox{$\bf u$},\mbox{$\bf u$}}Y(\mbox{$\bf s$})$, where
$\mbox{$\bf w$}=||\mbox{$\bf w$}||\mbox{$\bf u$}$ and $\bf u$ is a unit
direction, we henceforth only consider unit directions. First order
directional gradient processes, $D^{(1)}_{\mbox{$\bf u$}}Y(\mbox{$\bf s$})$,
are reviewed in Banerjee & Gelfand (2006) and in Section S1.1 of the
Supplement. Choosing a direction is emphasized with respect to interpreting
the directional curvature processes. Directional curvature is the change in
the normal to the surface $Y(\mbox{$\bf s$})$ at $\mbox{$\bf s$}_{0}$ when
moving along a slice of the surface in the direction $\bf w$. The associated
algebraic sign locally classifies the nature of curvature at $\mbox{$\bf
s$}_{0}$—for instance, convex or concave ellipsoids (see Stevens, 1981). A
detailed discussion, with illustration, is available in Section S2 of the
Supplement.

Since $\nabla^{2}Y(\mbox{$\bf s$})$ is a symmetric matrix, to avoid
singularities arising from duplication we modify $D^{(2)}_{\mbox{$\bf
u$},\mbox{$\bf v$}}Y(\mbox{$\bf s$})$ as follows. If $vech$ is the usual half-
vectorization operator for symmetric matrices and $\mathcal{D}_{d}$ is the
duplication matrix (Magnus & Neudecker, 1980) of order $d^{2}\times d(d+1)/2$
then, $D^{(2)}_{\mbox{$\bf u$},\mbox{$\bf v$}}Y(\mbox{$\bf s$})=\mbox{$\bf
c$}^{\top}_{\mbox{$\bf u$},\mbox{$\bf v$}}vech\left(\nabla^{2}Y(\mbox{$\bf
s$})\right)$ where $\mbox{$\bf c$}_{\mbox{$\bf u$},\mbox{$\bf
v$}}^{\top}=(\mbox{$\bf u$}\otimes\mbox{$\bf v$})^{\top}\mathcal{D}_{d}$ and
$\otimes$ is the Kronecker product for matrices. If $\mbox{$\bf
u$}=(u_{1},u_{2})^{\top},\mbox{$\bf
v$}=(v_{1},v_{2})^{\top}\in\mbox{$\mathbb{R}^{2}$}$, then $\mbox{$\bf
c$}_{\mbox{$\bf u$},\mbox{$\bf v$}}=(\mbox{$\bf u$}\otimes\mbox{$\bf
v$})^{\top}\mathcal{D}_{2}=\left(v_{1}u_{1},v_{1}u_{2}+v_{2}u_{1},v_{2}u_{2}\right)^{\top}$.
The process $vech\left(\nabla^{2}Y(\mbox{$\bf s$})\right)$ in
$\mathbb{R}^{d(d+1)/2}$ consists of the pure and mixed second order
derivatives in $\nabla^{2}Y(\mbox{$\bf s$})$. The distributions needed for
inference on directional curvature processes depend on
$vech\left(\nabla^{2}Y(\mbox{$\bf s$})\right)$ rather than
$\nabla^{2}Y(\mbox{$\bf s$})$. We refer to $(\nabla Y(\mbox{$\bf
s$})^{\top},vech(\nabla^{2}Y(\mbox{$\bf s$}))^{\top})^{\top}$ as the
differential process and $\\{\mbox{$\bf u$}^{\top}\nabla Y(\mbox{$\bf s$})$,
$\mbox{$\bf c$}_{\mbox{$\bf u$},\mbox{$\bf
u$}}^{\top}vech(\nabla^{2}Y(\mbox{$\bf s$}))\\}$ as the directional
differential processes induced by $Y(\mbox{$\bf s$})$ along $\bf u$.

Inference for differential processes requires $(Y(\mbox{$\bf s$}),\nabla
Y(\mbox{$\bf s$})^{{}^{\top}},vech(\nabla^{2}Y(\mbox{$\bf s$}))^{{}^{\top}})$
to be a valid multivariate process. Its existence is derived from the limit of
corresponding finite difference approximations, which yields the cross-
covariance matrix depending on fourth (and lower) order derivatives of $K$. We
investigate the parent and differential processes using a differential
operator ${\cal L}:\mathbb{R}^{1}\to\mathbb{R}^{m}$, $m=1+d+d(d+1)/2$, where
${\cal L}Y=\left(Y,\nabla Y^{\top},vech(\nabla^{2}Y)^{\top}\right)^{\top}$.
The resulting process ${\cal L}Y(\mbox{$\bf s$})$ is also stationary with a
zero mean and a cross-covariance matrix

$V_{{\cal L}Y}(\Delta)=\begin{pmatrix}K(\Delta)&-(\nabla
K(\Delta))^{\top}&vech(\nabla^{2}K(\Delta))^{\top}\\\ \nabla
K(\Delta)&-\nabla^{2}K(\Delta)&\nabla^{3}K(\Delta)^{\top}\\\
vech(\nabla^{2}K(\Delta))&-\nabla^{3}K(\Delta)&\nabla^{4}K(\Delta)\end{pmatrix}\;,$
(1)

where $\Delta=\mbox{$\bf s$}-\mbox{$\bf s$}^{\prime}$, $\nabla K(\Delta)$ is
the $d\times 1$ gradient, $\nabla^{2}K(\Delta)$ is the $d\times d$ Hessian,
$\nabla^{3}K(\Delta)$ is the $d(d+1)/2\times d$ matrix of third derivatives
and $\nabla^{4}K(\Delta)$ is the $d(d+1)/2\times d(d+1)/2$ matrix of fourth
order derivatives associated with $K(\Delta)$. Under isotropy, $\nabla
K(\Delta)=\frac{\nabla\widetilde{K}(||\Delta||)}{||\Delta||}\Delta$, if
$A_{0}=\left(\nabla^{2}\widetilde{K}(||\Delta||)-\frac{\nabla\widetilde{K}(||\Delta||)}{||\Delta||}\right)$
then,
$\nabla^{2}K(\Delta)=\frac{\nabla\widetilde{K}(||\Delta||)}{||\Delta||}I_{d}+A_{0}\frac{\Delta\Delta^{\top}}{||\Delta||^{2}}$,
$\nabla^{3}K(\Delta)=A_{0}\bigg{\\{}\frac{vech(I_{d})^{\top}\otimes\Delta}{||\Delta||^{2}}-3\frac{vech(\Delta\Delta^{\top})^{\top}\otimes\Delta}{||\Delta||^{4}}+\frac{1}{||\Delta||^{2}}\left(\frac{\partial
vech(\Delta\Delta^{\top})}{\partial\Delta}\right)\bigg{\\}}+\nabla^{3}\widetilde{K}(||\Delta||)\cdot\frac{vech(\Delta\Delta^{\top})^{\top}\otimes\Delta}{||\Delta||^{3}}$,
where $\widetilde{K}(\Delta)$ and its derivatives are analytically computed
for our covariance functions of interest in Section S3 of the Supplement. Let
$A_{1}=\frac{\partial\Delta\otimes vech(I_{d})^{\top}}{\partial\Delta}$,
$A_{2}=\frac{\partial\Delta\otimes
vech(\Delta\Delta^{\top})^{\top}}{\partial\Delta}$,
$A_{3}=\frac{\partial}{\partial\Delta}\left(\frac{\partial
vech(\Delta\Delta^{\top})}{\partial\Delta}\right)$ be reordered tensors
(matrices) of order $d(d+1)/2\times d(d+1)/2$ conforming to the order of
corresponding elements in $vech$. Let $A_{4}$ be the element-wise product of
$\Delta$ with $\left(\frac{\partial
vech(\Delta\Delta^{\top})}{\partial\Delta}\right)$ in the same order,
$B_{1}=vech(\Delta\Delta^{\top})vech(I_{d})^{\top}$ and
$B_{2}=vech(\Delta\Delta^{\top})vech(\Delta\Delta^{\top})^{\top}$. Then,
$\nabla^{4}K(\Delta)$ is,

$\begin{split}A_{0}\left\\{\frac{A_{1}}{||\Delta||^{2}}-3\frac{A_{2}}{||\Delta||^{4}}+\frac{A_{3}}{||\Delta||^{2}}-(1+A_{4})\left(\frac{2B_{1}}{||\Delta||^{4}}+\frac{B_{1}}{||\Delta||^{3}}\right)+3\left(\frac{4B_{2}}{||\Delta||^{6}}+\frac{B_{2}}{||\Delta||^{5}}\right)\right\\}\\\
+\nabla^{3}\widetilde{K}(||\Delta||)\Bigg{(}\frac{B_{1}}{||\Delta||^{3}}+\frac{A_{2}}{||\Delta||^{3}}+\frac{A_{4}}{||\Delta||^{3}}-6\frac{B_{2}}{||\Delta||^{5}}\Bigg{)}+\nabla^{4}\widetilde{K}(||\Delta||)\frac{B_{2}}{||\Delta||^{4}}\;.\end{split}$
(2)

The resulting multivariate differential process, ${\cal L}Y$, is stationary
but not isotropic. Evidently, for the differential operator to be well-defined
under isotropy, $\nabla^{4}K({\bf 0})$ must exist since
$var(D^{(2)}_{\mbox{$\bf u$},\mbox{$\bf u$}}Y(\mbox{$\bf
s$}))=\nabla^{4}\widetilde{K}({\bf 0})$ (analogous to results in Banerjee et
al., 2003, Section 3). The directional differential operator is defined
analogously as ${\cal L}_{\mbox{$\bf u$}}Y(\mbox{$\bf s$})$ such that ${\cal
L}_{\mbox{$\bf u$}}:\mbox{$\mathbb{R}$}\to\mbox{$\mathbb{R}^{3}$}$. If
$a_{0}=\left(1-\frac{(\mbox{$\bf
u$}^{\top}\Delta)^{2}}{||\Delta||^{2}}\right)$, then analogous to (2) the
covariance function of the directional curvature process,
$\mathrm{Cov}\left(\mbox{$\bf c$}_{\mbox{$\bf u$},\mbox{$\bf
u$}}^{\top}vech(\nabla^{2}Y(\mbox{$\bf s$})),\mbox{$\bf c$}_{\mbox{$\bf
u$},\mbox{$\bf u$}}^{\top}vech(\nabla^{2}Y(\mbox{$\bf
s$}^{\prime}))\right)=\frac{3}{||\Delta||^{2}}(5a_{0}-4)a_{0}A_{0}+\frac{6}{||\Delta||}(1-a_{0})a_{0}\nabla^{3}\widetilde{K}(||\Delta||)+(1-a_{0})^{2}\nabla^{4}\widetilde{K}(||\Delta||)$.
To characterize covariance functions that admit such processes, we turn to
spectral theory. Recall that for a positive definite function $K$ defined in
$\mathbb{R}$, Bochner’s theorem (see e.g., Williams & Rasmussen, 2006)
establishes the existence of a finite positive spectral measure ${\cal F}$ on
$\mathbb{R}$. $K$ can be expressed as the inverse Fourier transform of ${\cal
F}$, $K(t)=\int_{\mathbb{R}}e^{-i\lambda t}{\cal F}(d\lambda)$. In cases where
${\cal F}$ admits a spectral density, $K(t)=\int e^{-i\lambda
t}f(\lambda)\,d\lambda$. For $\nabla^{4}K$ to exist, a trivial extension of
the result in Wang et al. (2018) requires that $f$ possess a finite fourth
moment. Examples of covariance kernels that satisfy this condition are (a) the
squared exponential covariance kernel with $K(t)=\exp(-t^{2})$
($\sigma^{2}=\phi=1$), and $f(\lambda)=1/2\sqrt{\pi}\exp(-\lambda^{2}/4)$
then,
$\frac{1}{2\sqrt{\pi}}\int_{\mathbb{R}}\lambda^{4}\exp(-\lambda^{2}/4)\,d\lambda=3(\sqrt{2})^{4}=12$;
and (b) the Matérn class with fractal parameter, $\nu$; $f(\lambda)$ is known
to belong to the $t$-family (see e.g., Stein, 1999) with
$f(\lambda)=C(\phi,\nu)/(c(\phi,\nu)+\lambda^{2})^{\nu+1/2}$ then,
$\int_{\mathbb{R}}\lambda^{4}C(\phi,\nu)/(c(\phi,\nu)+\lambda^{2})^{\nu+1/2}\,d\lambda<\infty$,
for all $\nu>2$ (since the fourth central moment for the $t$-distribution
exists if $\nu>2$). Here, we consider formulating the directional differential
processes using these two classes of covariance functions (a) the squared
exponential,
$\widetilde{K}(||\Delta||)=\sigma^{2}\exp(-\phi||\Delta||^{\nu})$, $\nu=2$;
and (b) members of the Matérn class,
$\widetilde{K}(||\Delta||)=\sigma^{2}(\phi||\Delta||)^{\nu}K_{\nu}(\phi||\Delta||)$,
where $K_{\nu}$ is the modified Bessel function of order $\nu$ (see e.g.,
Abramowitz et al., 1988), and $\nu$ controls the smoothness of process
realizations. We are particularly interested in $\nu=5/2$.

The multivariate process, ${\cal L}Y(\mbox{$\bf s$})$, is valid under the
above assumptions without any further specific parametric assumptions over
what has been outlined above. To facilitate inference for ${\cal
L}Y(\mbox{$\bf s$})$, a probability distribution is specified for the parent
process. We assume that $Y(\mbox{$\bf s$})\sim GP(\mu(\mbox{$\bf
s$},\mbox{\boldmath${\beta}$}),K(\cdot;\sigma^{2},\phi))$ is a stationary
process specified on $\mathbb{R}^{d}$. In what follows we also assume that
$K=K(\cdot;\sigma^{2},\phi)$ admits four derivatives. There are some immediate
implications of a Gaussian assumption on the parent process. If $Y_{1}(\cdot)$
and $Y_{2}(\cdot)$ are zero mean, independent stationary Gaussian processes on
$\mathbb{R}^{d}$, then (i) the differential processes ${\cal L}Y_{1}$ and
${\cal L}Y_{2}$ are independent of each other; (ii) if
$c_{1},c_{2}\in\mbox{$\mathbb{R}$}$ are scalars, then ${\cal
L}(c_{1}Y_{1}+c_{2}Y_{2})=c_{1}{\cal L}Y_{1}+c_{2}{\cal L}Y_{2}$ is stationary
and (iii) any sub-vector of ${\cal L}Y$, for example $Y$ or $(Y,\nabla
Y^{\top})^{\top}$, is a stationary Gaussian processes.

If $K$ is $k$-times mean square differentiable (i.e. $\nabla^{2k}K$ exists),
the proposed differential operator can be extended to include higher order
derivatives of $\nabla^{k}Y(\mbox{$\bf s$})$ (Mardia et al., 1996).
Differential operators characterizing change in the response (and gradient)
surface also follow valid stationary Gaussian processes. For instance, at an
arbitrary location $\mbox{$\bf s$}_{0}$ the divergence operator, ${\rm
div}{(Y(\mbox{$\bf
s$}_{0}))}=\displaystyle{\sum_{i=1}^{d}{\frac{\partial}{\partial\mbox{$\bf
e$}_{i}}Y(\mbox{$\bf s$}_{0})}}=\mbox{$\bf c$}_{1}^{\top}{\cal L}{Y(\mbox{$\bf
s$}_{0})}$, where $\mbox{$\bf c$}_{1}$ is an $m\times 1$ vector with 0’s in
all places except for first order derivatives where it takes a value of 1, and
the Laplacian, defined as the divergence operator for gradients,
$\Delta{(Y(\mbox{$\bf
s$}_{0}))}=\displaystyle{\sum_{i=1}^{d}(\nabla^{2}{Y(\mbox{$\bf
s$}_{0})})_{ii}=\sum_{i=1}^{d}{\mfrac{\partial^{2}}{\partial\mbox{$\bf
e$}_{i}^{2}}Y(\mbox{$\bf s$}_{0})}}=\mbox{$\bf c$}_{2}^{\top}{\cal
L}{Y(\mbox{$\bf s$}_{0})}$, where $\mbox{$\bf c$}_{2}$ is a $m\times 1$ vector
with 0’s in all places except for pure second order derivatives where it takes
a value of 1. Furthermore, they follow valid Gaussian processes with ${\rm
var}({\rm div}(Y(\mbox{$\bf s$}_{0})))=\mbox{$\bf c$}_{1}^{\top}V_{\cal
L}\mbox{$\bf c$}_{1}$ and ${\rm var}(\Delta(Y(\mbox{$\bf s$}_{0})))=\mbox{$\bf
c$}_{2}^{\top}V_{\cal L}\mbox{$\bf c$}_{2}$.

Let $Y(\mbox{$\bf s$})$ be a Gaussian parent process with a twice-
differentiable mean function $\mu(\mbox{$\bf s$},\mbox{\boldmath${\beta}$})$,
i.e. $\nabla\mu(\mbox{$\bf s$},\mbox{\boldmath${\beta}$})$ and
$\nabla^{2}\mu(\mbox{$\bf s$},\mbox{\boldmath${\beta}$})$ exist, and let
$K(\cdot)$ be a covariance function with variance $\sigma^{2}$ and range
$\phi$. Let $\mbox{$\bf Y$}=(Y(\mbox{$\bf s$}_{1}),\ldots,Y(\mbox{$\bf
s$}_{L}))^{\top}$ be the observed realization over ${\cal S}$ with mean
$\mbox{\boldmath$\mu$}=(\mu(\mbox{$\bf
s$}_{1},\mbox{\boldmath${\beta}$}),\ldots,\mu(\mbox{$\bf
s$}_{L},\mbox{\boldmath${\beta}$}))^{\top}$ and $\Sigma_{\mbox{$\bf Y$}}$ be
the associated $L\times L$ covariance matrix with elements $K(\mbox{$\bf
s$}_{i},\mbox{$\bf s$}_{j})$, and $\mbox{$\bf s$}_{0}$ be an arbitrary
location. Let $\nabla\mbox{$\bf K$}_{1}=\left(\nabla
K(\delta_{1})^{\top},\ldots,\nabla K(\delta_{L})^{\top}\right)^{\top}$ and
$\nabla\mbox{$\bf
K$}_{2}=\left(vech(\nabla^{2}K(\delta_{1}))^{\top},\ldots,vech(\nabla^{2}K(\delta_{L}))^{\top}\right)^{\top}$
be $L\times d$ and $L\times d(d+1)/2$ matrices, respectively, and
$\delta_{i}=\mbox{$\bf s$}_{i}-\mbox{$\bf s$}_{0}$, $i=1,\ldots,L$. The
distribution $P(\mbox{$\bf Y$},\nabla Y(\mbox{$\bf
s$}_{0}),vech(\nabla^{2}Y(\mbox{$\bf s$}_{0}))\,|\;\mbox{\boldmath$\theta$})$,
where
$\mbox{\boldmath$\theta$}=\\{\mbox{\boldmath${\beta}$},\sigma^{2},\phi\\}$, is
the $m_{0}=L+d+d(d+1)/2$-dimensional Gaussian,

${\cal N}_{m_{0}}\left(\begin{pmatrix}\mbox{\boldmath$\mu$}\\\
\nabla\mu(\mbox{$\bf s$}_{0})\\\ vech(\nabla^{2}\mu(\mbox{$\bf
s$}_{0}))\end{pmatrix},\begin{pmatrix}\Sigma_{\mbox{$\bf
Y$}}&-\nabla\mbox{$\bf K$}_{1}&\nabla\mbox{$\bf K$}_{2}\\\ \nabla\mbox{$\bf
K$}_{1}^{\top}&-\nabla^{2}K({\bf 0})&\nabla^{3}K({\bf 0})\\\ \nabla\mbox{$\bf
K$}_{2}^{\top}&-\nabla^{3}K({\bf 0})^{\top}&\nabla^{4}K({\bf
0})\end{pmatrix}\right)\;,$ (3)

which is well-defined as long as the fourth order derivative of $K$ exists.
The posterior predictive distribution for the differential process at
$\mbox{$\bf s$}_{0}$ is

$P(\nabla Y(\mbox{$\bf s$}_{0}),vech(\nabla^{2}Y(\mbox{$\bf
s$}_{0}))\,|\;\mbox{$\bf Y$})=\int P(\nabla Y(\mbox{$\bf
s$}_{0}),vech(\nabla^{2}Y(\mbox{$\bf s$}_{0}))\,|\;\mbox{$\bf
Y$},\mbox{\boldmath$\theta$})P(\mbox{\boldmath$\theta$}\,|\;\mbox{$\bf
Y$})\,d\mbox{\boldmath$\theta$}\;.$ (4)

Posterior inference for curvature proceeds by sampling from
$P(vech(\nabla^{2}Y(\mbox{$\bf s$}_{0}))\mid\mbox{$\bf Y$})=\int
P(vech(\nabla^{2}Y(\mbox{$\bf s$}_{0}))\mid\nabla Y(\mbox{$\bf
s$}_{0}),\mbox{$\bf Y$},\mbox{\boldmath$\theta$})P(\nabla Y(\mbox{$\bf
s$}_{0})\mid\mbox{$\bf
Y$},\mbox{\boldmath$\theta$})P(\mbox{\boldmath$\theta$}\mid\mbox{$\bf
Y$})\,d\mbox{\boldmath$\theta$}\,d\nabla Y$. We sample from (4) by drawing one
instance of $(\nabla Y(\mbox{$\bf s$}_{0}),vech(\nabla^{2}Y(\mbox{$\bf
s$}_{0}))$ for each sample of $\theta$ obtained from
$P(\mbox{\boldmath$\theta$}\,|\;\mbox{$\bf Y$})$. The conditional predictive
distribution of the differential process is given by $\nabla Y(\mbox{$\bf
s$}_{0}),vech(\nabla^{2}Y(\mbox{$\bf s$}_{0}))\,|\;\mbox{$\bf
Y$},\mbox{\boldmath$\theta$}\sim{\cal
N}_{m_{1}}\left(\mbox{\boldmath$\mu$}_{1},\Sigma_{1}\right)$ where
$m_{1}=d+d(d+1)/2$, and

$\displaystyle\small\mbox{\boldmath$\mu$}_{1}$
$\displaystyle=\begin{pmatrix}\nabla\mu(\mbox{$\bf s$}_{0})\\\
vech(\nabla^{2}\mu(\mbox{$\bf
s$}_{0}))\end{pmatrix}-\begin{pmatrix}\nabla\mbox{$\bf K$}_{1}\\\
\nabla\mbox{$\bf K$}_{2}\end{pmatrix}^{\top}\Sigma_{\mbox{$\bf
Y$}}^{-1}(\mbox{$\bf Y$}-\mbox{\boldmath$\mu$})\;,$ (5)
$\displaystyle\Sigma_{1}$ $\displaystyle=\begin{pmatrix}-\nabla^{2}K({\bf
0})&\nabla^{3}K({\bf 0})^{\top}\\\ -\nabla^{3}K({\bf 0})&\nabla^{4}K({\bf
0})\end{pmatrix}-\begin{pmatrix}\nabla\mbox{$\bf K$}_{1}\\\ \nabla\mbox{$\bf
K$}_{2}\end{pmatrix}^{\top}\Sigma_{\mbox{$\bf
Y$}}^{-1}\begin{pmatrix}-\nabla\mbox{$\bf K$}_{1}\\\ \nabla\mbox{$\bf
K$}_{2}\end{pmatrix}\;.$ (6)

Analogous results follow for posterior predictive inference on the curvature
process.

If $\mu(\mbox{$\bf s$},\mbox{\boldmath${\beta}$})=\mu$ is a constant, as in
simple “kriging”, then $\nabla\mu(\mbox{$\bf s$})=\nabla^{2}\mu(\mbox{$\bf
s$})=0$. More generally, if $\mu(\mbox{$\bf
s$},\mbox{\boldmath${\beta}$})=\mbox{$\bf x$}(\mbox{$\bf
s$})^{\top}\mbox{\boldmath${\beta}$}$, where $\mbox{$\bf x$}(\mbox{$\bf s$})$
is a vector of spatially indexed covariates and $\mbox{$\bf x$}(\mbox{$\bf
s$})^{\top}\mbox{\boldmath${\beta}$}$ produces a twice differentiable trend
surface then explicit calculation of $\nabla\mu(\mbox{$\bf s$}_{0})$ and
$\nabla^{2}\mu(\mbox{$\bf s$}_{0})$ are possible. In case $Y(\mbox{$\bf
s$})=\mu(\mbox{$\bf s$},\mbox{\boldmath${\beta}$})+Z(\mbox{$\bf
s$})+\epsilon(\mbox{$\bf s$})$, where $Z(\mbox{$\bf s$})\sim GP({\bf
0},K(\cdot;\sigma^{2},\phi))$ and $\epsilon(\mbox{$\bf s$})\sim N(0,\tau^{2})$
is a white noise process, inference on gradients for the residual spatial
process, $Z(\mbox{$\bf s$})$, can be performed from the posterior predictive
distribution, $P(\nabla Z(\mbox{$\bf s$}_{0}),vech(\nabla^{2}Z(\mbox{$\bf
s$}_{0}))\,|\;\mbox{$\bf Y$})$. We address this in Section 4 in the context of
curvature wombling.

##  3 Wombling with Curvature Processes

Bayesian wombling deals with inference for line integrals

${\Gamma}(C)=\int_{C}g\left({\cal L}Y\right)\,{\rm d}\ell~{}~{}\text{
or},~{}~{}\overline{\Gamma}(C)=\frac{1}{\ell(C)}\int_{C}g\left({\cal
L}Y\right)\,{\rm d}\ell\;,$ (7)

where $C$ is a geometric structure of interest, such as lines or planar
curves, residing within the spatial domain of reference, $\ell$ is an
appropriate measure, often taken to be the arc-length measure, $g$ is a linear
function (or functional) of the differential operator ${\cal L}Y$. $\Gamma$
and $\overline{\Gamma}$ are referred to as the total and average wombling
measures respectively. The structure $C$ is defined to be a wombling boundary
if it yields a large total (or average) wombling measure. Depending on the
spatial domain, geometric structures of interest constructed within them may
vary. For example, if we are dealing with surfaces in $\mathbb{R}^{3}$,
choices of $C$ are curves and lines within the surface, with the local co-
ordinate being $\mathbb{R}^{2}$. In higher dimensions they would be planes
(curves) or hyperplanes (hypercurves). Specifically, Bayesian curvilinear
wombling involves estimating integrals in (7) over curves, which tracks rapid
change over the spatial domain by determining boundaries (curves) with large
gradients normal to the curve (see for e.g., Banerjee & Gelfand, 2006). The
integrand in (7) inherently involves a direction, in particular change
measured is always in a direction normal to $C$. Hence, $g({\cal L}Y)$ can
equivalently be expressed as a linear function (functional) of
$\displaystyle{{\cal L}_{\mbox{$\bf n$}}Y(\mbox{$\bf s$})}$, where $\mbox{$\bf
n$}=\mbox{$\bf n$}(\mbox{$\bf s$})$ denotes the unit normal vector to $C$ at
$\bf s$. The next few paragraphs provide more detail.

With wombling measures for directional gradients discussed the Supplement,
Section S1.2, we construct wombling measures for curvature. Given $C$,
depending on the smoothness of the surface, the rate at which gradients change
along the curve may present sufficient heterogeneity while traversing the
curve. If $C$ forms a wombling boundary with respect to the gradient, then
wombling boundaries for curvature are subsets of $C$ that feature segments
with large positive (negative) directional curvature along a normal direction
to the curve. Leveraging only gradients, we develop wombling measures for
curvature that further characterize such boundaries located for gradients. The
wombling measure for curvature in $Y(\mbox{$\bf s$})$ along $C$ ascertains
whether $C$ also forms a wombling boundary with respect to curvature. We
associate a directional curvature to each $\mbox{$\bf s$}\in C$,
$\displaystyle{g({\cal L}Y(\mbox{$\bf s$}))=D^{(2)}_{\mbox{$\bf n$},\mbox{$\bf
n$}}Y(\mbox{$\bf s$})=\mbox{$\bf c$}_{\mbox{$\bf n$},\mbox{$\bf
n$}}^{\top}vech(\nabla^{2}Y(\mbox{$\bf s$}))}$ (a linear function of
$\displaystyle{{\cal L}_{\mbox{$\bf n$}}Y(\mbox{$\bf s$})}$) along the
direction of a unit normal $\mbox{$\bf n$}=\mbox{$\bf n$}(\mbox{$\bf s$})$ to
$C$ at $\bf s$. Using (7) we define wombling measures for total and average
curvature as,

$\Gamma^{(2)}(C)=\int_{C}D^{(2)}_{\mbox{$\bf n$},\mbox{$\bf n$}}Y(\mbox{$\bf
s$})d\ell=\int_{C}\mbox{$\bf n$}(\mbox{$\bf s$})^{\top}\nabla^{2}Y(\mbox{$\bf
s$})\mbox{$\bf n$}(\mbox{$\bf
s$})d\ell\;,\qquad\overline{\Gamma}^{(2)}(C)=\Gamma^{(2)}(C)/\ell(C)\;,$ (8)

respectively, where $\ell(C)$ denotes the arc-length of $C$. Parameterized
curves, $C=\\{\mbox{$\bf s$}(t)=(s_{1}(t),s_{2}(t)):t\in{\cal
T}\subset\mathbb{R}\\}$, offer further insights. As $t$ varies over its
domain, $\mbox{$\bf s$}(t)$ outlines the curve $C$. Implicitly assuming that
$C$ is regular, i.e., $||\mbox{$\bf s$}^{\prime}(t)||\neq 0$, allows the
tangent and normal to exist at all points on the curve. The unit tangent and
normal at each point of the curve are $\mbox{$\bf s$}^{\prime}(t)/||\mbox{$\bf
s$}^{\prime}(t)||$ and $\mbox{$\bf n$}=\mbox{$\bf n$}(\mbox{$\bf
s$}(t))=(s^{\prime}_{2}(t),-s^{\prime}_{1}(t))^{\top}/||\mbox{$\bf
s$}^{\prime}(t)||$, respectively, while $\displaystyle\mbox{$\bf
c$}_{\mbox{$\bf n$},\mbox{$\bf n$}}=\mbox{$\bf c$}_{\mbox{$\bf n$}(\mbox{$\bf
s$}(t)),\mbox{$\bf n$}(\mbox{$\bf s$}(t))}=\left(\mbox{$\bf n$}(\mbox{$\bf
s$}(t))\otimes\mbox{$\bf n$}(\mbox{$\bf s$}(t))\right)^{\top}{\cal D}_{d}$
from Section 2.

The arc-length of $C$ is $\ell(C)=\int_{\cal T}||\mbox{$\bf
s$}^{\prime}(t)||\,dt$ or ${\rm d}\ell=||\mbox{$\bf s$}^{\prime}(t)||\,dt$. If
${\cal T}=[t_{0},t_{1}]$, then $\ell(C)=\int_{t_{0}}^{t_{1}}||\mbox{$\bf
s$}^{\prime}(t)||\,dt$ and $\Gamma^{(2)}(C)=\int_{t_{0}}^{t_{1}}\mbox{$\bf
n$}(\mbox{$\bf s$}(t))^{\top}\nabla^{2}Y(\mbox{$\bf s$}(t))\mbox{$\bf
n$}(\mbox{$\bf s$}(t))||\mbox{$\bf s$}^{\prime}(t)||\,dt$. If $C$ is an open
curve, then $\ell(C)^{-1}\int_{C}\mbox{$\bf n$}(\mbox{$\bf
s$})^{\top}\nabla^{2}Y(\mbox{$\bf s$})\mbox{$\bf n$}(\mbox{$\bf
s$})d\mbox{$\bf s$}=\ell(C)^{-1}\int_{C}\mbox{$\bf n$}(\mbox{$\bf
s$}(t))^{\top}\nabla^{2}Y(\mbox{$\bf s$}(t))\mbox{$\bf n$}(\mbox{$\bf
s$}(t))||\mbox{$\bf s$}^{\prime}(t)||\,dt$ is the average directional
curvature. For example, $C=\\{\mbox{$\bf s$}(t)=(r\cos t,r\sin
t),t\in[0,\pi/4]\\}$ is the arc of a parameterized circle of radius $r$. It
follows that $||\mbox{$\bf s$}^{\prime}(t)||=r$, $\mbox{$\bf n$}(\mbox{$\bf
s$}(t))=(\cos t,\sin t)^{\top}$ and ${\ell(C)}^{-1}\int_{0}^{\pi/4}\mbox{$\bf
n$}(\mbox{$\bf s$}(t))^{\top}\nabla^{2}Y(\mbox{$\bf s$}(t))\mbox{$\bf
n$}(\mbox{$\bf s$}(t))r\,dt=\frac{4}{\pi}\int_{0}^{\pi/4}\mbox{$\bf
n$}(\mbox{$\bf s$}(t))^{\top}\nabla^{2}Y(\mbox{$\bf s$}(t))\mbox{$\bf
n$}(\mbox{$\bf s$}(t))\,dt$. The average curvature in the tangential direction
of $C$ is $\frac{1}{\ell(C)}\int_{C}\mbox{$\bf u$}(\mbox{$\bf
s$}(t))^{\top}\nabla^{2}Y(\mbox{$\bf s$}(t))\mbox{$\bf u$}(\mbox{$\bf
s$}(t))||\mbox{$\bf
s$}^{\prime}(t)||\,dt=\displaystyle{{\ell(C)}^{-1}\int_{t_{0}}^{t_{1}}\frac{\mbox{$\bf
s$}^{\prime}(t)}{||\mbox{$\bf s$}^{\prime}(t)||}^{\top}\nabla^{2}Y(\mbox{$\bf
s$}(t))\frac{\mbox{$\bf s$}^{\prime}(t)}{||\mbox{$\bf
s$}^{\prime}(t)||}||\mbox{$\bf s$}^{\prime}(t)||\,dt}=\mbox{$\bf
u$}(\mbox{$\bf s$}(t_{1}))^{\top}\nabla Y(\mbox{$\bf s$}(t_{1}))-\mbox{$\bf
u$}(\mbox{$\bf s$}(t_{0}))^{\top}\nabla Y(\mbox{$\bf s$}(t_{0}))$. Hence, the
average directional curvature remains path independent and is the difference
of directional gradient at the end points of $C$.

For a closed curve $C$, $\oint_{C}\mbox{$\bf n$}(\mbox{$\bf
s$})^{\top}\nabla^{2}Y(\mbox{$\bf s$})\mbox{$\bf n$}(\mbox{$\bf
s$})d\mbox{$\bf s$}=\oint_{C}\mbox{$\bf n$}(\mbox{$\bf
s$}(t))^{\top}\nabla^{2}Y(\mbox{$\bf s$}(t))\mbox{$\bf n$}(\mbox{$\bf
s$}(t))||\mbox{$\bf s$}^{\prime}(t)||\,dt$. If the surface admits up to three
derivatives, i.e. $\nabla^{3}Y(\mbox{$\bf s$})$ exists, the average curvature
of the region, ${\cal D}$, enclosed by $C$, is free of $t$. If $\mbox{$\bf
F$}(\mbox{$\bf s$})=\nabla^{2}Y(\mbox{$\bf s$})=(F_{ij}(\mbox{$\bf
s$}))_{i,j=1,2}$, with $F_{12}(\mbox{$\bf s$})=F_{21}(\mbox{$\bf s$})$ and
$F_{ij}=F_{ij}(\mbox{$\bf s$})=\mfrac{\partial^{2}}{\partial s_{i}\partial
s_{j}}Y(\mbox{$\bf s$})$ then, $\oint_{C}\mbox{$\bf n$}(\mbox{$\bf
s$})^{\top}\nabla^{2}Y(\mbox{$\bf s$})\mbox{$\bf n$}(\mbox{$\bf
s$})d\mbox{$\bf s$}=\oint_{C}\mbox{$\bf n$}(\mbox{$\bf s$})^{\top}\mbox{$\bf
F$}(\mbox{$\bf s$})\mbox{$\bf n$}(\mbox{$\bf s$})d\mbox{$\bf
s$}=\oint_{C}\mbox{$\bf n$}(\mbox{$\bf s$}(t))^{\top}\mbox{$\bf F$}(\mbox{$\bf
s$})\mbox{$\bf n$}(\mbox{$\bf s$}(t))||\mbox{$\bf
s$}^{\prime}(t)||\,dt=\oint_{C}||\mbox{$\bf
s$}^{\prime}(t)||^{-1}\left(F_{11}s_{2}^{\prime}(t)^{2}-2F_{12}s_{1}^{\prime}(t)s_{2}^{\prime}(t)+F_{22}s_{1}^{\prime}(t)^{2}\right)\,dt=\iint\limits_{\cal
D}\left\\{\frac{\partial F_{11}n_{2}}{\partial s_{1}}+\left(\frac{\partial
F_{12}n_{2}}{\partial s_{2}}+\frac{\partial F_{12}n_{1}}{\partial
s_{1}}\right)+\frac{\partial F_{22}n_{1}}{\partial
s_{2}}\right\\}ds_{1}ds_{2}$. The last equality is obtained using Green’s
theorem (see for e.g., Rudin, 1976). This can be interpreted as “flux” in the
gradient within ${\cal D}$. Since, $F_{ij}(\mbox{$\bf
s$})=\nabla^{2}_{ij}Y(\mbox{$\bf s$})$, the integrand in the last equality
require the existence of $\nabla^{3}_{ijk}Y(\mbox{$\bf s$})$, $i,j,k=1,2$.
Denoting, $\widetilde{\nabla}^{3}Y(\mbox{$\bf
s$})=(\nabla^{3}_{ijk}Y(\mbox{$\bf s$}))_{i,j,k=1,2}^{\top}$, vector of unique
third derivatives, and $\mbox{$\bf n$}_{0}(\mbox{$\bf s$})=(n_{2}(\mbox{$\bf
s$}),n_{2}(\mbox{$\bf s$}),n_{1}(\mbox{$\bf s$}),n_{1}(\mbox{$\bf
s$}))^{\top}$ then,

$\frac{1}{\ell(C)}\oint_{C}\mbox{$\bf c$}_{\mbox{$\bf n$},\mbox{$\bf
n$}}^{\top}vech(\nabla^{2}Y(\mbox{$\bf s$}))\,d\mbox{$\bf
s$}=\frac{1}{\ell(C)}\iint\limits_{\cal D}\mbox{$\bf n$}_{0}(\mbox{$\bf
s$})^{\top}\widetilde{\nabla}^{3}Y(\mbox{$\bf s$})\,d\mbox{$\bf s$}.$ (9)

This extends the development in Section 3.2 of Banerjee & Gelfand (2006) to
study the behavior of spatial curvature over closed curves on surfaces in
$\mathbb{R}^{3}$. Sampling along $C$ is generally harder than sampling inside
${\cal D}$. Hence, the computational implications of (9) are more appealing.
When studying the same behavior along a tangential direction to $C$ with
$\mbox{$\bf s$}(t_{0})=\mbox{$\bf s$}(t_{1})=\mbox{$\bf s$}_{0}$,
$\displaystyle\oint_{C}\mbox{$\bf u$}(\mbox{$\bf
s$})^{\top}\nabla^{2}Y(\mbox{$\bf s$})\mbox{$\bf u$}(\mbox{$\bf
s$})d\mbox{$\bf s$}=\oint_{\mbox{$\bf s$}(t_{0})}^{\mbox{$\bf
s$}(t_{1})}F_{11}(\mbox{$\bf s$})n_{1}ds_{1}+F_{12}(\mbox{$\bf
s$})n_{1}ds_{2}+F_{21}(\mbox{$\bf s$})n_{2}ds_{1}+F_{22}(\mbox{$\bf
s$})n_{2}ds_{2}=\mbox{$\bf u$}(\mbox{$\bf s$}(t_{1}))^{\top}\nabla
Y(\mbox{$\bf s$}(t_{1}))-\mbox{$\bf u$}(\mbox{$\bf s$}(t_{0}))^{\top}\nabla
Y(\mbox{$\bf s$}(t_{0}))=0$, again a consequence of path independence. This
validates the choice of a normal direction to $C$ when measuring change in the
gradient. Using the rectilinear approximation to curvature wombling, as
discussed later, provides a more computationally tractable and simpler
approach, where double integrals manifest when computing variances of the
wombling measures.

Curvature wombling requires predictive inference performed using gradient
measures on the interval ${\cal T}$, to include $\Gamma^{(2)}({C})$ (or
$\overline{\Gamma}^{(2)}({C})$) in (8). Leveraging inference for differential
processes in Section 2, we obtain joint inference on the wombling measures.
Suppose $C=\\{\mbox{$\bf s$}(t):t\in[0,T]\\}$ is generated over ${\cal
T}=[0,T]$. For any $t^{*}\in[0,T]$, let $C_{t^{*}}$ denote the curve
restricted to $[0,t^{*}]$ and $\ell({C_{t^{*}}})$ its arc-length. Line
integrals for curvilinear gradient and curvature wombling measures are
$\Gamma^{(1)}({C_{t^{*}}})=\int_{0}^{t^{*}}D^{(1)}_{\mbox{$\bf
n$}}Y(\mbox{$\bf s$}(t))||\mbox{$\bf
s$}^{\prime}(t)||\,dt,~{}\overline{\Gamma}^{(1)}({C_{t^{*}}})=\frac{1}{\ell({C_{t^{*}}})}\Gamma^{(1)}({C_{t^{*}}})$,
$\Gamma^{(2)}({C_{t^{*}}})=\int_{0}^{t^{*}}D^{(2)}_{\mbox{$\bf n$},\mbox{$\bf
n$}}Y(\mbox{$\bf s$}(t))||\mbox{$\bf s$}^{\prime}(t)||\,dt$ and
$\overline{\Gamma}^{(2)}({C_{t^{*}}})=\frac{1}{\ell({C_{t^{*}}})}\Gamma^{(2)}({C_{t^{*}}})$.
Since $D^{(1)}_{\mbox{$\bf n$}}Y(\mbox{$\bf s$}(t))$ and $D^{(2)}_{\mbox{$\bf
n$},\mbox{$\bf n$}}Y(\mbox{$\bf s$}(t))$ are Gaussian processes on ${\cal T}$,
$\Gamma^{(1)}({C_{t^{*}}})$ and $\Gamma^{(2)}({C_{t^{*}}})$ are valid
dependent Gaussian processes on ${\cal T}$. Therefore,
$\mbox{\boldmath$\Gamma$}({C_{t^{*}}})=(\Gamma^{(1)}({C_{t^{*}}}),\Gamma^{(2)}({C_{t^{*}}}))^{\top}\sim{\cal
N}_{2}\big{(}\mbox{\boldmath$\mu$}_{\mbox{\boldmath$\Gamma$}}(t^{*}),\mbox{$\bf
K$}_{\mbox{\boldmath$\Gamma$}}(t^{*},t^{*})\big{)}$, where
$\mbox{\boldmath$\mu$}_{\mbox{\boldmath$\Gamma$}}(t^{*})=\left(\int_{0}^{t^{*}}D^{(1)}_{\mbox{$\bf
n$}}\mu(\mbox{$\bf s$}(t))||\mbox{$\bf
s$}^{\prime}(t)||\,dt\;,\int_{0}^{t^{*}}D^{(2)}_{\mbox{$\bf n$},\mbox{$\bf
n$}}\mu(\mbox{$\bf s$}(t))||\mbox{$\bf
s$}^{\prime}(t)||\,dt\right)^{\top}=(m_{1}(t^{*}),m_{2}(t^{*}))^{\top}$ and
$\mbox{$\bf
K$}_{\mbox{\boldmath$\Gamma$}}(t^{*},t^{*})=\\{k_{ij}(t^{*},t^{*})\\}_{i,j=1,2}$
whose elements are evaluated as

$k_{ij}(t^{*},t^{*})=(-1)^{j}\int_{0}^{t^{*}}\int_{0}^{t^{*}}\mbox{$\bf
a$}_{i}^{\top}(t_{1})\nabla^{i+j}K(\Delta(t_{1},t_{2}))\mbox{$\bf
a$}_{j}(t_{2})||\mbox{$\bf s$}^{\prime}(t_{1})||||\mbox{$\bf
s$}^{\prime}(t_{2})||\,dt_{1}\,dt_{2}\;,$ (10)

where $\mbox{$\bf a$}_{1}(t)=\mbox{$\bf n$}(\mbox{$\bf s$}(t))$ and
$\mbox{$\bf a$}_{2}(t)=\mbox{$\bf c$}_{\mbox{$\bf n$}(\mbox{$\bf
s$}(t)),\mbox{$\bf n$}(\mbox{$\bf s$}(t))}$. Simplifications arise in $d=2$.
For example, $\mbox{$\bf c$}_{\mbox{$\bf n$},\mbox{$\bf
n$}}(t)=(s^{\prime}_{2}(t)^{2},-2s^{\prime}_{2}(t)s^{\prime}_{1}(t),s^{\prime}_{1}(t)^{2})^{\top}$,
while $\nabla^{k}K$, for $k=2,3,4$, are matrices of orders $2\times 2$,
$2\times 3$ and $3\times 3$, respectively, of partial and mixed second, third
and fourth derivatives of $K$ and $\Delta(t_{1},t_{2})=\mbox{$\bf
s$}(t_{2})-\mbox{$\bf s$}(t_{1})$. For any two points
$t_{1}^{*},t_{2}^{*}\in{\cal T}$, the dependence is specified through
$\begin{pmatrix}\mbox{\boldmath$\Gamma$}({C_{t_{1}^{*}}})\\\
\mbox{\boldmath$\Gamma$}({C_{t_{2}^{*}}})\end{pmatrix}\sim{\cal
N}_{4}\left(\begin{pmatrix}\mbox{$\bf m$}_{1}\\\ \mbox{$\bf
m$}_{2}\end{pmatrix},\begin{pmatrix}\mbox{$\bf k$}_{11}&\mbox{$\bf k$}_{12}\\\
\ \mbox{$\bf k$}_{21}&\mbox{$\bf k$}_{22}\end{pmatrix}\right)$, where
$\mbox{$\bf m$}_{i}=(m_{i}(t_{1}^{*}),m_{i}(t_{2}^{*}))^{\top}$, $\mbox{$\bf
k$}_{ij}=\begin{pmatrix}k_{ij}(t_{1}^{*},t_{1}^{*})&k_{ij}(t_{1}^{*},t_{2}^{*})\\\
k_{ij}(t_{2}^{*},t_{1}^{*})&k_{ij}(t_{2}^{*},t_{2}^{*})\end{pmatrix}$,
$i,j=1,2$. Generally, for $n_{P}$ points partitioning ${\cal T}$ the above can
be analogously extended. Clearly, $\mbox{\boldmath$\Gamma$}({C_{t^{*}}})$ is a
mean squared continuous process. However, stationarity of $Y(\mbox{$\bf s$})$
does not imply stationarity of $\mbox{\boldmath$\Gamma$}({C_{t^{*}}})$. For
any $\mbox{$\bf s$}_{j}\in{\cal S}$ with $\mathrm{Cov}(Y(\mbox{$\bf
s$}_{j}),\mbox{\boldmath$\Gamma$}({C_{t^{*}}}))=\mbox{\boldmath$\gamma$}_{j}(t^{*})$
and $\Delta_{j}(t)=\mbox{$\bf s$}(t)-\mbox{$\bf s$}_{j}$ we have,

$\mbox{\boldmath$\gamma$}_{j}(t^{*})=\left(\int_{0}^{t^{*}}D^{(1)}_{\mbox{$\bf
n$}}K(\Delta_{j}(t))||\mbox{$\bf
s$}^{\prime}(t)||\,dt,\int_{0}^{t^{*}}D^{(2)}_{\mbox{$\bf n$},\mbox{$\bf
n$}}K(\Delta_{j}(t))||\mbox{$\bf s$}^{\prime}(t)||\,dt\right)^{\top}\;.$ (11)

A valid joint distribution can be specified over ${\cal T}$ by,

$\begin{pmatrix}\mbox{$\bf Y$}\\\
\mbox{\boldmath$\Gamma$}({C_{t^{*}}})\end{pmatrix}\sim{\cal
N}_{L+2}\left(\begin{pmatrix}\mbox{\boldmath$\mu$}\\\
\mu_{\mbox{\boldmath$\Gamma$}}(t^{*})\end{pmatrix},\begin{pmatrix}\Sigma_{\mbox{$\bf
Y$}}&\mbox{\boldmath$\gamma$}_{\mbox{\boldmath$\Gamma$}}(t^{*})\\\
\mbox{\boldmath$\gamma$}^{\top}_{\mbox{\boldmath$\Gamma$}}(t^{*})&\mbox{$\bf
K$}_{\mbox{\boldmath$\Gamma$}}(t^{*},t^{*})\end{pmatrix}\right)\;,$ (12)

where
$\mbox{\boldmath$\gamma$}^{\top}_{\mbox{\boldmath$\Gamma$}}(t^{*})=\left[\mbox{\boldmath$\gamma$}_{1}(t^{*})\;\mbox{\boldmath$\gamma$}_{2}(t^{*})\;\cdots\;\mbox{\boldmath$\gamma$}_{L}(t^{*})\right]$
is the $2\times L$ cross-covariance matrix.

In practical applications curvilinear wombling is performed by approximating
the curve $C$ using linear segments. These measures at the segment level are
then aggregated to produce a wombling measure for the curve. The curve is
segmented using a partition. Consequently, the accuracy of estimated wombling
measures for the curve depend on the choice of partition. Figures S2 and S3 in
the online Supplement illustrate this concept. Explicitly, let $C$ be a
regular rectifiable curve and $[a,b]\subset{\cal T}$ be a compact interval.
Let $g$ be a uniformly continuous function. For any partition, $P$ of $[a,b]$,
$a=t^{\prime}_{0}<t^{\prime}_{1}<\ldots<t^{\prime}_{n_{P}}=b$, with its norm
defined as
$|P|=\max\limits_{i=1,\ldots,n_{P}}(t^{\prime}_{i}-t^{\prime}_{i-1})$. A
polygonal (piecewise-linear) approximation to the curve is,
$\widetilde{C}_{P}=\bigcup\limits_{i=1}^{n_{P}}C_{t_{i}}$, where
$C_{t_{i}}=\\{\mbox{$\bf s$}(t^{\prime}_{i-1})+t\mbox{$\bf
u$}_{i},t\in[0,t_{i}]\\}$, $t_{i}=||\mbox{$\bf s$}(t^{\prime}_{i})-\mbox{$\bf
s$}(t^{\prime}_{i-1})||$ and $\mbox{$\bf u$}_{i}=||\mbox{$\bf
s$}(t^{\prime}_{i})-\mbox{$\bf s$}(t^{\prime}_{i-1})||^{-1}(\mbox{$\bf
s$}(t^{\prime}_{i})-\mbox{$\bf s$}(t^{\prime}_{i-1}))^{\top}$. Note that
$\mbox{$\bf s$}(t)=\mbox{$\bf s$}(t^{\prime}_{i-1})+t\mbox{$\bf u$}_{i}$ for
$t\in[0,t_{i}]$ and, hence, $||\mbox{$\bf s$}^{\prime}(t)||=||\mbox{$\bf
u$}_{i}||=1$. Wombling measure for $\widetilde{C}_{P}$ is,
$\Gamma(\widetilde{C}_{P})=\sum\limits_{i=1}^{n_{P}}\int_{C_{t_{i}}}g\left({\cal
L}Y(\mbox{$\bf s$}(t))\right)||\mbox{$\bf s$}^{\prime}(t)||\,dt$. As $|P|\to
0$ we have,
$\Gamma(\widetilde{C}_{P})\overset{a.s.}{\longrightarrow}\Gamma(C)=\int_{a}^{b}g\left({\cal
L}Y(\mbox{$\bf s$}(t))\right)||\mbox{$\bf s$}^{\prime}(t)||\,dt$. This
provides us with an estimate, $\Gamma(\widetilde{C}_{P})$ for curvilinear
wombling measures associated with any general curve $C$. Further details are
provided in the Supplement, at the end of Section S5.

The choices of $g$ for our wombling measures result in, $\mbox{$\bf
u$}^{\top}\nabla Y$ and $\mbox{$\bf c$}_{\mbox{$\bf u$},\mbox{$\bf
u$}}^{\top}vech(\nabla^{2}Y)$, which are linear and therefore uniformly
continuous over any compact interval. Since predictive inference is performed
iteratively on individual line segments, it is sufficient to show the
inferential procedure for an arbitrary curve segment $C_{t_{i}}$. The normal
to $C_{t_{i}}$ is free of $t$ and denoted as, $\mbox{$\bf u$}_{i}^{\perp}$,
which is the normal to $\mbox{$\bf u$}_{i}$. The associated wombling measures
with $C_{t_{i}}$ are
$\displaystyle\mbox{\boldmath$\Gamma$}(t_{i})=\left(\int_{0}^{t_{i}}D^{(1)}_{\mbox{$\bf
u$}_{i}^{\perp}}Y(\mbox{$\bf s$}(t))\,dt,\int_{0}^{t_{i}}D^{(2)}_{\mbox{$\bf
u$}_{i}^{\perp},\mbox{$\bf u$}_{i}^{\perp}}Y(\mbox{$\bf
s$}(t))\,dt\right)^{\top}$. For a point $\mbox{$\bf s$}_{j}$ define
$\Delta_{i-1,j}=\mbox{$\bf s$}_{i-1}-\mbox{$\bf s$}_{j}$, $j=1,2,\ldots,L$.
Their joint distribution is specified by (12), where
$\mbox{\boldmath$\gamma$}_{j}(t_{i})$ is obtained from (11) by replacing
$\Delta_{j}(t)$ with $\Delta_{i-1,j}+t\mbox{$\bf u$}_{i}$ and $\mbox{$\bf
K$}_{\mbox{\boldmath$\Gamma$}}(t_{i},t_{i})$ is obtained from (10) replacing
$\Delta(t_{1},t_{2})=(t_{2}-t_{1})\mbox{$\bf u$}_{i}$ in the integrand. The
analytic tractability of the line integrals in
$\mbox{\boldmath$\gamma$}_{j}(t_{i})$ is not a concern. Given choices of
$\mu(\cdot)$ and $K(\cdot)$, they are all one or two dimensional integrals
which are efficiently computed using simple quadrature. For example, let
$Y(\mbox{$\bf s$})$ be the isotropic Gaussian process with mean
$\mu(\mbox{$\bf s$})=\mu$ and
$K(||\Delta||;\sigma^{2},\phi)=\sigma^{2}\exp(-\phi||\Delta||^{2})$, where
$\Delta=(\delta_{1},\delta_{2})^{\top}$. $\nabla^{k}K(\Delta)$, $k=2,3,4$ is
obtained from (2) and related results.
$\mbox{\boldmath$\gamma$}_{j}(t_{i})=\mbox{\boldmath$\gamma$}_{j}(t_{i};\sigma^{2},\phi)=\left\\{\Phi\left(\sqrt{2\phi}\left(t_{i}+\mbox{$\bf
u$}_{i}^{\top}\Delta_{i-1,j}\right)\right)-\Phi\left(\sqrt{2\phi}\mbox{$\bf
u$}_{i}^{\top}\Delta_{i-1,j}\right)\right\\}(c_{1},c_{2})^{\top}$ where,
$c_{1}=c_{1}(\sigma^{2},\phi,\mbox{$\bf
u$}_{i}^{\perp},\Delta_{i-1,j})=-2\sigma^{2}\sqrt{\pi\phi}{\mbox{$\bf
u$}_{i}^{\perp}}^{\top}\Delta_{i-1,j}e^{-\phi\left({\mbox{$\bf
u$}_{i}^{\perp}}^{\top}\Delta_{i-1,j}\right)^{2}}$,
$c_{2}=c_{2}(\sigma^{2},\phi,\mbox{$\bf
u$}_{i}^{\perp},\Delta_{i-1,j})=c_{1}(1-2\phi{\mbox{$\bf
u$}_{i}^{\perp}}^{\top}\Delta_{i-1,j}\Delta_{i-1,j}^{\top}\mbox{$\bf
u$}_{i}^{\perp})$, and $\Phi(\cdot)$ denotes the standard Gaussian cumulative
distribution function. These are simple computations with quadrature required
only for computing $K_{\mbox{\boldmath$\Gamma$}}(t_{i},t_{i})$.

##  4 Bayesian Hierarchical Model

We operate under a Bayesian hierarchical model, which is specified as

$Y(\mbox{$\bf s$})=\mu(\mbox{$\bf s$},\mbox{\boldmath${\beta}$})+Z(\mbox{$\bf
s$})+\epsilon(\mbox{$\bf s$})\;,$ (13)

where $Z(\mbox{$\bf s$})\sim GP(0,K(\cdot;\sigma^{2},\phi))$ is a Gaussian
process, and $\epsilon(\mbox{$\bf s$})\sim N(0,\tau^{2})$ is a white noise
process, termed as the nugget (see Banerjee et al., 2014, and references
therein). The process parameters are
$\mbox{\boldmath$\theta$}=\\{\mbox{\boldmath${\beta}$},\sigma^{2},\phi,\tau^{2}\\}$.
More generally, we can consider a latent specification for response arising
from exponential families, $\alpha(\mbox{\boldmath$\medmath{\eta}$}(\mbox{$\bf
s$}))=\mbox{$\bf x$}^{\top}(\mbox{$\bf
s$})\mbox{\boldmath${\beta}$}+Z(\mbox{$\bf s$})+\epsilon(\mbox{$\bf s$})$,
$Z(\mbox{$\bf s$})\sim GP(0,K(\cdot;\sigma^{2},\phi))$ and $Y(\mbox{$\bf
s$})\sim\pi\left(\mbox{\boldmath$\medmath{\eta}$}(\mbox{$\bf
s$}),\cdot\right)$, where $\alpha$ is a monotonic link function, $\pi$ is a
member of the exponential family and $\medmath{\eta}$ is the natural
parameter. Predictive inference on differential processes and curvature
wombling proceeds on the latent surface through $P({\cal L}Z\,|\;\mbox{$\bf
Y$})$. The joint posterior for differential processes is obtained through,
$P(\nabla Z^{\top},vech(\nabla^{2}Z)^{\top}\,|\;\mbox{$\bf Y$})=\int P(\nabla
Z^{\top},vech(\nabla^{2}Z)^{\top}\,|\;\mbox{$\bf
Z$},\mbox{\boldmath$\theta$})P(\mbox{$\bf Z$}\,|\;\mbox{$\bf
Y$},\mbox{\boldmath$\theta$})P(\mbox{\boldmath$\theta$}\,|\;\mbox{$\bf
Y$})\,d\mbox{\boldmath$\theta$}\,d\mbox{$\bf Z$}$, while wombling measures
$\mbox{\boldmath$\Gamma$}_{Z}(t^{*})$ for a curve $C_{t^{*}}$ within the
estimated posterior surface for $\bf Z$, are sampled from the posterior,
$P(\mbox{\boldmath$\Gamma$}_{Z}(t^{*})\,|\;\mbox{$\bf Y$})=\int
P(\mbox{\boldmath$\Gamma$}_{Z}(t^{*})\,|\;\mbox{$\bf
Z$},\mbox{\boldmath$\theta$})P(\mbox{$\bf Z$}\,|\;\mbox{$\bf
Y$},\mbox{\boldmath$\theta$})P(\mbox{\boldmath$\theta$}\,|\;\mbox{$\bf
Y$})\,d\mbox{\boldmath$\theta$}\,d\mbox{$\bf Z$}$. Customary prior
specifications for $\theta$ yield

$\begin{split}P(\mbox{\boldmath$\theta$},\mbox{$\bf Z$}\,|\;\mbox{$\bf
Y$})&\propto U(\phi\,|\;a_{\phi},b_{\phi})\times
IG(\sigma^{2}\,|\;a_{\sigma},b_{\sigma})\times
IG(\tau^{2}\,|\;a_{\tau},b_{\tau})\times\mathcal{N}_{L}(\mbox{$\bf
Z$}\,|\;\mathbf{0},\sigma^{2}\mbox{$\bf R$}_{Z})\\\
&\qquad\times\mathcal{N}_{p}(\mbox{\boldmath${\beta}$}\,|\;\mu_{\beta},\Sigma_{\beta})\times\prod\limits_{l=1}^{L}\mathcal{N}_{1}\big{(}Y(\mbox{$\bf
s$}_{l})\,|\;\mbox{$\bf x$}(\mbox{$\bf
s$}_{l})^{\top}\mbox{\boldmath${\beta}$}+Z(\mbox{$\bf
s$}_{l}),\tau^{2}\big{)}\;,\end{split}$ (14)

where $IG$ denotes the inverse-gamma distribution with a shape-rate
parameterization, $U$ is a uniform distribution and $\mbox{$\bf R$}_{Z}$ is
the correlation matrix corresponding to $K(\cdot;\sigma^{2},\phi)$. The
resulting full conditionals are
$\mbox{\boldmath${\beta}$}\,|\;\tau^{2},\mbox{$\bf Z$},\mbox{$\bf
Y$}\sim\mathcal{N}_{p}(M_{\beta}m_{\beta},M_{\beta})$,
$\sigma^{2}\,|\;\phi,\mbox{$\bf Z$}\sim
IG(a_{\sigma}+\frac{L}{2},b_{\sigma}+\frac{1}{2}\mbox{$\bf
Z$}^{\top}\mbox{$\bf R$}^{-1}_{Z}(\cdot;\phi)\mbox{$\bf Z$})$,
$\tau^{2}\,|\;\mbox{\boldmath${\beta}$},\mbox{$\bf Z$},\mbox{$\bf Y$}\sim
IG\left(a_{\tau}+\frac{L}{2},b_{\tau}+\frac{1}{2}||\mbox{$\bf Y$}-\mbox{$\bf
X$}\mbox{\boldmath${\beta}$}-\mbox{$\bf Z$}||_{2}^{2}\right)$, $\mbox{$\bf
Z$}\,|\;\mbox{$\bf
Y$},\mbox{\boldmath$\theta$}\sim\mathcal{N}_{L}(M_{Z}m_{Z},\tau^{2}M_{Z})$,
where $\bf X$ is the $L\times p$ matrix with $\mbox{$\bf x$}(\mbox{$\bf
s$}_{i})^{\top}$ as rows,
$M_{\beta}^{-1}=\Sigma_{\beta}^{-1}+\tau^{-2}\mbox{$\bf X$}^{\top}\mbox{$\bf
X$}$, $m_{\beta}=\Sigma_{\beta}^{-1}\mu_{\beta}+\tau^{-2}\mbox{$\bf
X$}^{\top}(\mbox{$\bf Y$}-\mbox{$\bf Z$})$,
$M_{Z}^{-1}=\tau^{-2}\big{(}\tau^{-2}I_{L}+\sigma^{-2}\mbox{$\bf
R$}_{Z}^{-1}(\cdot;\phi)\big{)}$, and $m_{Z}=\mbox{$\bf Y$}-\mbox{$\bf
X$}\mbox{\boldmath${\beta}$}$. $\phi$ is updated using Metropolis steps with a
normal proposal and an adaptive variance.

Under this setup posterior samples for the differential processes and wombling
measures result from (5) and (6). For each posterior sample of $\\{\mbox{$\bf
Z$},\mbox{\boldmath$\theta$}\\}$, we draw
$\mbox{\boldmath$\Gamma$}_{Z}(t^{*})\,|\;\mbox{$\bf
Z$},\mbox{\boldmath$\theta$}\sim{\cal
N}_{2}\big{(}\mu_{\mbox{\boldmath$\Gamma$}_{Z}}(t^{*})-\mbox{\boldmath$\gamma$}^{\top}_{\mbox{\boldmath$\Gamma$}_{Z}}(t^{*})\Sigma_{\mbox{$\bf
Z$}}^{-1}\mbox{$\bf
Z$},K_{\mbox{\boldmath$\Gamma$}_{Z}}(t^{*},t^{*})-\mbox{\boldmath$\gamma$}^{\top}_{\mbox{\boldmath$\Gamma$}_{Z}}(t^{*})\Sigma_{\mbox{$\bf
Z$}}^{-1}\mbox{\boldmath$\gamma$}_{\mbox{\boldmath$\Gamma$}_{Z}}(t^{*})\big{)}$,
where $\mu_{\mbox{\boldmath$\Gamma$}_{Z}}(t^{*})$,
$\mbox{\boldmath$\gamma$}_{\mbox{\boldmath$\Gamma$}_{Z}}(t^{*})$, and
$K_{\mbox{\boldmath$\Gamma$}_{Z}}(t^{*},t^{*})$ are computed from (10) and
(11). Algorithms 1 and 2 in the Supplement, Section S4, present further
details for posterior sampling. Next, we turn to numerical experiments and
data analyses. Codes required for reproducing and emulating the analyses
presented in the manuscript are produced for the R statistical programming
environment and available for download in the public domain at
https://github.com/arh926/spWombling .

##  5 Simulation Experiments

###  5.1 Data generation

The proposed differential processes are not observed in reality, but are
induced by an observed spatially indexed parent process. To evaluate
statistical learning of the curvature process we perform simulation
experiments within a setup where true values of the differential process and
wombling measures are available. We consider locations $\mbox{$\bf
s$}=(s_{1},s_{2})^{\top}\in\mathbb{R}^{2}$ over the unit square
$[0,1]\times[0,1]\subset\mathbb{R}^{2}$. We generate synthetic data from two
distributions: (a) Pattern 1: $y_{1}(\mbox{$\bf s$})\sim N(10[\sin(3\pi
s_{1})+\cos(3\pi s_{2})],\tau^{2})$; (b) Pattern 2: $y_{2}(\mbox{$\bf s$})\sim
N(10[\sin(3\pi s_{1})\cdot\cos(3\pi s_{2})],\tau^{2})$,

Figure 1: Spatial plots for synthetic patterns, from Pattern 1 (left) and
Pattern 2 (right). Scales are shown in the legend alongside.

where the value of $\tau^{2}=1$. Figure 1 presents spatial plots of the
generated synthetic response from these patterns. The rationale behind
selecting these distributions is: (i) synthetic data is more practical and not
from the model in (13), and (ii) true gradient and curvature can be computed
at every location $\bf s$.

The synthetic patterns chosen feature two different scenarios that may arise.
In the first pattern expressions for differentials along the principal
directions $\mbox{$\bf e$}_{1}=(1,0)^{\top}$ and $\mbox{$\bf
e$}_{2}=(0,1)^{\top}$ are functions of either $s_{1}$ or $s_{2}$,
$\nabla\mu_{1}(\mbox{$\bf s$})=30\pi(\cos(3\pi s_{1}),-\sin(3\pi
s_{2}))^{\top}$, $\nabla^{2}\mu_{1}(\mbox{$\bf s$})=-90\pi^{2}{\rm
diag}\\{\sin(3\pi s_{1}),\cos(3\pi s_{2})\\}$. The curvature along $s_{1}$
does not influence curvature along $s_{2}$,
$\left(\nabla^{2}\mu_{1}(\mbox{$\bf s$})\right)_{12}=0$ for all $\bf s$. While
$\nabla\mu_{2}(\mbox{$\bf s$})=30\pi(\cos(3\pi s_{1})\cos(3\pi
s_{2}),-\sin(3\pi s_{1})\sin(3\pi s_{2}))^{\top}$,
$\nabla^{2}\mu_{2}(\mbox{$\bf s$})=-90\pi^{2}M(\mbox{$\bf s$})$, where
$M(\mbox{$\bf s$})$ is a $2\times 2$ matrix with, $m_{11}=\sin(3\pi
s_{1})\cos(3\pi s_{2})$, $m_{12}=m_{21}=\cos(3\pi s_{1})\sin(3\pi s_{2})$ and
$m_{22}=\sin(3\pi s_{1})\cos(3\pi s_{2})$ with differentials being functions
of both $s_{1}$ and $s_{2}$ and $\left(\nabla^{2}\mu(\mbox{$\bf
s$})\right)_{12}\neq 0$ for some $\bf s$. While setting up the experiments we
vary $L\in\\{100,500,1000\\}$ with 10 replicated instances under each setting.

###  5.2 Bayesian model fitting

We fit the model in (14) with only an intercept allowing the spatial process
to learn the functional patterns in the synthetic response. We use the
following hyper-parameter values in (14): $a_{\phi}=3/\max||\Delta||$,
$b_{\phi}=30$, $a_{\sigma}=2$, $b_{\sigma}=1$, $a_{\tau}=2$, $b_{\tau}=0.1$
$\mu_{\beta}=0$ and $\Sigma_{\beta}=10^{6}I_{p}$. These choices comprise
reasonable weakly informative priors. While a $\rm Uniform(2,3)$ prior on
$\nu$ can be specified (and was implemented as part of this experiment) to
ensure the existence of the curvature process, here our choice of scales in
the data generating patterns ensured that $\nu=5/2$ provided the best model
fit when compared with values of $\nu\in\\{1/2,3/2,5/2\\}$. Hence, we present
the results with $\nu=5/2$.

The parameter estimates for $\theta$ are computed using posterior medians and
their highest posterior density (HPD) intervals (Chen & Shao, 1999; Plummer et
al., 2015). For each replicate, we assess our ability to estimate the local
geometry of the resulting posterior surface. For this we overlay a grid
spanning the unit square. We perform posterior predictive inference for the
differential processes at each grid location following Section 2. Posterior
predictive medians (accompanied by $95$% HPD intervals) summarize inference
for the differential processes over the grid locations (Section 5.4 offers
supplementary analysis).

###  5.3 Bayesian wombling with curvature processes

For wombling with curvature processes, or _curvature wombling_ , we focus on
locating curves that track rapid change within the simulated random surfaces.
For example, consider the surface produced by the first pattern. If a curve is
provided to us, we can evaluate the posterior distribution of the average or
total curvature wombling measures to assess their statistical significance. On
the other hand, without a given curve, we consider three different approaches
for constructing them from a boundary analysis or wombling perspective: (a)
level curves: $C_{y_{0}}=\\{\mbox{$\bf s$}:Y(\mbox{$\bf s$})=y_{0}\\}$:
Bayesian wombling literature finds that curves parallel to contours often form
wombling boundaries (see, e.g., Banerjee & Gelfand, 2006) and level curves on
a surface are parallel to local contours by definition; (b) smooth curves:
produces a smooth curve using Bézier splines (see, e.g., Gallier & Gallier,
2000) from a set of annotated points that are of interest within the surface;
and (c) rectilinear curves: produces a rectilinear curve joining adjacent
annotated points of interest within the surface using straight lines, performs
curvature wombling using a Riemann sum approximation (see (S1) in the
Supplement). Curves of types (b) and (c) allow the investigator to specify a
region of interest that house possible wombling boundaries. For the surface
realization produced by Pattern 1, we consider four different types of curves
on the response surface, (A) a closed curve enclosing a trough corresponding
to a level curve, $C_{y_{0}=-18}$, (B) a closed curve enclosing a peak
corresponding to a level curve, $C_{y_{0}=+18}$, (C) a closed curve that
outlines a contour corresponding to a level curve, $C_{y_{0}=+15}$ and (D) an
open curve along a contour constructed using a Bézier spline. These curves are
marked in Figure 2c.

Figure 2: (left) shows color coded directional gradients for segments,
(center) shows color coded directional curvature for segments in the direction
normal to the curve, (right) shows curves selected for performing curvature
wombling. green indicates a positive significance, cyan indicates negative
significance and white indicates no significance.

Curvature wombling is performed using methods outlined in Section 3. Referring
to the discussion on rectilinear approximation, for each curve, given a
partition, we compute $t_{i}$ and $\mbox{$\bf u$}_{i}$. Combining the segments
produces a vector $\bf t$ and a matrix of directions, $\bf U$ that represents
the curve. Algorithm 2 in the Supplement, Section S4 devises efficient
computation using $\bf t$ and $\bf U$. The total (and average) wombling
measures $\overline{\mbox{\boldmath$\Gamma$}}(C)$ are sampled from their
posteriors using (12). For curves A, B, C and D, we use partitions with
sufficiently small norms ($|P|$) to achieve accuracy ($3.99\times 10^{-3}$,
$3.97\times 10^{-3}$, $4.42\times 10^{-3}$ and $2.66\times 10^{-2}$
respectively). One and two dimensional line integrals (refer to (10) and (11))
are computed via quadrature using grids of size 10 on $[0,t_{i}]$, and size
100 on $[0,t_{i}]\times[0,t_{i}]$ respectively, for $i=1,2,\ldots,n_{P}$. The
median of sampled $\overline{\mbox{\boldmath$\Gamma$}}(\widetilde{C}_{P})$ is
our estimated wombling measure for the curve. Significance at the curve-
segment level is assessed based on the inclusion of 0 within the HPD
intervals. Our design allows us to compute true values of average wombling
measures for each rectilinear segment in the curve. They are computed using,
$\mu_{\small\mbox{\boldmath$\Gamma$}}^{true}(\widetilde{C}_{P})=\left(\sum_{i=1}^{n_{P}}t_{i}\right)^{-1}\left(\sum_{i=1}^{n_{P}}\int_{0}^{t_{i}}{\mbox{$\bf
u$}_{i}^{\perp}}^{\top}\nabla\mu_{1}(\mbox{$\bf
s$}(t))\,dt,\sum_{i=1}^{n_{P}}\int_{0}^{t_{i}}{\mbox{$\bf
u$}_{i}^{\perp}}^{\top}\nabla^{2}\mu_{1}(\mbox{$\bf s$}(t))\mbox{$\bf
u$}_{i}^{\perp}\,dt\right)^{\top}$. We compute HPD intervals for the wombling
measures at the segment level. Coverage probabilities (CPs) are then
constructed by aggregating coverage of true values by HPD intervals over
segments.

Table 1: Results from curvature wombling performed on curves A, B, C and D as
shown in Figure 2. The estimated average directional gradient and curvature
are accompanied by their respective HPD intervals in brackets. HPD intervals
containing 0 are marked in bold.

Curves $(C)$ | Average Gradient ($\Gamma^{(1)}(C)$) | Average Curvature ($\Gamma^{(2)}(C)$)
---|---|---
True | Estimated | True | Estimated
Curve A | -61.54 | -64.97 | 731.94 | 768.03
(-92.37, -38.57) | (599.30, 913.70)
Curve B | 40.85 | 49.19 | -808.04 | -850.84
(20.45, 73.12) | (-1066.98, -630.09)
Curve C | 84.03 | 85.65 | -558.58 | -504.98
(59.81, 109.97) | (-767.55, -241.61)
Curve D | -110.84 | -113.27 | 11.32 | -94.64
(-153.23, -77.01) | (-386.94, 233.78)

Curve A encloses a trough and a local minima for the surface, while B and C
enclose peaks and local maximums (referring to corresponding locations in
Figures S8c and S9c). Along all segments of A we expect negative gradients
owing to the decreasing nature of the response in that region, while for B and
C we expect positive gradients. Each of them would be expected to yield
significant wombling measures for gradients. Referring to the Laplacian
surface (see Supplement, Figures S8e and S9e) A, B, and C are located in
regions manifesting rapid change in the gradient surface, implying they should
yield large positive (curve A) or negative (curves B and C) curvature, forming
curvature wombling boundaries. These are all aligned with our findings
presented in Table 1, which presents measures of quality assessment for
wombling. The magnitude and sign of wombling measures also allow us to
differentiate between the type of curvature for the different wombling
boundaries. For instance, B is located in a region of higher convexity
compared to C, while the nature of convexity for regions enclosed by them are
different compared to A. Plots in Figure 2 (left and center) show line segment
level inference for average wombling measures. Arrows indicate segments which
were not significant with respect to gradient or curvature, while regions of
significance are color coded. D is located in a “relatively flat” region of
the surface (see Figures S8e and S9e) and is expected to have gradients but no
curvature, which aligns with results shown in Table 1. We conclude by noting
that the true values, $\mu_{\small\mbox{\boldmath$\Gamma$}}^{true}(C)$ of the
wombling measures for the curves considered, are all covered by the estimated
HPD intervals for respective curves. Additionally, at the line segment level
we achieved a CP of 1.0 across all curves.

###  5.4 Supplementary analysis

We present additional results in the online supplement. Tables S1 and S2
present parameter estimates, measures of goodness of fit for the fitted
process, and assessment of derivative process characteristics for each pattern
considered. We compute root mean square errors (RMSE) across observed
locations averaged over 10 replicates for each sample size setting for the
fitted process $\widehat{Y}(\mbox{$\bf
s$})=\widehat{\beta_{0}}+\widehat{Z}(\mbox{$\bf s$})$, and $\widehat{\nabla
Y(\mbox{$\bf s$})}$, $\widehat{vech(\nabla^{2}Y(\mbox{$\bf s$})})$. We report
standard deviations across replicates. With increasing number of observed
locations we are able to effectively learn the underlying process and induced
differential processes. Figures S4, S5, S6 and S7 present spatial plots of
posterior medians of gradient and curvature processes, for $L=100$ locations.
These plots demonstrate the effectiveness of our methods in learning about the
differential processes from the underlying patterns. Similarly plots shown in
Figures S8, S9, S10 and S11 demonstrate the same for derived quantities and
operators of ${\cal L}Y(\mbox{$\bf s$})$—principal curvature (eigenvalues),
Gaussian curvature (determinant) (see, e.g., Spivak, 1999; Do Carmo, 2016),
divergence and Laplacian, which pertain to geometric analysis of curvature for
the random surface resulting from the underlying patterns. Statistical
significance is assessed at every grid point by checking the inclusion of 0 in
their HPD intervals. Significantly positive (negative) points are color coded.
We compute average CPs at every grid location to measure the accuracy of our
assessment. These CPs are then averaged over replicates. We observed high CPs
across the grid for parent and differential processes. Figures S12 and S13
compare observed against estimated differential processes coupled with their
HPD regions.

##  6 Applications

Frameworks developed for differential assessment and boundary analysis in
spatially indexed response are applied to multiple data sets with the aim of
locating curvature wombling boundaries that track rapid change in response.
The chosen data arise from varied areas of scientific interest, we briefly
describe the origin and significance of each with respect to our methods
before performing our analysis. Response is modeled using the hierarchical
model in (13). Prior specifications used in (14) are, $\phi\sim{\rm
Unif}\left(3/\max_{\mbox{$\bf s$}\in{\cal S}}{||\Delta||},300\right)$,
$\sigma^{2}\sim IG(2,1)$, $\tau^{2}\sim IG(2,1)$ (mean 1, infinite variance),
$\mbox{\boldmath${\beta}$}\sim N(0,10^{6}I_{p})$, $p$ being the number of
covariates and $\nu=5/2$ for the Matérn kernel ensuring existence of the
differential processes.

Boston Housing: The Boston housing data (see, e.g., Harrison Jr & Rubinfeld,
1978) was collected by the United States Census Service featuring median house
prices for tracts and towns in Boston, Massachusetts area. The purpose was to
study heterogeneity in the market caused by the need for residents to have
clean air. To study such heterogeneity, modern equitable housing policies are
incorporating statistical modeling to quantify such behavior. Often they are a
result of unobserved effects of rapidly shifting socioeconomic conditions
(see, e.g., Hu et al., 2019). Within a spatial map this manifests as
neighboring regions of disparity. Figure 3 shows two such regions: high priced
including Downtown Boston, Cambridge, Newton, Wellesley, Brookline etc. and
low priced including South and East End. For effective policy implementation,
identifying such regions becomes crucial. Spatial variation in the median
house prices is evidenced in Figure 4. Curvature wombling effected on the
house price surface would locate regions that feature such change.

Figure 3: Plots showing (left) probability density of median house prices (in
USD 1000) (right) spatial plot of median owner occupied house prices in
Boston.

The data contains median house price values for 506 census-tracts along with
demographic data. Latitude-longitude centers of the census-tracts are used for
spatial referencing. To allow $Z(\mbox{$\bf s$})$ to capture all the spatial
variation, we include only an intercept in the model. Table 2 shows posterior
estimates and HPD intervals for process parameters. We observe that
$\frac{\sigma^{2}}{\sigma^{2}+\tau^{2}}\approx 78.75\%$—larger portion of
total variance being explained by varying location.

Table 2: Posterior Estimates from the hierarchical linear model in (13) to
Boston housing

Parameters ($\theta$) | Posterior Estimates ($\widehat{\mbox{\boldmath$\theta$}}$) | HPD
---|---|---
$\phi$ | 0.96 | (0.83, 1.11)
$\sigma^{2}$ | 55.18 | (43.91, 68.06)
$\tau^{2}$ | 14.89 | (11.77, 18.71)
$\beta_{0}$ | 25.58 | (24.29, 27.34)

Figure 4: Plots (left to right) showing fitted process, divergence and
Laplacian for the median house price surface.

Modeled spatial variation in the response is shown in Figure 3 (left).
Significance for the estimate, $\widehat{Z}(\mbox{$\bf s$})$, is assessed
using the inclusion of 0 in its posterior HPD. Using posterior samples we
estimate the derivative processes for $Z(\mbox{$\bf s$})$. A grid, (${\cal
G}=\\{\mbox{$\bf s$}_{g}:\mbox{$\bf s$}_{g}\in{\tt convex-hull}({\cal S})\\}$,
containing 1229 equally spaced locations) is overlaid over the region with the
same purpose. To effect posterior surface analysis on the estimated surface we
use posterior predictive distributions of ${\rm div}(Z)$ and $\Delta(Z)$
revealing zones that manifest rapid change in response and gradients
respectively. They are shown in Figures 4 (center and right).

Figure 5: Curvature wombling on the Boston Housing Data.

Next, we focus on performing curvature wombling on the estimated surface.
Strategic posterior surface analysis is used to locate level-sets of interest
within the surface that could possibly contain wombling boundaries. We start
with contours shown in Figure 5 (left column). Boundary 1 (2) bounds a region
where the fitted process has positive (negative) significant estimates.
Evidently, the chosen curves should house significant gradients along most
segments, but significant curvature should only be detected for segments
located at the center (lat-long: $(42.18,42.23)\times(-71.05,-70.05)$) of the
surface in Figures 4 (center and right). Estimated average wombling measures
for these curves are shown in Table 3. Figures 5 (center and right) correspond
to segment level posterior inference for the curves, line segments with
significant directional differentials are indicated in bold. Summarizing, we
observe that the gradient, curvature and posterior surface analysis allow us
to highlight towns (with census-tracts) within Boston that exhibit
heterogeneity in prices. Curvature wombling performed on the surface allows us
to delineate zones that house such heterogeneity. For instance, towns located
within boundaries 3 (South and East End) and 6 (Newton and Brookline) show
significant change in price gradients, compared to towns within boundaries 4
(Lincoln and Weston) and 5 (Wellesley and Dover). These findings can be
verified referring back to price dynamics for real estate in Boston during
1978 (see e.g., Schnare & Struyk, 1976). The same regions are scrutinized for
studying segmentation—towns within curves 1 and 3 are accessible to lower
income groups willing to sacrifice air quality.

Table 3: Curvature wombling measures for boundaries in Boston housing
accompanied by corresponding HPD intervals in brackets below. Estimates
corresponding to HPD intervals containing 0 are marked in bold.

Curve ($C$) | Average Gradient ($\overline{\Gamma}^{(1)}(C)$) | Average Curvature ($\overline{\Gamma}^{(2)}(C)$)
---|---|---
Boundary 1 | -8.91 | 10.14
(-11.31, -6.65) | (2.84, 18.34)
Boundary 2 | 6.18 | -0.09
(4.75, 7.49) | (-3.45, 3.35)
Boundary 3 | -6.47 | 12.69
(-9.74, -3.27) | (2.65, 22.48)
Boundary 4 | 6.92 | 1.26
(4.63, 9.19) | (-5.04, 7.14)
Boundary 5 | 5.47 | 1.36
(2.95, 7.86) | (-4.33, 7.42)
Boundary 6 | 11.82 | -16.27
(7.28, 16.14) | (-26.68,-6.57)

Meuse River Data: The Meuse river data features in Pebesma et al. (2012). It
provides locations of topsoil heavy metal concentrations, along with soil and
landscape variables at the observed locations, collected in a flood plain of
the river Meuse, near the village of Stein, Netherlands. The heavy metal
concentrations recorded include Cadmium (Cd), Copper (Cu), Lead (Pb) and Zinc
(Zn). A distinguishing feature is the naturally occurring boundary—the Meuse.
From a boundary analysis standpoint we are interested in examining
differentials in heavy metal concentrations along the flood plain of the river
to understand the heterogeneous effect of the river on the topsoil.

Figure 6: Plots showing heavy metal concentrations in the topsoil of a flood
plain at 155 locations for (from left to right) Cadmium (Cd), Copper (Cu),
Lead (Pb) and Zinc (Zn) (in mg/kg of soil).

The soils of the floodplain are commonly used for agriculture. Crops grown on
the floodplain of the river banks of the Meuse may be consumed by man and/or
livestock.

Table 4: Posterior estimates of process parameters and covariates for the
Meuse river data accompanied by their corresponding HPD intervals in brackets
below. Effects with HPDs containing 0 are marked in bold.

Parameters ($\theta$) | Cadmium (Cd) | Copper (Cu) | Lead (Pb) | Zinc (Zn)
---|---|---|---|---
$\phi$ | 0.0379 | 0.1138 | 0.0399 | 0.0472
(0.0207, 0.0618) | (0.0871, 0.1471) | (0.0131, 0.1900) | (0.0230, 0.0744)
$\sigma^{2}$ | 2.9566 | 3.2044 | 0.9303 | 38.3538
(1.2803, 5.2227) | (2.3955, 4.0892) | (0.2763, 1.7641) | (16.7815, 65.1450)
$\tau^{2}$ | 1.7771 | 0.0067 | 0.8555 | 23.2226
(0.9107, 2.6328) | (0.0012, 0.0244) | (0.0010, 1.2280) | (9.3743, 35.6867)
Intercept | 9.4973 | 4.8503 | 6.1120 | 37.1315
(5.9750, 13.3704) | (3.1392, 6.8308) | (3.4615, 8.1910) | (25.0903, 53.0870)
elev | -0.7672 | -0.4065 | -0.5413 | -2.8781
(-1.2531, -0.3574) | (-0.7418, -0.1656) | (-0.7853, -0.1442) | (-4.7805, -1.2834)
om | 0.4011 | 0.4293 | 0.3434 | 0.8606
(0.2616, 0.5233) | (0.3276, 0.4728) | (0.2490, 0.4253) | (0.3681, 1.3166)
dist | -0.0033 | -0.0025 | -0.0011 | -0.0081
(-0.0061, 0.0000) | (-0.0043, -0.0014) | (-0.0029, 0.0006) | (-0.0197, 0.0038)
ffreq (=2) | -1.4176 | -2.4727 | -0.8483 | -4.3182
(-2.3202, -0.3432) | (-3.1794, -1.6716) | (-1.6109, -0.2598) | (-7.9184, -0.6220)
ffreq (=3) | -0.7322 | -1.4298 | -0.1865 | -3.3159
(-2.0520, 0.6248) | (-2.4443, -0.5157) | (-1.2972, 0.6861) | (-7.9307, 1.9128)
soil (=2) | -0.3337 | 0.2236 | 0.5988 | -2.2213
(-1.4661, 0.7491) | (-0.7248, 0.9799) | (-0.0345, 1.2956) | (-6.1446, 2.0831)
soil (=3) | -0.3884 | 0.6344 | 0.3707 | -2.9922
(-2.0891, 1.2628) | (-0.2309, 1.8474) | (-0.7108, 1.4029) | (-9.0918, 3.6289)
lime (=1) | 0.5752 | 1.3223 | 0.7759 | -0.4759
(-0.3509, 1.4341) | (0.7152, 1.9427) | (0.1173, 1.4645) | (-3.9057, 2.6510)

The spatial variation in heavy metal concentration can be seen in Figure 6.
The path of the Meuse river is shown in each of the spatial plots. Evidently,
the heavy metal concentrations decreases with increasing distance from the
river. We model the concentrations as independent Gaussian processes.
Covariates used are relative elevation above local river bed (elev, measured
in meters), organic matter (om measured in kg/(100kg) of soil), distance to
Meuse (dist), frequency of flooding, soil type (soil), and lime content in
soil ($p=9$). Table 4 shows the posterior estimates of process parameters and
model coefficients, ${\beta}$ for each of the heavy metals in question. We
observe that $\sigma^{2}/(\sigma^{2}+\tau^{2})\approx$ 62.45%, 99.79%, 52.09%,
62.29% for Cd, Cu, Pb and Zn respectively, indicating larger portions of total
variation being explained by spatial heterogeneity, except for Pb. Variation
in Cd and Zn concentration is significantly affected by elevation, organic
matter and flooding frequency, while variations in Cu and Pb concentration is
significantly affected by elevation, organic matter and flooding frequency and
lime content. The estimated residual surface is shown in Figure 7 (left) for
Cd concentrations. We observe significant positive gradients with varying
curvature depending on segments of the river bed for all heavy metals. We
perform curvature wombling on the Meuse using the residual surface, $\bf Z$.
The results of curvature wombling for cadmium are shown in Figure 7. Results
and plots for other metals can be found in the Supplement, Section S7, Figure
S14. The accompanying wombling measures are shown in Table 5. We observe
sufficient heterogeneity in the signs of the wombling measures, yielding
contiguous positive (negative) segments. For example, in Cd concentration,
boundaries located for average gradients in the northern and southern region
are positive, as opposed to boundaries located in the north western region.
Therefore, while displaying the wombling measures, in Table 5, we separate
them by their sign.

Figure 7: Plots showing results for curvature wombling on the Meuse river for
Cadmium (Cd) concentration. Plots showing (left) the resulting fitted process
(center) the contiguous segments that display significant gradients (right)
the contiguous segments with significant curvature.

We conclude that effects of river Meuse on regions of the flood plain exhibit
significant heterogeneity when considered across heavy metals. Compared to
other metals, Pb concentrations are limited to northern regions of the flood
plain. Concentrations of Cd and Zn concentrations along the river are similar.
Compared to the northern region, in the northwestern region Zn concentrations
decrease significantly as we move inland. Studies corroborating such evidence
can be found in Leenaers et al. (1988) and Albering et al. (1999).

##  7 Discussion and Future Work

We developed a fully model-based Bayesian inferential framework for
differential process assessment and curvature-based boundary analysis for
spatial processes. Introducing the directional curvature process and its
associated inferential framework supplements the directional gradients with
inference for their rates of change, while its induction into the folds of
Bayesian curvilinear wombling allows for further characterization of
difference boundaries. Adopting a Bayesian hierarchical model allows for
Gaussian calibration when characterizing points, regions and boundaries within
a surface. This framework is widely applicable; our applications arise from
selected disciplines indicating the utilities of mapping curvature process
boundaries to understand spatial data generating patterns. Substantive case
studies will be reported separately.

Table 5: Curvature wombling measures for the Meuse, separated by zones of
positive and negative signs, they are accompanied by their corresponding HPD
intervals in brackets below.

Wombling Measures | Cd | Cu | Pb | Zn
---|---|---|---|---
$\overline{\Gamma}^{(1)}(>0)$ | 0.0510 | 0.1273 | 0.0375 | 0.1984
(0.0298, 0.07401) | (0.0913, 0.1729) | (0.0019, 0.1561) | (0.0876, 0.3162)
$\overline{\Gamma}^{(1)}(<0)$ | -0.0400 | -0.2561 | – | -0.1890
(-0.0635, -0.0170) | (-0.3187, -0.1997 ) | – | (-0.2967, -0.0669)
$\overline{\Gamma}^{(2)}(>0)$ | 0.0074 | – | – | 0.0422
(0.0019, 0.0158) | – | – | (0.0111, 0.0879)
$\overline{\Gamma}^{(2)}(<0)$ | -0.0078 | -0.1247 | -0.0039 | -0.0473
(-0.0223, -0.0024) | (-0.1979, -0.0860) | (-0.1076, -0.0006) | (-0.1114, -0.0095)

Several avenues hold scope for future developments. A more generalized
theoretical framework can be developed for studying joint behavior of the
principal curvature (direction of maximum (or minimum) curvature) and the
aspect (direction of maximum gradient) (see, e.g., Wang et al., 2018)
leveraging dependent circular uniform distributions (see, e.g., Kent et al.,
2008). We offer some brief remarks. To obtain the direction of maximum
curvature for a spatial surface, we solve $\displaystyle{\max_{\mbox{$\bf
u$}\in\mathbb{R}^{2}}\left|\mbox{$\bf u$}^{\top}\nabla^{2}Y(\mbox{$\bf
s$})\mbox{$\bf u$}\right|}$, such that $||\mbox{$\bf u$}||=1$, at an arbitrary
point $\bf s$. Using Lagrange multipliers and denoting $\kappa(\mbox{$\bf
u$})=|\mbox{$\bf u$}^{\top}\nabla^{2}Y(\mbox{$\bf s$})\mbox{$\bf u$}|$, define
${\cal O}(\mbox{$\bf u$})=\kappa(\mbox{$\bf u$})-\lambda(||\mbox{$\bf
u$}||^{2}-1)$ hence, $\partial{\cal O}(\mbox{$\bf u$})/\partial
u_{i}=\kappa(\mbox{$\bf u$})^{-1}(\nabla^{2}_{ii}Y(\mbox{$\bf
s$})u_{i}+\nabla^{2}_{ij}Y(\mbox{$\bf s$})u_{j})-\lambda u_{i}=0$, $i,j=1,2$.
With $u_{2}/u_{1}=\tan\theta_{pc}$, eliminating $\lambda$ we get
$\displaystyle{\tan\theta_{pc}=\frac{\nabla^{2}_{22}Y(\mbox{$\bf
s$})\tan\theta_{pc}+\nabla^{2}_{12}Y(\mbox{$\bf
s$})}{\nabla^{2}_{11}Y(\mbox{$\bf s$})+\nabla^{2}_{12}Y(\mbox{$\bf
s$})\tan\theta_{pc}}}$. Defining $h_{1}=h_{1}(\mbox{$\bf
s$})=(\nabla^{2}_{11}Y(\mbox{$\bf s$})-\nabla^{2}_{22}Y(\mbox{$\bf
s$}))/\nabla^{2}_{12}Y(\mbox{$\bf s$})$ given $\nabla^{2}_{12}Y(\mbox{$\bf
s$})\neq 0$ and solving
$\theta_{pc}=\tan^{-1}\frac{1}{2}\left[-h_{1}\pm\sqrt{h_{1}^{2}+4}\right]$. If
$\nabla^{2}_{12}Y(\mbox{$\bf s$})=0$ then, $\nabla^{2}Y(\mbox{$\bf s$})$ is
diagonal and $\theta_{pc}$ corresponds to the direction of
$\max\\{\nabla^{2}_{11}Y(\mbox{$\bf s$}),\nabla^{2}_{22}Y(\mbox{$\bf s$})\\}$.
We propose that $\Theta=(\theta_{asp},\theta_{pc})^{\top}$ follows a dependent
circular uniform distribution over $[0,2\pi]\times[0,2\pi]$. Further
developments with circular regression methods can proceed to examine the
effect of covariates on $\Theta$. Multivariate extensions would involve
formulating these differential processes on arbitrary manifolds. This requires
simulating a Gaussian process on manifolds and inspecting the covariant
derivative. Bayesian curvilinear wombling could then be implemented on curves
of interest to the investigator. This would not only involve an inferential
framework for normal curvature, but also geodesic curvature for such curves.
Spatiotemporal curvature processes can build upon Quick et al. (2015) to study
evolutionary behavior of the curvature processes with respect to variations in
the response across time. Finally, we remark that while there have been
substantial recent developments in scalable spatial processes for massive data
sets—a comprehensive review is beyond the scope of the current article (see,
e.g., Heaton et al., 2019)—not all scalable processes admit the correct degree
of smoothness for curvature processes to exist. Constructing scalable
processes for curvilinear wombling, and subsequent inference, remains a
problem of interest in the wombling community.

## Supplementary Materials

The following supplement includes additional theoretical derivations,
computing details, additional simulation experiments and wombling for
Northeastern US temperatures.

Supplementary Materials for
“Bayesian Modeling with Spatial Curvature Processes”

##  8 Review of Directional Gradients and Wombling

###  8.1 Directional Gradients

For the scalar $h$ and unit vector $\bf u$ we define $Y^{(1)}_{\mbox{$\bf
u$},h}=\left(Y(\mbox{$\bf s$}+h\mbox{$\bf u$})-Y(\mbox{$\bf s$})\right)/h$ to
be the first order finite difference processes at location $\bf s$ in the
directions of $\bf u$. Being a linear function of stationary processes this is
well-defined. Passing to limits, we define $D^{(1)}_{\mbox{$\bf
u$}}Y(\mbox{$\bf s$})=\lim_{h\to 0}Y^{(1)}_{\mbox{$\bf u$},h}(\mbox{$\bf
s$})$. Provided the limit exist, $D^{(1)}_{\mbox{$\bf u$}}Y(\mbox{$\bf s$})$
is defined as the directional gradient process. If $Y(\mbox{$\bf s$})$ is a
mean square differentiable process in $\mathbb{R}^{d}$ for every $\mbox{$\bf
s$}_{0}\in\mbox{$\mathbb{R}^{d}$}$ then $D^{(1)}_{\mbox{$\bf u$}}Y(\mbox{$\bf
s$})=\mbox{$\bf u$}^{\top}\nabla Y(\mbox{$\bf s$})$. Then, if $\mbox{$\bf
u$}=\sum_{i=1}^{d}u_{i}\mbox{$\bf e$}_{i}$, we can compute
$D^{(1)}_{\mbox{$\bf u$}}Y(\mbox{$\bf
s$})=\sum_{i=1}^{d}u_{i}D^{(1)}_{\mbox{$\bf e$}_{i}}Y(\mbox{$\bf s$})$. The
directional gradient process is linear in $\bf u$, hence $D^{(1)}_{-\mbox{$\bf
u$}}Y(\mbox{$\bf s$})=-D^{(1)}_{\mbox{$\bf u$}}Y(\mbox{$\bf s$})$ and for any
vector $\mbox{$\bf w$}=||\mbox{$\bf w$}||\mbox{$\bf u$}$, $D^{(1)}_{\mbox{$\bf
w$}}Y(\mbox{$\bf s$})=||\mbox{$\bf w$}||D^{(1)}_{\mbox{$\bf u$}}Y(\mbox{$\bf
s$})$. The directional gradient at $\mbox{$\bf s$}_{0}$ in the direction $\bf
w$ is the slope at $\mbox{$\bf s$}_{0}$ of the curve traced out by slicing
$Y(\mbox{$\bf s$})$ in the direction $\bf w$ (see e.g., Banerjee & Gelfand,
2006, Section 2, for more details).

###  8.2 Wombling Measures

Wombling measures constructed from (7) for total and average gradient are
associated with curves to characterize the magnitude of change. To each point
$\mbox{$\bf s$}\in C$, a directional gradient is associated, $g\left({\cal
L}Y(\mbox{$\bf s$})\right)=D^{(1)}_{\mbox{$\bf n$}}Y(\mbox{$\bf
s$})=\mbox{$\bf n$}(\mbox{$\bf s$})^{\top}\nabla Y(\mbox{$\bf s$})$ (also a
linear function of $\displaystyle{{\cal L}_{\mbox{$\bf n$}}Y(\mbox{$\bf
s$})}$), along the direction of a unit normal $\mbox{$\bf u$}=\mbox{$\bf
n$}(\mbox{$\bf s$})$ to the curve. For a curve tracking rapid change in the
surface, choice of the normal direction to a curve is motivated by sharp
directional gradients orthogonal to the curve; $\ell$ is chosen to be the arc-
length measure. The rationale behind this choice is to measure change in
response with respect to distance traversed on the curve. With reference to
(7) the total and average gradients are,
${{\Gamma}^{(1)}(C)=\int_{t_{0}}^{t_{1}}\nabla Y(\mbox{$\bf
s$}(t))^{\top}\mbox{$\bf n$}(\mbox{$\bf s$}(t))||\mbox{$\bf
s$}^{\prime}(t)||dt}$ and
$\overline{\Gamma}^{(1)}({C})={\Gamma}^{(1)}(C)/\ell({C})$ respectively (see
e.g., Banerjee & Gelfand, 2006, Section 3, for more details).

##  9 Interpretation of Spatial Curvature

Pursuits in geospatial analysis generally encounter surfaces which have
canonical coordinate systems (e.g. latitude-longitude, easting-northing etc.).
This facilitates a parameterization for the surface that leverages the
coordinate system, commonly known as the Monge parameterization (also called a
Monge patch, named after Gaspard Monge, see e.g., O’Neill, 2006; Pressley,
2010)—a surface, $S$, embedded in $\mathbb{R}^{3}$, is parameterized by giving
its height $Y$ over some plane as a function of the orthonormal co-ordinates
$s_{1}$ and $s_{2}$ in the plane, $S=\\{{\cal
S}\subset\mbox{$\mathbb{R}^{2}$}\mapsto\mbox{$\mathbb{R}^{3}$}:\mbox{$\bf
s$}=(s_{1},s_{2})\mapsto Y(s_{1},s_{2})=Y(\mbox{$\bf s$})\\}$. A point is
then, $(\mbox{$\bf s$},Y(\mbox{$\bf s$}))=(s_{1},s_{2},Y(s_{1},s_{2}))$. The
two tangent vectors at $\bf s$ are, $\mbox{$\bf E$}_{1}(\mbox{$\bf
s$})=(1,0,\nabla_{1}Y(\mbox{$\bf s$}))^{\top}$ and $\mbox{$\bf
E$}_{2}(\mbox{$\bf s$})=(0,1,\nabla_{2}Y(\mbox{$\bf s$}))^{\top}$, where
$\nabla_{i}Y(\mbox{$\bf s$})=\displaystyle{\frac{\partial}{\partial\mbox{$\bf
s$}_{i}}Y(\mbox{$\bf s$})},\;i=1,2$. Let $\nabla Y(\mbox{$\bf
s$})=(\nabla_{1}Y(\mbox{$\bf s$}),\nabla_{2}Y(\mbox{$\bf s$}))^{\top}$ denote
the gradient vector, consider a unit direction vector, $\mbox{$\bf
u$}=(u_{1},u_{2})^{\top}\in S\subset\mbox{$\mathbb{R}^{2}$}$, then
$u_{1}\mbox{$\bf E$}_{1}(\mbox{$\bf s$})+u_{2}\mbox{$\bf E$}_{2}(\mbox{$\bf
s$})=(\mbox{$\bf u$}^{\top},\mbox{$\bf u$}^{\top}\nabla Y(\mbox{$\bf
s$}))^{\top}\in T_{S}(\mbox{$\bf s$})$ corresponds to the directional
derivative of $Y$ along the direction $\bf u$, where $T_{S}(\mbox{$\bf s$})$
is the local tangent plane at $\bf s$, that is generated by $\\{\mbox{$\bf
E$}_{1}(\mbox{$\bf s$}),\mbox{$\bf E$}_{2}(\mbox{$\bf s$})\\}$. The outward
pointing normal to the surface $S$, denoted by $\mbox{$\bf N$}(\mbox{$\bf
s$})=\mbox{$\bf E$}_{1}(\mbox{$\bf s$})\times\mbox{$\bf E$}_{2}(\mbox{$\bf
s$})=(-\nabla_{1}Y(\mbox{$\bf s$}),-\nabla_{2}Y(\mbox{$\bf s$}),1)^{\top}$,
where $\times$ denotes the usual cross-product of vectors. Evidently,
$\mbox{$\bf N$}(\mbox{$\bf s$})$ is orthogonal to the local tangent plane at
that point, $T_{S}(\mbox{$\bf s$})$. Quantifying the local geometry of a
surface, we are interested in how $\mbox{$\bf N$}(\mbox{$\bf s$})$ changes
(“tips”) as we move in the direction $\bf u$ from the point $\bf s$ on the
surface—derivatives for $\mbox{$\bf N$}(\mbox{$\bf s$})$ at the point $\bf s$,
which lie in $T_{S}(\mbox{$\bf s$})$. This is quantified by the normal
curvature of $S$ along a direction $\bf u$. Before defining normal curvature
for a surface, we digress briefly to investigate effects of surface curvature
on curves—for a curve parameterized by $t$, $C=\\{\mbox{$\bf
s$}(t)=(s_{1}(t),s_{2}(t)):t\in[a,b]\\}$, passing through $\bf s$, if the
curvature of $C$ is, $\kappa$, the tangent to $C$ at $\bf s$, $\mbox{$\bf
t$}(\mbox{$\bf s$})$, and the principal unit normal, i.e. the normal to $C$ on
$S$, $\mbox{$\bf n$}(\mbox{$\bf s$})=\mbox{$\bf t$}^{\prime}(\mbox{$\bf
s$})/\kappa$, then we have, $\mbox{$\bf n$}(\mbox{$\bf s$})\cdot\mbox{$\bf
N$}(\mbox{$\bf s$})=\cos(\theta)$, $\mbox{$\bf t$}^{\prime}(\mbox{$\bf
s$})=\kappa\mbox{$\bf n$}(\mbox{$\bf s$})$, which implies
$\kappa\cos(\theta)=\mbox{$\bf t$}^{\prime}(\mbox{$\bf s$})\cdot\mbox{$\bf
N$}(\mbox{$\bf s$})$. We observe that,

$s_{1}$
---
$s_{2}$
---
$Y(s_{1},s_{2})$
---
$S=\left({\bf s},Y({\bf s})\right)$
---
${\bf N}({\bf s})$
---
${\bf T}_{\bf S}({\bf s})$
---
$\theta$
---
${\bf n}({\bf s})$
---
${\bf C}$
---
Figure 8: A Monge patch, $S=(s_{1},s_{2},Y(s_{1},s_{2}))^{\top}=(\mbox{$\bf
s$}^{\top},Y(\mbox{$\bf s$}))^{\top}$, showing a point $\bf s$, a curve $C$
passing through $\bf s$, normal to the surface, $\mbox{$\bf N$}(\mbox{$\bf
s$})$, normal to the curve $\mbox{$\bf n$}(s)$, $\theta$ is the angle between
them, and the local tangent plane to the surface, $T_{S}(\mbox{$\bf s$})$. The
thin perpendicular pink arrows are tangent vectors, $\mbox{$\bf
E$}_{1}(\mbox{$\bf s$})$ and $\mbox{$\bf E$}_{2}(\mbox{$\bf s$})$. The thin
outward pointing black arrows around $\mbox{$\bf N$}(\mbox{$\bf s$})$
demonstrate change in $\bf N$ as we move along the direction (dotted line) on
the surface.

$\mbox{$\bf t$}^{\prime}(\mbox{$\bf s$}(t))=\mfrac{\partial^{2}}{\partial
t^{2}}\left(\mbox{$\bf s$}(t),Y(\mbox{$\bf
s$}(t))\right)^{\top}=\frac{\partial}{\partial t}\mbox{$\bf E$}_{i}(\mbox{$\bf
s$}(t))s_{i}^{\prime}(t)=\mbox{$\bf E$}_{ij}(\mbox{$\bf
s$}(t))s_{i}^{\prime}(t)s_{j}^{\prime}(t)+\mbox{$\bf E$}_{i}(\mbox{$\bf
s$}(t))s_{i}^{\prime\prime}(t)\;,$

since $\mbox{$\bf E$}_{i}(\mbox{$\bf s$})\cdot\mbox{$\bf N$}(\mbox{$\bf
s$})=0$, $\kappa\cos(\theta)=\mbox{$\bf t$}^{\prime}(\mbox{$\bf
s$})\cdot\mbox{$\bf N$}(\mbox{$\bf s$})=(\mbox{$\bf E$}_{ij}(\mbox{$\bf
s$})\cdot\mbox{$\bf N$}(\mbox{$\bf s$}))\partial s_{i}\partial s_{j}$,
$i,j=1,2$, where the expression in the parenthesis is a property of the
surface, independent of curve $C$, and is defined as the second fundamental
form,

$\Pi(\mbox{$\bf s$})=\begin{pmatrix}\mbox{$\bf E$}_{11}(\mbox{$\bf
s$})\cdot\mbox{$\bf N$}(\mbox{$\bf s$})&\mbox{$\bf E$}_{12}(\mbox{$\bf
s$})\cdot\mbox{$\bf N$}(\mbox{$\bf s$})\\\ \mbox{$\bf E$}_{21}(\mbox{$\bf
s$})\cdot\mbox{$\bf N$}(\mbox{$\bf s$})&\mbox{$\bf E$}_{22}(\mbox{$\bf
s$})\cdot\mbox{$\bf N$}(\mbox{$\bf
s$})\end{pmatrix}=\begin{pmatrix}\nabla^{2}_{11}Y(\mbox{$\bf
s$})&\nabla^{2}_{12}Y(\mbox{$\bf s$})\\\ \nabla^{2}_{21}Y(\mbox{$\bf
s$})&\nabla^{2}_{22}Y(\mbox{$\bf s$})\end{pmatrix}=\nabla^{2}Y(\mbox{$\bf
s$})\;,$

where $\mbox{$\bf E$}_{ij}$ and $\nabla^{2}_{ij}Y$ are the partial
differentiation of $\bf E$ and $Y$ with respect to $s_{i},s_{j}$ respectively,
$\cdot$ is the usual dot product for vectors and
$\nabla^{2}_{12}Y=\nabla^{2}_{21}Y$. The second to last equality is obtained
under the Monge parameterization. The second fundamental form is invariant
with respect to transformations of the local co-ordinate which preserves the
sense of $\bf N$, i.e. the transformation does not change an outward (inward)
pointing normal to an inward (outward) pointing normal for $S$. Such surfaces
are termed as orientable surfaces. The Möbius transformation is an example of
non-orientable surfaces. The individual terms of $\Pi$, quantify the local
geometry of a surface (or curvature) along orthonormal coordinates.

Curvature of the $C$ can be attributed to (a) the curvature of the curve
itself, and (b) the curvature of the surface on which $C$ lies. $\kappa$ is
the curvature of $C$, termed as geodesic curvature. The curvature of the
surface is termed as normal curvature, computed along a direction $\mbox{$\bf
u$}=\mbox{$\bf u$}(\mbox{$\bf s$})$ is denoted by $\kappa_{n}(\mbox{$\bf
u$})$. The normal curvature, which is an intrinsic property of the surface
independent of $C$, is of primary interest to us, $\kappa_{n}(\mbox{$\bf
u$})=\mbox{$\bf u$}^{\top}\Pi(\mbox{$\bf s$})\mbox{$\bf u$}/\mbox{$\bf
u$}^{\top}\mbox{$\bf u$}=\mbox{$\bf u$}^{\top}\Pi\mbox{$\bf u$}$, under
$\mbox{$\bf u$}^{\top}\mbox{$\bf u$}=1$. $\kappa_{n}(\mbox{$\bf u$})$ is also
the directional curvature of the $Y$ along $\bf u$. For our exploits,
$\mbox{$\bf u$}=\mbox{$\bf n$}$, the normal direction to $C$. The sign of
$\kappa_{n}(\mbox{$\bf u$})$, or equivalently eigen-values of $\Pi$ inform
about the nature of curvature at $\bf s$—for example, if $\kappa_{1}$ and
$\kappa_{2}$ denote eigenvalues of $\Pi$, with
$K=\det\Pi=\kappa_{1}\kappa_{2}$, if $K>0$, it implies that the surface is
bending away from $T_{S}(\mbox{$\bf s$})$; depending on whether
$\kappa_{1},\kappa_{2}<0$ (or $>0$), $\bf s$ can be locally classified as a
concave (convex) ellipsoid (for more details see Stevens, 1981).

For a purely differential geometric treatment of this discussion see—Gauss
(1902); Spivak (1999); Do Carmo (2016); Kreyszig (2019). Figure 8 illustrates
this discussion.

##  10 Examples for selected Covariance Functions

The detailed calculations for closed form expression of selected covariance
functions are presented. We start with the power exponential family of
isotropic covariance functions,
$\widetilde{K}(||\Delta||)=\alpha\exp(-\phi||\Delta||^{\nu})$, $0<\nu\leq 2$.
It is clear that $\nabla^{4}\widetilde{K}(||\Delta||)$ and exists only for
$\nu=2$, we have the following form,
$\widetilde{K}(||\Delta||)=\sigma^{2}\exp(-\phi||\Delta||^{2})$. For this
choice we have,

$\displaystyle\left(\nabla\widetilde{K}(\Delta)\right)_{i}=-2\sigma^{2}\phi\exp(-\phi||\Delta||^{2})\delta_{i},$
$\displaystyle\left(\nabla^{2}\widetilde{K}(\Delta)\right)_{ii}=-2\sigma^{2}\phi\exp(-\phi||\Delta||^{2})(1-2\phi\delta_{i}^{2}),$
$\displaystyle\left(\nabla^{2}\widetilde{K}(\Delta)\right)_{ij}=4\sigma^{2}\phi^{2}\exp(-\phi||\Delta||^{2})\delta_{i}\delta_{j},$
$\displaystyle\left(\nabla^{3}\widetilde{K}(\Delta)\right)_{iii}=4\sigma^{2}\phi^{2}\exp(-\phi||\Delta||^{2})(3-2\phi\delta_{i}^{2})\delta_{i},$
$\displaystyle\left(\nabla^{3}\widetilde{K}(\Delta)\right)_{iij}=4\sigma^{2}\phi^{2}\exp(-\phi||\Delta||^{2})(1-2\phi\delta_{i}^{2})\delta_{j},$
$\displaystyle\left(\nabla^{3}\widetilde{K}(\Delta)\right)_{ijk}=-8\sigma^{2}\phi^{3}\exp(-\phi||\Delta||^{2})\delta_{i}\delta_{j}\delta_{k},$
$\displaystyle\left(\nabla^{4}\widetilde{K}(\Delta)\right)_{iiii}=4\sigma^{2}\phi^{2}\exp(-\phi||\Delta||^{2})(3-12\phi\delta_{i}^{2}+4\phi^{2}\delta_{i}^{4}),$
$\displaystyle\left(\nabla^{4}\widetilde{K}(\Delta)\right)_{iiij}=-8\sigma^{2}\phi^{3}\exp(-\phi||\Delta||^{2})(3-2\phi\delta_{i}^{2})\delta_{i}\delta_{j},$
$\displaystyle\left(\nabla^{4}\widetilde{K}(\Delta)\right)_{iijj}=4\sigma^{2}\phi^{2}\exp(-\phi||\Delta||^{2})(1-2\phi\delta_{i}^{2})(1-2\phi\delta_{j}^{2}),$
$\displaystyle\left(\nabla^{4}\widetilde{K}(\Delta)\right)_{ijkl}=16\sigma^{2}\phi^{4}\exp(-\phi||\Delta||^{2})\delta_{i}\delta_{j}\delta_{k}\delta_{l},$

where $i,j,k,l=1,2,\ldots,d$. The squared exponential or the Gaussian
covariance kernel is the only member of its class that admits such
derivatives, although they have been critiqued to produce realizations that
are too smooth to be of practical use in modeling (see Stein (1999)).

Turning to the Matérn class we see that with,
$\widetilde{K}||\Delta||=\alpha(\phi||\Delta||)^{\nu}K_{\nu}(\phi||\Delta||)$,
where $\nu$ is a parameter controlling the smoothness of realizations, that is
mean square differentiability and $K_{\nu}$ is the modified Bessel function of
order $\nu$. At $\nu=3/2$ and $\nu=5/2$, $\widetilde{K}(||\Delta||)$ takes the
forms,

$\displaystyle\widetilde{K}(||\Delta||)=\begin{cases}\sigma^{2}(1+\sqrt{3}\phi||\Delta||)e^{-\sqrt{3}\phi||\Delta||},&\nu=3/2\\\
\sigma^{2}\left(1+\sqrt{5}\phi||\Delta||+\mfrac{5}{3}\phi^{2}||\Delta||^{2}\right)e^{-\sqrt{5}\phi||\Delta||},&\nu=5/2\end{cases},$

where $\sigma^{2}$ is the overall process variance. Matérn with $\nu=3/2$ is
once mean square differentiable, where as Matérn with $\nu=5/2$ is twice mean
square differentiable at 0. As $\nu\to\infty$, Matérn covariances tend to the
Gaussian covariance. Unlike the Gaussian covariance, they do not yield overly
smoothed process realizations. For $\nu=3/2$ we have,

$\displaystyle\left(\nabla\widetilde{K}(\Delta)\right)_{i}=-3\sigma^{2}\phi^{2}e^{-\sqrt{3}\phi||\Delta||}\delta_{i},$
$\displaystyle\left(\nabla^{2}\widetilde{K}(\Delta)\right)_{ii}=-3\sigma^{2}\phi^{2}e^{-\sqrt{3}\phi||\Delta||}\left(1-\sqrt{3}\phi\mfrac{\delta_{i}^{2}}{||\Delta||}\right),$
$\displaystyle\left(\nabla^{2}\widetilde{K}(\Delta)\right)_{ij}=3\sqrt{3}\sigma^{2}\phi^{3}e^{-\sqrt{3}\phi||\Delta||}\mfrac{\delta_{i}\delta_{j}}{||\Delta||}.$

where $i,j=1,2,\ldots,d$. Since the process is just once mean square
differentiable higher order derivatives do not exist. However, for $\nu=5/2$
we have,

$\displaystyle\left(\nabla\widetilde{K}(\Delta)\right)_{i}=-\mfrac{5}{3}\sigma^{2}\phi^{2}e^{-\sqrt{5}\phi||\Delta||}\left(1+\sqrt{5}\phi||\Delta||\right)\delta_{i},$
$\displaystyle\left(\nabla^{2}\widetilde{K}(\Delta)\right)_{ii}=-\mfrac{5}{3}\sigma^{2}\phi^{2}e^{-\sqrt{5}\phi||\Delta||}\left(1+\sqrt{5}\phi||\Delta||-5\phi^{2}\delta_{i}^{2}\right),$
$\displaystyle\left(\nabla^{2}\widetilde{K}(\Delta)\right)_{ij}=\mfrac{25}{3}\sigma^{2}\phi^{4}e^{-\sqrt{5}\phi||\Delta||}\delta_{i}\delta_{j},$
$\displaystyle\left(\nabla^{3}\widetilde{K}(\Delta)\right)_{iii}=\mfrac{25}{3}\sigma^{2}\phi^{4}e^{-\sqrt{5}\phi||\Delta||}\left(3-\sqrt{5}\phi\mfrac{\delta_{i}^{2}}{||\Delta||}\right)\delta_{i},$
$\displaystyle\left(\nabla^{3}\widetilde{K}(\Delta)\right)_{iij}=\mfrac{25}{3}\sigma^{2}\phi^{4}e^{-\sqrt{5}\phi||\Delta||}\left(1-\sqrt{5}\phi\mfrac{\delta_{i}^{2}}{||\Delta||}\right)\delta_{j},$
$\displaystyle\left(\nabla^{3}\widetilde{K}(\Delta)\right)_{ijk}=-\mfrac{25\sqrt{5}}{3}\sigma^{2}\phi^{5}e^{-\sqrt{5}\phi||\Delta||}\mfrac{\delta_{i}\delta_{j}\delta_{k}}{||\Delta||},$
$\displaystyle\left(\nabla^{4}\widetilde{K}(\Delta)\right)_{iiii}=\mfrac{25}{3}\sigma^{2}\phi^{4}e^{-\sqrt{5}\phi||\Delta||}\left[3-6\sqrt{5}\phi\mfrac{\delta_{i}^{2}}{||\Delta||}+\sqrt{5}\phi\left(\sqrt{5}\phi+\mfrac{1}{||\Delta||}\right)\mfrac{\delta_{i}^{4}}{||\Delta||^{2}}\right],$
$\displaystyle\left(\nabla^{4}\widetilde{K}(\Delta)\right)_{iiij}=\mfrac{25\sqrt{5}}{3}\sigma^{2}\phi^{5}e^{-\sqrt{5}\phi||\Delta||}\left[\mfrac{\delta_{i}}{||\Delta||^{3}}-\mfrac{1}{||\Delta||}\left(3-\sqrt{5}\phi\mfrac{\delta_{i}^{2}}{||\Delta||}\right)\right]\delta_{i}^{2}\delta_{j},$
$\displaystyle\left(\nabla^{4}\widetilde{K}(\Delta)\right)_{iijj}=\mfrac{25}{3}\sigma^{2}\phi^{4}e^{-\sqrt{5}\phi||\Delta||}\left[\left(1-\sqrt{5}\phi\mfrac{\delta_{i}^{2}}{||\Delta||}\right)\left(1-\sqrt{5}\phi\mfrac{\delta_{j}^{2}}{||\Delta||}\right)+\sqrt{5}\phi\mfrac{\delta_{i}^{2}\delta_{j}^{2}}{||\Delta||^{3}}\right],$
$\displaystyle\left(\nabla^{4}\widetilde{K}(\Delta)\right)_{ijkl}=-\mfrac{25\sqrt{5}}{3}\sigma^{2}\phi^{5}e^{-\sqrt{5}\phi||\Delta||}\left(\sqrt{5}\phi+\mfrac{1}{||\Delta||}\right)\mfrac{\delta_{i}\delta_{j}\delta_{k}\delta_{l}}{||\Delta||^{2}},$

where $i,j,k,l=1,2,\ldots,d$.

The above expressions for entries of the cross-covariance matrices correspond
to the joint process, ${\cal L}Y(\mbox{$\bf s$})=(Y(\mbox{$\bf s$}),\nabla
Y(\mbox{$\bf s$})^{\top},vech(\nabla^{2}Y(\mbox{$\bf s$}))^{\top})^{\top}$
with respect to our kernel choices. The following cross-covariance matrices
are evaluated at $||\Delta||\to 0$.

  1. 1.

Squared Exponential: In $\mathbb{R}^{d}$, we have,

$\displaystyle V_{{\cal L}Y}({\bf 0})=\sigma^{2}\begin{pmatrix}1&{\bf
0}^{\top}&-2\phi vech(I_{d})^{\top}\\\ {\bf 0}&2\phi I_{d}&{\bf O}\\\ -2\phi
vech(I_{d})&{\bf
O}&4\phi^{2}\mathrm{diag}\\{3,1\ldots,1,3,1,\ldots,1,\ldots,3\\}\end{pmatrix},$

  2. 2.

Mateŕn $(\nu=3/2)$: For this kernel the existence of only the gradient process
is guaranteed, therefore the covariance is for the process, ${\cal
L}Y(\mbox{$\bf s$})=(Y(\mbox{$\bf s$}),\nabla Y(\mbox{$\bf
s$})^{\top})^{\top}$. In $\mathbb{R}^{d}$ we have, $V_{{\cal L}Y}({\bf
0})=\sigma^{2}\begin{pmatrix}1&{\bf 0}^{\top}\\\ {\bf
0}&3\phi^{2}I_{d}\end{pmatrix}$.

  3. 3.

Mateŕn $(\nu=5/2)$: In $\mathbb{R}^{d}$ we have,

$\displaystyle V_{{\cal L}Y}({\bf 0})=\sigma^{2}\begin{pmatrix}1&{\bf
0}^{\top}&-\frac{5\phi^{2}}{3}vech(I_{d})^{\top}\\\ {\bf
0}&\frac{5\phi^{2}}{3}I_{d}&{\bf O}\\\ -\frac{5\phi^{2}}{3}vech(I_{d})&{\bf
O}&\frac{25\phi^{4}}{3}\mathrm{diag}\\{3,1\ldots,1,3,1,\ldots,1,\ldots,3\\}\end{pmatrix}.$

##  11 Algorithms

In what follows, we provide the required algorithms for sampling gradients and
wombling measures. Although listed separately to highlight the requirement of
only posterior samples, the required steps could be included within the MCMC
subroutine devised for spatial learning (or fitting the model) of
$Y(\mbox{$\bf s$})$.

#### Sampling Gradients and Curvature:

The choice for $K$ varies between Gaussian, Matérn with $\nu=3/2$ and
$\nu=5/2$. There is scope for parallel computation across grid locations.
Additionally if the inverse of estimated covariance matrices are stored for
the MCMC runs from the model fit, sufficient gains in compilation can be
achieved while sampling gradients. If $(\nu=3/2)$ is chosen the $\nabla^{2}$
terms are not computed.

Input: ${\cal S}$, A Grid ${\cal G}$ spanning ${\cal S}$, posterior MCMC
samples ${\beta}$, $\mbox{\boldmath$\theta$}_{K}=\\{\sigma^{2},\phi\\}$, $\bf
Z$

Result: Posterior samples for gradients $\nabla Y(\mbox{$\bf s$}_{g})$ and
curvature $\nabla^{2}Y(\mbox{$\bf s$}_{g})$ for $\mbox{$\bf s$}_{g}\in{\cal
G}$

for  _$i=1,2,\ldots,L$_ do

for  _$j=1,2,\ldots,n_{G}$_ do

$\Delta[i,j]=\mbox{$\bf s$}_{g}[j]-\mbox{$\bf s$}[i]$ $\triangleright$ Compute
distances of grid locations to observed process

for  _$i=1,2,\ldots,n_{\rm MCMC}$_ do

$K[i]=K(\cdot;\mbox{\boldmath$\theta$}_{K}[i])$

$K.inv[i]=(K(\cdot;\mbox{\boldmath$\theta$}_{K}[i]))^{-1}$

for  _$j=1,2,\ldots,n_{G}$_ do

$\nabla K[i,j]=(\nabla
K(\Delta[,j];\mbox{\boldmath$\theta$}_{K}[i])^{\top},vech(\nabla^{2}K(\Delta[,j];\mbox{\boldmath$\theta$}_{K}[i]))^{\top})^{\top}$

$V[i,j]=V_{{\cal L}Y}({\bf 0})$

$\mu[i,j]=\nabla\mu(\mbox{$\bf s$};\mbox{\boldmath${\beta}$}[i])-\nabla
K[i,j]^{\top}K.inv[i]\mbox{$\bf Z$}[i]$ $\triangleright$ *[r]$\mu(\mbox{$\bf
s$};\mbox{\boldmath${\beta}$}[i])=\mbox{$\bf X$}\mbox{\boldmath${\beta}$}[i]$
$\Sigma[i,j]=V[i,j]-\nabla K[i,j]^{\top}K.inv[i]\nabla K[i,j]$

${\cal L}Y[i,j]={\cal N}(\mu[i,j],\Sigma[i,j])$ $\triangleright$ Posterior
sample of Gradients and Curvature

return _${\cal L}Y$_

Algorithm 1 Algorithm for Sampling Gradients and Curvature

Input: ${\cal S}$, A curve $C$, posterior MCMC samples ${\beta}$,
$\mbox{\boldmath$\theta$}_{K}=\\{\sigma^{2},\phi\\}$, $\bf Z$

Result: Posterior samples for wombling measures
$\mbox{\boldmath$\Gamma$}(\widetilde{C}_{P})$.

for  _$j=1,2,\ldots,(n_{P}-1)$_ do

$t[j]=||C[j]-C[j+1]||$

$\mbox{$\bf u$}[j]=(C[j]-C[j+1])/t[j]$ $\triangleright$ Compute $\bf t$ and
$\bf U$

for  _$i=1,2,\ldots,L$_ do

for  _$j=1,2,\ldots,n_{P}$_ do

$\Delta[i,j]=C[j]-\mbox{$\bf s$}[i]$ $\triangleright$ Compute distances of
points in curve to observed process and norms:
$||\Delta||[i,j]=||\Delta[i,j]||$

for  _$i=1,2,\ldots,n_{\rm MCMC}$_ do

$K[i]=K(\cdot;\mbox{\boldmath$\theta$}_{K}[i])$

$K.inv[i]=(K(\cdot;\mbox{\boldmath$\theta$}_{K}[i]))^{-1}$

for  _$j=1,2,\ldots,n_{P}$_ do

$\nabla K[i,j]=\begin{pmatrix}{\tt
q}_{1}(0,t[j],D^{(1)}K(\Delta[,j]+t\mbox{$\bf
u$}[j];\mbox{\boldmath$\theta$}_{K}[i]))\\\ {\tt
q}_{1}(0,t[j],D^{(2)}K(\Delta[,j]+t\mbox{$\bf
u$}[j];\mbox{\boldmath$\theta$}_{K}[i]))\end{pmatrix}^{\top}$,

$V[i,j]=\begin{pmatrix}{\tt
q}_{2}(0,0,t[j],t[j],k_{11}(\mbox{\boldmath$\theta$}_{K}[i]))&{\tt
q}_{2}(0,0,t[j],t[j],k_{12}(\mbox{\boldmath$\theta$}_{K}[i]))\\\ {\tt
q}_{2}(0,0,t[j],t[j],k_{21}(\mbox{\boldmath$\theta$}_{K}[i]))&{\tt
q}_{2}(0,0,t[j],t[j],k_{22}(\mbox{\boldmath$\theta$}_{K}[i]))\end{pmatrix}$

$\mu[i,j]=\mu_{\mbox{\boldmath$\Gamma$}}(t[j])-\nabla
K[i,j]^{\top}K.inv[i]\mbox{$\bf Z$}[i]$ $\Sigma[i,j]=V[i,j]-\nabla
K[i,j]^{\top}K.inv[i]\nabla K[i,j]$

$\mbox{\boldmath$\Gamma$}[i,j]={\cal N}(\mu[i,j],\Sigma[i,j])$
$\triangleright$ Posterior sample of Wombling Measures

return _$\Gamma$_

Algorithm 2 Algorithm for Sampling Wombling Measures

#### Sampling Wombling Measures

The choices for $K$ are again between Gaussian, Matérn with $\nu=3/2$ and
$\nu=5/2$. Choices for curves to be evaluated for wombling boundaries range
from those outlined in Section 5.2. In case $\nu=3/2$ wombling measures for
curvature are not computed. Choices for approximations include computing
Riemann sums replacing quadrature for line integrals. There is scope for
parallel computation with the curve being broken into segments evaluated in
parallel for wombling boundaries.

The functions ${\tt q}_{1}$ and ${\tt q}_{2}$ denote one and two-dimensional
quadrature respectively. In case a Riemann sum (see (15), Section 12)
approximation is chosen, the points partitioning $C$ are treated as grid
points and the algorithm for sampling gradients and curvature is used for
predictive inference on the differential process. The Riemann sums are
computed using $\bf t$ and $\bf U$ and returned.

##  12 Proofs and Discussion

For the curvature process formulated in Section 2, we aim to show that the
covariance matrix associated with the process $D^{(2)}_{\mbox{$\bf
u$},\mbox{$\bf v$}}Y(\mbox{$\bf s$})=\mbox{$\bf c$}_{\mbox{$\bf u$},\mbox{$\bf
v$}}^{\top}vech\nabla^{2}Y(\mbox{$\bf s$})$ is valid (pg. 5 last paragraph).
We obtain the expression for the covariance matrix by leveraging the
directional finite difference process $Y^{(2)}_{\mbox{$\bf u$},\mbox{$\bf
v$},h}(\mbox{$\bf s$})$. For points $\mbox{$\bf s$},\mbox{$\bf s$}^{\prime}$,
we denote $\Delta=\mbox{$\bf s$}-\mbox{$\bf s$}^{\prime}$ and $\bf u$, $\bf v$
are unit vectors specifying direction and
$\mbox{\boldmath${\delta}$}(x,y)=\Delta+x\mbox{$\bf u$}+y\mbox{$\bf v$}$ as a
map from $\mbox{$\mathbb{R}^{2}$}\to\mbox{$\mathbb{R}^{d}$}$, after
suppressing dependence on $\Delta$ and $\bf u$, $\bf v$, let
$g(x,y)=K(\Delta(x,y))=K(\Delta+x\mbox{$\bf u$}+y\mbox{$\bf v$})$ denote a map
from $\mbox{$\mathbb{R}^{2}$}\to\mbox{$\mathbb{R}$}$ we compute the covariance
,

$\displaystyle C^{(2)}_{\mbox{$\bf u$},\mbox{$\bf v$}}(\mbox{$\bf
s$},\mbox{$\bf s$}^{\prime})$ $\displaystyle=\lim_{h\to 0}\lim_{k\to
0}E\left[Y^{(2)}_{\mbox{$\bf u$},\mbox{$\bf v$},h}(\mbox{$\bf
s$})Y^{(2)}_{\mbox{$\bf u$},\mbox{$\bf v$},k}(\mbox{$\bf
s$}^{\prime})\right],$ $\displaystyle=\lim_{h\to 0}\lim_{k\to
0}\frac{1}{h^{2}k^{2}}[g(h-k,h-k)-g(h-k,h)-g(h,h-k)+g(h,h)$
$\displaystyle\hskip 85.35826pt-g(h-k,-k)+g(h-k,0)+g(h,-k)-g(h,0)$
$\displaystyle\hskip 85.35826pt-g(-k,h-k)+g(-k,h)+g(0,h-k)-g(0,h)$
$\displaystyle\hskip 85.35826pt+g(-k,-k)-g(-k,0)-g(0,-k)+g(0,0)],$
$\displaystyle=\lim_{h\to
0}\frac{g^{\prime\prime}(h,h)-g^{\prime\prime}(h,0)-g^{\prime\prime}(0,h)+g^{\prime\prime}(0,0)}{h^{2}}=g^{(iv)}(0,0).$

On repeated application of the chain rule and noting that
$\mbox{\boldmath${\delta}$}_{x}(x,y)=\mbox{$\bf u$}^{\top}$,
$\mbox{\boldmath${\delta}$}_{y}(x,y)=\mbox{$\bf v$}^{\top}$ with all other
higher order derivatives being 0 we have,

$\displaystyle g_{x}(x,y)$
$\displaystyle=\mbox{\boldmath${\delta}$}_{x}(x,y)\nabla
K(\Delta(x,y))=\mbox{$\bf u$}^{\top}\nabla K(\Delta(x,y)),$ $\displaystyle
g_{y}(x,y)$ $\displaystyle=\mbox{\boldmath${\delta}$}_{y}(x,y)\nabla
K(\Delta(x,y))=\mbox{$\bf v$}^{\top}\nabla K(\Delta(x,y)),$ $\displaystyle
g_{xx}(x,y)$
$\displaystyle=\mbox{\boldmath${\delta}$}_{x}(x,y)\nabla^{2}K(\Delta(x,y))\mbox{\boldmath${\delta}$}_{x}(x,y)^{\top}=\mbox{$\bf
u$}^{\top}\nabla^{2}K(\Delta(x,y))\mbox{$\bf u$},$ $\displaystyle g_{xy}(x,y)$
$\displaystyle=\mbox{\boldmath${\delta}$}_{x}(x,y)\nabla^{2}K(\Delta(x,y))\mbox{\boldmath${\delta}$}_{y}(x,y)^{\top}=\mbox{$\bf
u$}^{\top}\nabla^{2}K(\Delta(x,y))\mbox{$\bf v$},$ $\displaystyle g_{yy}(x,y)$
$\displaystyle=\mbox{\boldmath${\delta}$}_{y}(x,y)\nabla^{2}K(\Delta(x,y))\mbox{\boldmath${\delta}$}_{y}(x,y)^{\top}=\mbox{$\bf
v$}^{\top}\nabla^{2}K(\Delta(x,y))\mbox{$\bf v$},$ $\displaystyle g_{xxx}(x)$
$\displaystyle=\sum_{i=1}^{d}\delta_{i,xxx}(x,y)\left(\frac{\partial
K}{\partial\delta_{i}}\right)+3\sum\limits_{i,j=1}^{d}\delta_{i,xx}(x,y)\frac{\partial^{2}K}{\partial\delta_{i}\partial\delta_{j}}\delta_{j,x}(x,y)$
$\displaystyle\hskip
56.9055pt+\sum\limits_{i,j,k=1}^{d}\delta_{i,x}(x,y)\delta_{j,x}(x,y)\delta_{k,x}(x,y)\frac{\partial^{3}K}{\partial\delta_{i}\partial\delta_{j}\partial\delta_{k}},$
$\displaystyle=\sum\limits_{i,j,k=1}^{d}\delta_{i,x}(x,y)\delta_{j,x}(x,y)\delta_{k,x}(x,y)\frac{\partial^{3}K}{\partial\delta_{i}\partial\delta_{j}\partial\delta_{k}}=\mbox{$\bf
c$}_{\mbox{$\bf u$},\mbox{$\bf u$}}^{\top}\nabla^{3}K(\Delta(x,y))\mbox{$\bf
u$}.$

Similarly $g_{xxy}(x)=\mbox{$\bf c$}_{\mbox{$\bf u$},\mbox{$\bf
u$}}^{\top}\nabla^{3}K(\Delta(x,y))\mbox{$\bf v$}$, $g_{yyx}(x)=\mbox{$\bf
c$}_{\mbox{$\bf v$},\mbox{$\bf v$}}^{\top}\nabla^{3}K(\Delta(x,y))\mbox{$\bf
u$}$ and
$g_{yyy}(x)=\mbox{$\bf c$}_{\mbox{$\bf v$},\mbox{$\bf
v$}}^{\top}\nabla^{3}K(\Delta(x,y))\mbox{$\bf v$}$. Next,

$\displaystyle g_{xxxx}(x,y)$
$\displaystyle=\sum_{i=1}^{d}\delta_{i,xxxx}(x,y)\left(\frac{\partial
K}{\partial\delta_{i}}\right)+\sum\limits_{i,j=1}^{d}\delta_{i,xxx}(x,y)\frac{\partial^{2}K}{\partial\delta_{i}\partial\delta_{j}}\delta_{j,x}(x,y)$
$\displaystyle+3\bigg{[}\sum\limits_{i,j=1}^{d}\delta_{i,xxx}(x,y)\frac{\partial^{2}K}{\partial\delta_{i}\partial\delta_{j}}\delta_{j,x}(x,y)+\sum\limits_{i,j=1}^{d}\delta_{i,xx}(x,y)\frac{\partial^{2}K}{\partial\delta_{i}\partial\delta_{j}}\delta_{j,xx}(x,y)$
$\displaystyle\hskip
28.45274pt+\sum_{i,j,k=1}^{d}\delta_{i,xx}(x,y)\delta_{j,x}(x,y)\delta_{k,x}(x,y)\frac{\partial^{3}K}{\partial\delta_{i}\partial\delta_{j}\partial\delta_{k}}\bigg{]}$
$\displaystyle\hskip
14.22636pt+3\sum\limits_{i,j,k=1}^{d}\delta_{i,xx}(x,y)\delta_{j,x}(x,y)\delta_{k,x}(x,y)\frac{\partial^{3}K}{\partial\delta_{i}\partial\delta_{j}\partial\delta_{k}}$
$\displaystyle+\sum_{i,j,k,l=1}^{d}\delta_{i,x}(x,y)\delta_{j,x}(x,y)\delta_{k,x}(x,y)\delta_{l}^{\prime}(x)\frac{\partial^{4}K}{\partial\delta_{i}\partial\delta_{j}\partial\delta_{k}\partial\delta_{l}},$
$\displaystyle=\sum_{i,j,k,l=1}^{d}\delta_{i,x}(x,y)\delta_{j,x}(x,y)\delta_{k,x}(x,y)\delta_{l}^{\prime}(x)\frac{\partial^{4}K}{\partial\delta_{i}\partial\delta_{j}\partial\delta_{k}\partial\delta_{l}}=\mbox{$\bf
c$}_{\mbox{$\bf u$},\mbox{$\bf
u$}}^{\top}\nabla^{4}K(\mbox{\boldmath${\delta}$}(x))\mbox{$\bf
c$}_{\mbox{$\bf u$},\mbox{$\bf u$}}.$

Similarly, $g_{xxxy}(x,y)=\mbox{$\bf c$}_{\mbox{$\bf u$},\mbox{$\bf
u$}}^{\top}\nabla^{4}K(\mbox{\boldmath${\delta}$}(x))\mbox{$\bf
c$}_{\mbox{$\bf u$},\mbox{$\bf v$}}$, $g_{yyyx}(x,y)=\mbox{$\bf
c$}_{\mbox{$\bf v$},\mbox{$\bf
v$}}^{\top}\nabla^{4}K(\mbox{\boldmath${\delta}$}(x))\mbox{$\bf
c$}_{\mbox{$\bf v$},\mbox{$\bf u$}}$, $g_{xxyy}(x,y)=\mbox{$\bf
c$}_{\mbox{$\bf u$},\mbox{$\bf
v$}}^{\top}\nabla^{4}K(\mbox{\boldmath${\delta}$}(x))\mbox{$\bf
c$}_{\mbox{$\bf u$},\mbox{$\bf v$}}$ and $g_{yyyy}(x,y)=\mbox{$\bf
c$}_{\mbox{$\bf v$},\mbox{$\bf
v$}}^{\top}\nabla^{4}K(\mbox{\boldmath${\delta}$}(x))\mbox{$\bf
c$}_{\mbox{$\bf v$},\mbox{$\bf v$}}$. Evaluated at $x,y=0$, i.e.
$\mbox{\boldmath${\delta}$}(0,0)=\Delta$, $g_{xxyy}$ provides us with the
required expression, $C^{(2)}_{\mbox{$\bf u$},\mbox{$\bf v$}}(\mbox{$\bf
s$},\mbox{$\bf s$}^{\prime})=\mbox{$\bf c$}_{\mbox{$\bf u$},\mbox{$\bf
v$}}^{\top}\nabla^{4}K(\Delta)\mbox{$\bf c$}_{\mbox{$\bf u$},\mbox{$\bf v$}}$
and $var(D^{(2)}_{\mbox{$\bf u$},\mbox{$\bf v$}}Y(\mbox{$\bf s$}))=\lim_{h\to
0}E(Y^{(2)}_{\mbox{$\bf u$},\mbox{$\bf v$},h}(\mbox{$\bf
s$}),Y^{(2)}_{\mbox{$\bf u$},\mbox{$\bf v$},k}(\mbox{$\bf s$}))=\mbox{$\bf
c$}_{\mbox{$\bf u$},\mbox{$\bf v$}}^{\top}\nabla^{4}K({\bf 0})\mbox{$\bf
c$}_{\mbox{$\bf u$},\mbox{$\bf v$}}$ which exists if $K^{(iv)}(\Delta)$ exists
for all $\Delta$, including $\Delta={\bf 0}$.

#### Note:

In the above proof we make some abuse of notation for brevity of mathematical
expressions involved. To clarify, $\delta_{x}=\frac{\partial}{\partial
x}\delta(x,y)$, $\delta_{xx}=\mfrac{\partial^{2}}{\partial x^{2}}\delta(x,y)$
and so on, $g_{x}=\frac{\partial}{\partial x}g(x,y)$,
$g_{xx}=\mfrac{\partial^{2}}{\partial x^{2}}g(x,y)$ etc.,
$\sum_{i,j=1}^{d}=\sum_{i=1}^{d}\sum_{j=1}^{d}$ etc.

To derive (2) and the covariance for the directional curvature process, we
assume that $Y(\mbox{$\bf s$})$ is isotropic, i.e.
$K(\Delta)=\widetilde{K}(||\Delta||)$ therefore,

$\displaystyle C^{(2)}_{\medmath\mbox{$\bf u$},\medmath\mbox{$\bf
u$}}(\mbox{$\bf s$},\mbox{$\bf s$}^{\prime})$ $\displaystyle=\lim_{h\to
0}\lim_{k\to 0}E\left[Y^{(2)}_{\mbox{$\bf u$},\mbox{$\bf v$},h}(\mbox{$\bf
s$})Y^{(2)}_{\mbox{$\bf u$},\mbox{$\bf v$},k}(\mbox{$\bf
s$}^{\prime})\right],$ $\displaystyle=\lim_{h\to 0}\lim_{k\to
0}\frac{1}{h^{2}k^{2}}\bigg{[}E\left(Y(\mbox{$\bf s$}+h(\mbox{$\bf
u$}+\mbox{$\bf v$}))Y^{(2)}_{\mbox{$\bf u$},\mbox{$\bf v$},k}(\mbox{$\bf
s$}^{\prime})\right)-E(Y(\mbox{$\bf s$}+h\mbox{$\bf u$})Y^{(2)}_{\mbox{$\bf
u$},\mbox{$\bf v$},k}(\mbox{$\bf s$}^{\prime}))$ $\displaystyle\hskip
85.35826pt-E(Y(\mbox{$\bf s$}+k\mbox{$\bf v$})Y^{(2)}_{\mbox{$\bf
u$},\mbox{$\bf v$},k}(\mbox{$\bf s$}^{\prime}))+E(Y(\mbox{$\bf
s$})Y^{(2)}_{\mbox{$\bf u$},\mbox{$\bf v$},k}(\mbox{$\bf
s$}^{\prime}))\bigg{]}$

where,

$\displaystyle E\left(Y(\mbox{$\bf s$}+h(\mbox{$\bf u$}+\mbox{$\bf
v$}))Y^{(2)}_{\mbox{$\bf u$},\mbox{$\bf v$},k}(\mbox{$\bf
s$}^{\prime})\right)$ $\displaystyle=\widetilde{K}(||\Delta+(h-k)(\mbox{$\bf
u$}+\mbox{$\bf v$})||)-\widetilde{K}(||\Delta+(h-k)\mbox{$\bf u$}+h\mbox{$\bf
v$}||)-$ $\displaystyle~{}~{}\widetilde{K}(||\Delta+h\mbox{$\bf
u$}+(h-k)\mbox{$\bf v$}||)+\widetilde{K}(||\Delta+h(\mbox{$\bf u$}+\mbox{$\bf
v$})||)$ $\displaystyle~{}~{}-\widetilde{K}(||\Delta+(h-k)\mbox{$\bf
u$}-k\mbox{$\bf v$}||),$ $\displaystyle E(Y(\mbox{$\bf s$}+h\mbox{$\bf
u$})Y^{(2)}_{\mbox{$\bf u$},\mbox{$\bf v$},k}(\mbox{$\bf s$}^{\prime}))$
$\displaystyle=\widetilde{K}(||\Delta+(h-k)\mbox{$\bf u$}-k\mbox{$\bf
v$}||)-\widetilde{K}(||\Delta+(h-k)\mbox{$\bf u$}||)$
$\displaystyle~{}~{}-\widetilde{K}(||\Delta+h\mbox{$\bf u$}-k\mbox{$\bf
v$}||)+\widetilde{K}(||\Delta+h\mbox{$\bf u$}||)$ $\displaystyle
E(Y(\mbox{$\bf s$}+k\mbox{$\bf v$})Y^{(2)}_{\mbox{$\bf u$},\mbox{$\bf
v$},k}(\mbox{$\bf s$}^{\prime}))$
$\displaystyle=\widetilde{K}(||\Delta-k\mbox{$\bf u$}+(h-k)\mbox{$\bf
v$}||)-\widetilde{K}(||\Delta-k\mbox{$\bf u$}+h\mbox{$\bf v$}||)$
$\displaystyle~{}~{}-\widetilde{K}(||\Delta+(h-k)\mbox{$\bf
v$}||)+\widetilde{K}(||\Delta+h\mbox{$\bf v$}||)$ $\displaystyle
E(Y(\mbox{$\bf s$})Y^{(2)}_{\mbox{$\bf u$},\mbox{$\bf v$},k}(\mbox{$\bf
s$}^{\prime}))$ $\displaystyle=\widetilde{K}(||\Delta-k(\mbox{$\bf
u$}+\mbox{$\bf v$})||)-\widetilde{K}(||\Delta-k\mbox{$\bf
u$}||)-\widetilde{K}(||\Delta-k\mbox{$\bf v$}||)+\widetilde{K}(||\Delta||)$

suppressing dependence on $\Delta$, $\bf u$ and $\bf v$ , we define
$\rho(h,k)=||\Delta(h,k)||=||\Delta+h\mbox{$\bf u$}+k\mbox{$\bf v$}||$ and let
$g(h,k)=K(\rho(h,k))$. Hence,

$\displaystyle C^{(2)}_{\mbox{$\bf u$},\mbox{$\bf v$}}(\mbox{$\bf
s$},\mbox{$\bf s$}^{\prime})$ $\displaystyle=\lim_{h\to 0}\lim_{k\to
0}\frac{1}{h^{2}k^{2}}[\rho(h-k,h-k)-\rho(h-k,h)-\rho(h,h-k)+\rho(h,h)$
$\displaystyle\hskip 85.35826pt-\rho(h-k,-k)+\rho(h-k,0)+\rho(h,-k)-\rho(h,0)$
$\displaystyle\hskip 85.35826pt-\rho(-k,h-k)+\rho(-k,h)+\rho(0,h-k)-\rho(0,h)$
$\displaystyle\hskip 85.35826pt+\rho(-k,-k)-\rho(-k,0)-\rho(0,-k)+\rho(0,0)],$
$\displaystyle=\lim_{h\to
0}\frac{\rho^{\prime\prime}(h,h)-\rho^{\prime\prime}(h,0)-\rho^{\prime\prime}(0,h)+\rho^{\prime\prime}(0,0)}{h^{2}}=\rho^{(iv)}(0,0)\;.$

Since, $\rho(0,0)=||\Delta||$, from the previous proof we can see that,
$\rho^{(iv)}(0,0)=\nabla\widetilde{K}(||\Delta||)\rho^{(iv)}(0,0)+\nabla^{2}\widetilde{K}(||\Delta||)\left(4\rho^{\prime\prime\prime}(0,0)\rho^{\prime}(0,0)+3[\rho^{\prime\prime}(0,0)]^{2}\right)+6\nabla^{3}\widetilde{K}(||\Delta||)\rho^{\prime\prime}(0,0)[\rho^{\prime}(0,0)]^{2}+\nabla^{4}\widetilde{K}(||\Delta||)[\rho^{\prime}(0,0)]^{4}$.
After some algebra, $\rho^{\prime}(0)=\frac{\mbox{$\bf
u$}^{\top}\Delta}{||\Delta||}$,
$\rho^{\prime\prime}(0)=\frac{1}{||\Delta||}\left(1-\frac{(\mbox{$\bf
u$}^{\top}\Delta)^{2}}{||\Delta||^{2}}\right)$,
$\rho^{\prime\prime\prime}(0)=-3\frac{\mbox{$\bf
u$}^{\top}\Delta}{||\Delta||^{3}}\left(1-\frac{(\mbox{$\bf
u$}^{\top}\Delta)^{2}}{||\Delta||^{2}}\right)$ and
$\rho^{(iv)}(0)=\frac{3}{||\Delta||^{3}}\left(5\frac{(\mbox{$\bf
u$}^{\top}\Delta)^{2}}{||\Delta||^{2}}-1\right)\left(1-\frac{(\mbox{$\bf
u$}^{\top}\Delta)^{2}}{||\Delta||^{2}}\right)$. Substituting and grouping
terms corresponding to
$\left(\nabla^{2}\widetilde{K}(||\Delta||)-\frac{\nabla\widetilde{K}(||\Delta||)}{||\Delta||}\right)$,
$\nabla^{3}\widetilde{K}(||\Delta||)$ and
$\nabla^{4}\widetilde{K}(||\Delta||)$ we get,

$\begin{split}g^{(iv)}(0)&=\frac{3}{||\Delta||^{2}}\left\\{1-5\frac{(\mbox{$\bf
u$}^{\top}\Delta)^{2}}{||\Delta||^{2}}\right\\}\left[1-\frac{(\mbox{$\bf
u$}^{\top}\Delta)^{2}}{||\Delta||^{2}}\right]\left(\nabla^{2}\widetilde{K}(||\Delta||)-\frac{\nabla\widetilde{K}(||\Delta||)}{||\Delta||}\right)\\\
&\hskip 14.22636pt+\frac{6}{||\Delta||}\frac{(\mbox{$\bf
u$}^{\top}\Delta)^{2}}{||\Delta||^{2}}\left[1-\frac{(\mbox{$\bf
u$}^{\top}\Delta)^{2}}{||\Delta||^{2}}\right]\nabla^{3}\widetilde{K}(||\Delta||)+\left(\frac{(\mbox{$\bf
u$}^{\top}\Delta)^{2}}{||\Delta||^{2}}\right)^{2}\nabla^{4}\widetilde{K}(||\Delta||),\end{split}$

which is the required expression. We discuss validity of the curvature process
for surfaces in $\mathbb{R}^{3}$. It can be easily extended to surfaces in
$\mathbb{R}^{d}$. Define,

$\displaystyle{\cal L}_{\mbox{$\bf e$}_{1},\mbox{$\bf e$}_{2},h}Y(\mbox{$\bf
s$})$ $\displaystyle=\begin{pmatrix}1&0&0&0&0&0\\\
-\frac{1}{h}&\frac{1}{h}&0&0&0&0\\\ -\frac{1}{h}&0&\frac{1}{h}&0&0&0\\\
\frac{1}{h^{2}}&-\frac{2}{h^{2}}&0&\frac{1}{h^{2}}&0&0\\\
\frac{1}{h^{2}}&-\frac{1}{h^{2}}&-\frac{1}{h^{2}}&0&\frac{1}{h^{2}}&0\\\
\frac{1}{h^{2}}&0&-\frac{2}{h^{2}}&0&0&\frac{1}{h^{2}}\end{pmatrix}\begin{pmatrix}Y(\mbox{$\bf
s$})\\\ Y(\mbox{$\bf s$}+h\mbox{$\bf e$}_{1})\\\ Y(\mbox{$\bf s$}+h\mbox{$\bf
e$}_{2})\\\ Y(\mbox{$\bf s$}+2h\mbox{$\bf e$}_{1})\\\ Y(\mbox{$\bf
s$}+h(\mbox{$\bf e$}_{1}+\mbox{$\bf e$}_{2}))\\\ Y(\mbox{$\bf s$}+2h\mbox{$\bf
e$}_{2})\end{pmatrix}$ $\displaystyle=\mbox{\boldmath$A$}_{h}\mbox{$\bf
L$}_{h}(\mbox{$\bf s$})=\left(Y(\mbox{$\bf s$}),Y^{(1)}_{\mbox{$\bf
e$}_{1},h}(\mbox{$\bf s$}),Y^{(1)}_{\mbox{$\bf e$}_{2},h}(\mbox{$\bf
s$}),Y^{(2)}_{\mbox{$\bf e$}_{1},\mbox{$\bf e$}_{1},h}(\mbox{$\bf
s$}),Y^{(2)}_{\mbox{$\bf e$}_{1},\mbox{$\bf e$}_{2},h}(\mbox{$\bf
s$}),Y^{(2)}_{\mbox{$\bf e$}_{2},\mbox{$\bf e$}_{2},h}(\mbox{$\bf
s$})\right)^{\top},$

as finite difference corresponding to the differential operator ${\cal L}Y$
using the expressions of $Y^{(1)}_{\mbox{$\bf e$}_{i},h}(\mbox{$\bf s$})$ and
$Y^{(2)}_{\mbox{$\bf e$}_{i},\mbox{$\bf e$}_{j},h}(\mbox{$\bf s$})$,
$i,j=1,2$. For every $h>0$ this defines a linear transformation, since the
determinant of $\mbox{\boldmath$A$}_{h}$ is $h^{-8}$. We denote the
differenced differential process on the right of $\mbox{\boldmath$A$}_{h}$,
suppressing dependence on $\mbox{$\bf e$}_{1},\mbox{$\bf e$}_{2}$ by
$\mbox{$\bf L$}_{h}(\mbox{$\bf s$})$. The associated covariance matrix is
given by, $Cov({\cal L}_{\mbox{$\bf e$}_{1},\mbox{$\bf e$}_{2}h}Y(\mbox{$\bf
s$}),{\cal L}_{\mbox{$\bf e$}_{1},\mbox{$\bf e$}_{2},k}Y(\mbox{$\bf
s$}^{\prime}))=\mbox{\boldmath$A$}_{h}{\cal K}_{\mbox{$\bf e$}_{1},\mbox{$\bf
e$}_{2},h,k}(\Delta)\mbox{\boldmath$A$}_{k}^{\top}$, where elements of ${\cal
K}_{\mbox{$\bf e$}_{1},\mbox{$\bf e$}_{2},h,k}(\Delta)$ are obtained from
$Cov(\mbox{$\bf L$}_{h}(\mbox{$\bf s$}),\mbox{$\bf L$}_{k}(\mbox{$\bf
s$}^{\prime}))$. Hence, as $h\downarrow 0$, ${\cal L}_{\mbox{$\bf
e$}_{1},\mbox{$\bf e$}_{2},h}Y(\mbox{$\bf s$})\to{\cal L}Y(\mbox{$\bf s$})$
and $\lim_{h\downarrow 0,k\downarrow 0}Cov({\cal L}_{\mbox{$\bf
e$}_{1},\mbox{$\bf e$}_{2}h}Y(\mbox{$\bf s$}),{\cal L}_{\mbox{$\bf
e$}_{1},\mbox{$\bf e$}_{2},k}Y(\mbox{$\bf s$}^{\prime}))=V_{{\cal
L}Y}(\Delta)$, where the limits operate element-wise on the matrix
$\mbox{\boldmath$A$}_{h}{\cal K}_{\mbox{$\bf e$}_{1},\mbox{$\bf
e$}_{2},h,k}(\Delta)\mbox{\boldmath$A$}_{k}^{\top}$ and, the expression for
each element is obtained from previous computations, by setting $\mbox{$\bf
u$}=\mbox{$\bf e$}_{1}$ and $\mbox{$\bf v$}=\mbox{$\bf e$}_{2}$.

For the directional operator, this can be established by observing that the
directional differential operator is obtained as follows,

$\displaystyle{\cal L}_{\mbox{$\bf u$},h}Y(\mbox{$\bf
s$})=\begin{pmatrix}1&0&0&0&0&0\\\ 0&u_{1}&u_{2}&0&0&0\\\
0&0&0&u_{1}^{2}&2u_{1}u_{2}&u_{2}^{2}\end{pmatrix}{\cal L}_{\mbox{$\bf
e$}_{1},\mbox{$\bf e$}_{2},h}Y(\mbox{$\bf s$})=\left(1\oplus\mbox{$\bf
u$}^{\top}\oplus\mbox{$\bf c$}_{\mbox{$\bf u$},\mbox{$\bf
u$}}^{\top}\right){\cal L}_{\mbox{$\bf e$}_{1},\mbox{$\bf
e$}_{2},h}Y(\mbox{$\bf s$}),$

where $\oplus$ denotes the direct sum for matrices, then as $h\to 0$, ${\cal
L}_{\mbox{$\bf u$},h}Y(\mbox{$\bf s$})\to{\cal L}_{\mbox{$\bf u$}}Y(\mbox{$\bf
s$})$. The covariance matrix is obtained following similar arguments presented
in the proof for the previous result. In case the covariance is isotropic we
have $K(\Delta)=\widetilde{K}(||\Delta||)$, on repeated differentiation and
noting that
$\frac{\partial}{\partial\Delta}||\Delta||=\frac{\Delta}{||\Delta||}$ we have,

$\displaystyle\nabla K(\Delta)$
$\displaystyle=\frac{\nabla\widetilde{K}(||\Delta||)}{||\Delta||}\Delta,$
$\displaystyle\nabla^{2}K(\Delta)$
$\displaystyle=\frac{\nabla\widetilde{K}(||\Delta||)}{||\Delta||}I-\frac{\nabla\widetilde{K}(||\Delta||)}{||\Delta||^{3}}\Delta\Delta^{\top}+\frac{\nabla^{2}\widetilde{K}(||\Delta||)}{||\Delta||^{2}}\Delta\Delta^{\top}$
$\displaystyle=\frac{\nabla\widetilde{K}(||\Delta||)}{||\Delta||}I+\left(\nabla^{2}\widetilde{K}(||\Delta||)-\frac{\nabla\widetilde{K}(||\Delta||)}{||\Delta||}\right)\frac{\Delta\Delta^{\top}}{||\Delta||^{2}}.$

Differentiating $\nabla^{2}K(\Delta)$ w.r.t. $\Delta$ we obtain,

$\displaystyle\nabla^{3}K(\Delta)$
$\displaystyle=\frac{\nabla^{2}\widetilde{K}(||\Delta||)}{||\Delta||^{2}}vech(I)^{\top}\otimes\Delta-\frac{\nabla\widetilde{K}(||\Delta||)}{||\Delta||^{3}}vech(I)^{\top}\otimes\Delta$
$\displaystyle~{}~{}+\left(\nabla^{2}\widetilde{K}(||\Delta||)-\frac{\nabla\widetilde{K}(||\Delta||)}{||\Delta||}\right)\frac{1}{||\Delta||^{2}}\frac{\partial
vech(\Delta\Delta^{\top})}{\partial\Delta}$
$\displaystyle~{}~{}-2\left(\nabla^{2}\widetilde{K}(||\Delta||)-\frac{\nabla\widetilde{K}(||\Delta||)}{||\Delta||}\right)\frac{1}{||\Delta||^{4}}\frac{\partial
vech(\Delta\Delta^{\top})}{\partial\Delta}$
$\displaystyle~{}~{}-\left(\nabla^{2}\widetilde{K}(||\Delta||)-\frac{\nabla\widetilde{K}(||\Delta||)}{||\Delta||}\right)\frac{1}{||\Delta||^{4}}\frac{\partial
vech(\Delta\Delta^{\top})}{\partial\Delta}$
$\displaystyle~{}~{}+\nabla^{3}\widetilde{K}(||\Delta||)\cdot\frac{vech(\Delta\Delta^{\top})^{\top}\otimes\Delta}{||\Delta||^{3}}.$

On grouping terms for
$\nabla^{2}\widetilde{K}(||\Delta||)-\frac{\nabla\widetilde{K}(||\Delta||)}{||\Delta||}$
and $\nabla^{3}\widetilde{K}(||\Delta||)$ we obtain the required expression,

$\displaystyle\nabla^{3}K(\Delta)$
$\displaystyle=\left(\nabla^{2}\widetilde{K}(||\Delta||)-\frac{\nabla\widetilde{K}(||\Delta||)}{||\Delta||}\right)\Bigg{\\{}\frac{vech(I)^{\top}\otimes\Delta}{||\Delta||^{2}}-3\frac{vech(\Delta\Delta^{\top})^{\top}\otimes\Delta}{||\Delta||^{4}}$
$\displaystyle\hskip 199.16928pt+\frac{1}{||\Delta||^{2}}\frac{\partial
vech(\Delta\Delta^{\top})}{\partial\Delta}\Bigg{\\}}$ $\displaystyle\hskip
85.35826pt+\nabla^{3}\widetilde{K}(||\Delta||)\cdot\frac{vech(\Delta\Delta^{\top})^{\top}\otimes\Delta}{||\Delta||^{3}}$

To obtain $\nabla^{4}K(\Delta)$ we differentiate $\nabla^{3}K(\Delta)$ w.r.t.
$\Delta$, we use notations $A_{1}=\frac{\partial\Delta\otimes
vech(I)^{\top}}{\partial\Delta}$, $A_{2}=\frac{\partial\Delta\otimes
vech(\Delta\Delta^{\top})^{\top}}{\partial\Delta}$,
$A_{3}=\frac{\partial}{\partial\Delta}\left(\frac{\partial
vech(\Delta\Delta^{\top})}{\partial\Delta}\right)$ for matricized tensors of
order $d(d+1)/2\times d(d+1)/2$, where the order of matricization conforms to
the listing order of the half-vectorization operator $vech$ and $A_{4}$
denotes the element-wise product of $\Delta$ with $\left(\frac{\partial
vech(\Delta\Delta^{\top})}{\partial\Delta}\right)$ in the same order as the
matricized tensor.

On differentiating the factor corresponding to the coefficient
$\nabla^{2}\widetilde{K}(||\Delta||)-\frac{\nabla\widetilde{K}(||\Delta||)}{||\Delta||}$
we obtain,

$\displaystyle\begin{aligned}
\frac{1}{||\Delta||^{2}}A_{1}-3\frac{1}{||\Delta||^{4}}A_{2}+\frac{1}{||\Delta||^{2}}A_{3}-\frac{2}{||\Delta||^{4}}vech(\Delta\Delta^{\top})vech(I)^{\top}\\\
+\frac{12}{||\Delta||^{6}}vech(\Delta\Delta^{\top})vech(\Delta\Delta^{\top})^{\top}-\frac{2}{||\Delta||^{4}}A_{4}.\end{aligned}$

Differentiating
$\nabla^{2}\widetilde{K}(||\Delta||)-\frac{\nabla\widetilde{K}(||\Delta||)}{||\Delta||}$,

$\displaystyle\frac{\nabla^{3}\widetilde{K}(||\Delta||)}{||\Delta||}\Delta-\left(\nabla^{2}\widetilde{K}(||\Delta||)-\frac{\nabla\widetilde{K}(||\Delta||)}{||\Delta||}\right)\frac{\Delta\Delta^{\top}}{||\Delta||^{2}}.$

Differentiating the factor corresponding to the coefficient
$\nabla^{3}\widetilde{K}(||\Delta||)$ we obtain,

$\displaystyle-3\frac{vech(\Delta\Delta^{\top})vech(\Delta\Delta^{\top})^{\top}}{||\Delta||^{5}}+\frac{1}{||\Delta||^{3}}A_{2},$

finally, differentiating
$\nabla^{3}\widetilde{K}(||\Delta||)=\nabla^{4}\widetilde{K}(||\Delta||)\frac{\Delta}{||\Delta||}$.
Grouping coefficients for
$\nabla^{2}\widetilde{K}(||\Delta||)-\frac{\nabla\widetilde{K}(||\Delta||)}{||\Delta||}$,
$\nabla^{3}\widetilde{K}(||\Delta||)$ and
$\nabla^{4}\widetilde{K}(||\Delta||)$ we obtain the required expression,
thereby completing the proof.

For proving the result focusing on spectral theory, note that since $f$ is
symmetric about 0 by hypothesis,

$\displaystyle K(t)=\int_{\mbox{$\mathbb{R}$}}e^{i\lambda
t}f(\lambda)d\lambda=\int_{\mbox{$\mathbb{R}$}}\cos(\lambda
t)f(\lambda)d\lambda+i\int_{\mbox{$\mathbb{R}$}}\sin(\lambda
t)f(\lambda)d\lambda=\int\cos(\lambda t)f(\lambda)d\lambda.$

Differentiating w.r.t. $t$ on both sides we have,

$\displaystyle\nabla K(t)=-\int\sin(\lambda t)\lambda f(\lambda)d\lambda,$

Since $|\sin(\lambda t)\lambda f(\lambda)|\leq|\lambda|f(\lambda)$ and
$\int|\lambda|f(\lambda)d\lambda<\infty$ under hypothesis, differentiation
under the integral sign is valid. We repeat the process to obtain,

$\displaystyle\nabla^{2}K(t)=-\int\cos(\lambda
t)\lambda^{2}f(\lambda)d\lambda$ ,
$\displaystyle\nabla^{3}K(t)=\int\sin(\lambda t)\lambda^{3}f(\lambda)d\lambda$
, $\displaystyle\nabla^{4}K(t)=\int\cos(\lambda
t)\lambda^{4}f(\lambda)d\lambda,$

Next we make the following observations for limits of these derivatives,

$\displaystyle\lim_{t\to 0}\left(\nabla^{2}K(t)-\frac{\nabla
K(t)}{t}\right)=\lim_{t\to 0}\int\left(\cos(\lambda t)-\frac{\sin(\lambda
t)}{\lambda t}\right)\lambda^{2}f(\lambda)d\lambda=0,$
$\displaystyle\lim_{t\to 0}\nabla^{2}K(t)=-\lim_{t\to 0}\int\cos(\lambda
t)\lambda^{2}f(\lambda)d\lambda=-\int\lambda^{2}f(\lambda)d\lambda<\infty,$
$\displaystyle\lim_{t\to 0}\nabla^{3}K(t)=\lim_{t\to 0}\int\sin(\lambda
t)\lambda^{3}f(\lambda)d\lambda=0,$ $\displaystyle\lim_{t\to
0}\nabla^{4}K(t)=\lim_{t\to 0}\int\cos(\lambda
t)\lambda^{4}f(\lambda)d\lambda=\int\lambda^{4}f(\lambda)d\lambda<\infty,$

We evaluate the results obtained above under these observations. Making note
of,

$\displaystyle\left\\{\frac{vech(I)^{\top}\otimes\Delta}{||\Delta||^{2}}-3\frac{vech(\Delta\Delta^{\top})^{\top}\otimes\Delta}{||\Delta||^{4}}+\frac{1}{||\Delta||^{2}}\frac{\partial
vech(\Delta\Delta^{\top})}{\partial\Delta}\right\\}$ ,
$\displaystyle\frac{vech(\Delta\Delta^{\top})^{\top}\otimes\Delta}{||\Delta||^{3}}$

stay bounded as $\Delta\to 0$ implying $\nabla^{3}K(\Delta)\to 0$ as
$\Delta\to 0$. For $\nabla^{4}K(\Delta)$, we observe that the factors
corresponding to
$\nabla^{2}\widetilde{K}(||\Delta||)-\frac{\nabla\widetilde{K}(||\Delta||)}{||\Delta||}$,
$\nabla^{3}\widetilde{K}(||\Delta||)$ and