Covariance functions motivated by spatial random field models with local interactions

Abstract

Random fields based on energy functionals with local interactions possess flexible covariance functions, lead to computationally efficient algorithms for spatial data processing, and have important applications in Bayesian field theory. In this paper we address the calculation of covariance functions for a family of isotropic local-interaction random fields in two dimensions. We derive explicit expressions for non-differentiable Spartan covariance functions in ${\mathbb{R}}^2$ that are based on the modified Bessel function of the second kind. We also derive a family of infinitely differentiable, Bessel-Lommel covariance functions that exhibit a hole effect and are valid in ${\mathbb{R}}^{d}$, where d > 2. Finally, we define a generalized spectrum of correlation scales that can be applied to both differentiable and non-differentiable random fields in contrast with the smoothness microscale.

Appendices

Appendix A: Proof of Proposition 1

Proof

Defining dimensionless wavevectors $u= k\,\xi$ and lag distances $h=r/\xi$, the spectral integral (9) is simplified as follows:

$$ C_{\mathrm{xx}}(h;{\varvec{\theta }}) = \frac{\eta _0 }{2\pi } \, \int \limits _{0}^{\infty } du \,\frac{ u \, J_{0}(u\,h) }{1+{\eta _{1}}u^2+ u^4}. $$

(27)

Equation (27) shows that the only non-trivial parameter is ${\eta _{1}}$; $\eta _0$ is a multiplicative scale factor, whereas the characteristic length $\xi$ is absorbed in the non-dimensional lag $h.$ The rational function $1/\varPi (u)$, where $\varPi (u)$ is the SSRF characteristic polynomial defined in (7b), admits the following expansion

$$ \frac{1}{\varPi (u)}= \left\{ \begin{array}{cc} \frac{1}{t_{+}^{*}-t_{-}^{*}} \, \left( \frac{1}{u^2-t_{+}^{*}} - \frac{1}{u^2 -t_{-}^{*}} \right) , &{} {\eta _{1}}\ne 2, \\ \frac{1}{(u^2+1)^2} &{} {\eta _{1}}= 2, \end{array} \right. $$

(28)

where $t_{\pm }^{*}={\left( -{\eta _{1}}\pm \Delta \right) }/{2}$ are the roots of $\varPi (t=u^2)$.

In light of (28), the integral (27) is evaluated using the Hankel–Nicholson formula (11.4.44) in (Abramowitz and Stegun (1972), p. 364):

$$\int _{0}^{\infty } du \frac{u^{\nu +1}J_{\nu }(h\,u)}{(u^2 + z^2)^{\mu +1}} = \frac{h^{\mu }z^{\nu -\mu }}{2^{\mu }\Gamma (\mu +1)} K_{\nu -\mu }(h\,z). $$

(29)

This equation is valid for $h>0, {\mathfrak{R}}(z)>0$, and $-1<{\mathfrak{R}}(\nu )< 2{\mathfrak{R}}(\mu )+\frac{3}{2}$. The above is applied to (27) with (i) ${\eta _{1}}\ne 2$, $\nu =0$, $\mu =0$, $z^{2}_{\pm }= - t_{\pm }^{*}$ and (ii) ${\eta _{1}}= 2$, $\nu =0$ and $\mu =1$. In case (ii) we obtain (10b) and in case (i) the following

$$ C_{\mathrm{xx}}(h ;{\varvec{\theta }}) =\frac{\eta _0 \left[ K_{0} (hz_{+}) - K_{0} (hz_{-}) \right] }{2\pi \sqrt{\eta _{1}^{2}-4}} , \quad {\eta _{1}}\ne 2. $$

The coefficients $z_{\pm }=\sqrt{-t_{\pm }^{*}}$ are plotted versus ${\eta _{1}}$

Fig. 10

Real (a) and imaginary (b) parts of the roots $z_{+}=\sqrt{-t_{+}^{*}}$ and $z_{-}=\sqrt{-t_{-}^{*}}$ of the characteristic polynomial Π(u)  = 1 + η ₁ u ² + u ⁴ according to (11)

in Fig. 10. For ${\eta _{1}}>2$ both $z_{+}$ and $z_{-}$ are real numbers, hence proving (10a). For $-2< {\eta _{1}}<2$ ${\mathfrak{R}}(z_{+})={\mathfrak{R}}(z_{-})$, whereas ${\mathfrak{I}}(z_{+})=-{\mathfrak{I}}(z_{-})$, i.e., $z_{-} =\overline{z_{+}}$. The analytic continuation property $K_{0}(\overline{z})=\overline{K_{0}(z)}$ ((Abramowitz and Stegun 1972, p. 377)) leads to (10c) which is explicitly real-valued.

1.1 Continuity

A stationary SRF is mean square continuous $\forall {\mathbf{s}}\in {\mathbb{R}}^d$ if and only if its covariance function is continuous at zero lag (Adler 1981; Abrahamsen 1997). This condition is satisfied for the SSRF covariance.

1.2 Differentiability

Differentiability of the SRF $X({\mathbf{s}},\omega )$ in the mean-square sense requires that all second-order partial derivatives of the covariance function at $\Vert {\mathbf{r}}\Vert =0$ exist ((Adler 1981, p. 27)). This requirement is equivalent to the convergence of the second-order spectral moment

$$ \varLambda _{d}^{(2)}:=\int _{{\mathbb{R}}^d} d{\mathbf{k}}\, k^2 \, \widetilde{C_{\mathrm{xx}}}({\mathbf{k}},{\varvec{\theta }}).$$

For the SSRF spectral density in $d=2$ the above becomes

$$ \varLambda _{2}^{(2)} \propto \lim _{{k_{c}}\rightarrow \infty } \int _{0}^{{k_{c}}} dk\, \frac{k^3 }{1 + {\eta _{1}}(k\xi )^2 + (k\xi )^4}. $$

This integral develops a logarithmic divergence as ${k_{c}}\rightarrow \infty$. Hence, the SSRF is mean-square non-differentiable. $\square$

Appendix B: Proof of Proposition 2

Proof

Let $J_{\nu }(\cdot )$ be the Bessel function of the first kind of order $\nu$, and define

$$ {\mathcal{A}}_{\mu ,\nu }(z)=\int _{0}^{1} dx \, x^{\mu } J_{\nu }(z\,x), $$

(30)

where $z=u_{c}\, h$, $\nu =d/2-1$, and $\mu > -(\nu +1)$. Then, ${\mathcal{A}}_{\mu ,\nu }(z)$ is evaluated using (Gradshteyn and Ryzhik (2007), eq. (6.561.13), p. 676) as follows

$$ {\mathcal{A}}_{\mu ,\nu }(z) =\frac{2^{\mu } \Gamma \left( \frac{\nu +\mu +1}{2}\right) }{z^{\mu +1} \, \Gamma \left( \frac{\nu -\mu +1}{2}\right) } + \frac{(\mu + \nu -1) J_{\nu }(z) S_{\mu -1,\nu -1}(z)}{z^{\mu }} - \frac{J_{\nu -1}(z) S_{\mu ,\nu }(z)}{z^{\mu }}.$$

(31)

Further, we use the normalizing variable transformations $x=k/{k_{c}}$, $h=r/\xi ,$ and $u_{c}={k_{c}}\xi$. In view of the dimensionless variables $x, h, u_c$, the integral (12) becomes

$$ C_{\mathrm{xx}}^{BL}(h;{\varvec{\theta }}) =\frac{ u_c^{1+d/2}\, h^{1-d/2}}{(2\pi )^{d/2} \eta _0\, \xi ^{2d} } \int _{0}^{1} d x \, {x^{d / 2} J_{d/2-1}(x h u_c)} \cdot \left[ 1+{\eta _{1}}(x u_c)^2+ (x u_c)^4 \right]. $$

(32)

In light of (30), (32) and using $z =u_{c} h = {k_{c}}r$ as the dimensionless distance, the function $C_{\mathrm{xx}}^{BL}({\mathbf{r}};{\varvec{\theta }})$ defined by (12) is given by

$$ C_{\mathrm{xx}}^{BL}(z;{\varvec{\theta }}) =\frac{g_{0}({\varvec{\theta }})}{ z^{\nu }} \left[ {\mathcal{A}}_{\nu +1,\nu }(z) + {\eta _{1}}u_{c}^2 {\mathcal{A}}_{\nu +3,\nu }(z) \right. \left. + \, u_{c}^4 {\mathcal{A}}_{\nu +5,\nu }(z)\right] ,$$

(33)

$$ g_{0}({\varvec{\theta }}) = \frac{ {k_{c}}^{d}}{(2\pi )^{d/2} \eta _0\, \xi ^{d} }. $$

(34)

For the three terms ${\mathcal{A}}_{\mu ,\nu }(z)$ $(\mu =\nu +1, \nu +3, \nu +5)$ included in $C_{\mathrm{xx}}^{BL}(z;{\varvec{\theta }})$, the parameters $\mu , \nu$ satisfy the relation

$$ \nu - \mu +1 = - 2\,l, \;\text{ where } \; l=0,1,2. $$

(35)

Equations (16) follow directly from (33) which expresses $C_{\mathrm{xx}}^{BL}(z;{\varvec{\theta }})$ in terms of ${\mathcal{A}}_{\mu ,\nu }(z)$, and from (31) which expresses the integrals ${\mathcal{A}}_{\mu ,\nu }(z)$ in terms of Lommel functions. In view of (35), the Gamma function contributions to ${\mathcal{A}}_{\nu +2l+1,\nu }(z)$ in (31) vanish due to the poles of $\Gamma (n)$ at $n \in \mathbb {Z}_{0,-}$.

1.1 Permissibility

The non-negative definiteness of $C_{\mathrm{xx}}^{BL}(z;{\varvec{\theta }})$ is based on Bochner’s theorem and the fact that, according to (13), $\widetilde{C_{\mathrm{xx}}^{BL}}(k;{\varvec{\theta }}) \ge 0$ for ${\eta _{1}}>-2$.

1.2 Differentiability

The existence of the $n$-th order partial derivatives of the Bessel-Lommel SRF in the mean-square sense requires that all the partial derivatives of order $2n$ of $C_{\mathrm{xx}}^{BL}(z;{\varvec{\theta }})$ exist at $z=0$. This condition is ensured by the convergence of the $2n$-th order spectral moment, i.e., of the integral

$$ \varLambda _{d}^{(2n)} = \int _{0}^{{k_{c}}} dk\, {k^{n+d-1} }{\, \left[ 1 + {\eta _{1}}(k\xi )^2 + (k\xi )^4 \right] } $$

for ${k_{c}}\in {\mathbb{R}}$. $\square$

Appendix C: Proof of Proposition 3

Proof

To find the supremum of $f(k):=k^{2\alpha } \, \widetilde{C_{\mathrm{xx}}}(k;{\varvec{\theta }})$ we consider the extremum condition $\mathrm{d} f(k)/\mathrm{d} k=0$, which admits the following two roots:

$$ \tilde{\kappa }_{1,2} = \sqrt{\frac{ \pm \sqrt{{\eta _{1}}^2 \, (1 - \alpha )^2 - 4 \alpha (\alpha -2)} - {\eta _{1}}\, (1 - \alpha )}{2(2 - \alpha ) \xi ^2}}. $$

For $0 \le \alpha < 1$ only $\tilde{\kappa }_{1} \in {\mathbb{R}}$ and $\sup f(k) = f(\tilde{\kappa }_{1})$.

According to (7) the denominator of (23) becomes

$${\mathcal{S}}_d \, \eta _{0} \, \xi ^{1-2\alpha } \, \int _{0}^{\infty } dx \,\phi _{\alpha }(x) $$

(36)

$$ \text{ where } \quad \phi _{\alpha }(x) =\frac{x^{1+2\alpha }}{ 1 + {\eta _{1}}\, \xi ^2 \, x^2 + \xi ^4 \, x^4 }. $$

(37)

To simplify notation we define

$$ I_{\alpha }(\phi ):=\int _{0}^{\infty } dx \phi _{\alpha }(x). $$

(38)

In order to calculate the integral (36) we use Lebesgue's dominated convergence theorem (Schwartz 2008) expressed as follows:

Theorem 2

Let $\phi _{\alpha }(x)$ be a real-valued function $\forall x \in {\mathbb{R}}$ which is integrable $\forall \alpha \in [0, 1]$. If there is a real-valued function $g_{n}(x)$ such that (i) $\lim _{n \rightarrow \infty } \phi _{\alpha }(x) \, g_{n}(x) = \phi _{\alpha }(x), \, \forall x \in {\mathbb{R}}$ and (ii) $|\phi _{\alpha }(x) \, g_{n}(x)| \le \phi ^{*}(x), \forall x \in {\mathbb{R}}$, where $\phi ^{*}(x)$ is an integrable function, then

$$ I_{\alpha }(\phi )= \int _{0}^{\infty } dx \lim _{n \rightarrow \infty } g_{n}(x) \phi _{\alpha }(x) = \lim _{n \rightarrow \infty } \int _{0}^{\infty } dx g_{n}(x) \phi _{\alpha }(x). $$

We define the following auxiliary function

$$ g_{n}(x) : = \frac{ 2^{\alpha }\, \Gamma (\alpha +1) \, J_{\alpha }(x/n)}{ \left( \frac{x}{n}\right) ^\alpha }. $$

(39)

Condition (i) of Theorem 2 is satisfied because $\lim _{n \rightarrow \infty } g_{n}(x)=1$ based on the infinite series expansion of the Bessel function of the first kind around zero ((Watson 1995, p. 40)).

To prove the condition (ii) we apply the following steps.

• 1. For condition (ii) it suffices to show that $|\phi _{\alpha }(x) \, g_{n}(x)| \le \phi _{\alpha }(x)$, because $\phi _{\alpha }(x)$ is integrable. Given that $\phi _{\alpha }(x) >0$, the latter is equivalent to $| g_{n}(x) | \le 1$.

• 2. We use the integral representation of $J_{\alpha }(z)$ given by (Gradshteyn and Ryzhik (2007), 8.411.4), where $z \in {\mathbb{R}}$:

  $$ J_{\alpha }(z) = \frac{2\left( \frac{z}{2} \right) ^{\alpha }}{\Gamma (\alpha +1/2)\Gamma (1/2)} \, \int _{0}^{\pi /2} d\theta \sin ^{2\alpha }\theta \cos \left( z \cos \theta \right) $$

• 3. Since $|\sin ^{2\alpha }(\theta ) \cos \left( z \cos (\theta )|\right) \le 1$ and $\Gamma (1/2) = \sqrt{\pi }$ it follows from the above that $|J_{\alpha }(z)| \le \frac{\left( \frac{z}{2} \right) ^{\alpha } \sqrt{\pi }}{\Gamma (\alpha +1/2)}.$

• 4. In light of this inequality and (39), proving that $|g_{n}(x)| \le 1$ is equivalent to showing that $\mu _{\alpha }:= \Gamma (\alpha +1) / \Gamma (\alpha +1/2) \le \sqrt{\pi }.$

• 5. Based on the inequality $\mu _{\alpha } < \sqrt{\alpha +1/2}$ (valid for $\alpha > -1/4$) the maximum upper bound of $\mu _{\alpha }$ for $0 \le \alpha \le 1$ is $\sqrt{3/2} < \sqrt{\pi }$. Hence, in light of the previous step $|g_{n}(x)| \le 1$. This concludes the proof of condition (ii).

In light of the above, we can use dominated convergence to calculate $I_{\alpha }(\phi )$ as follows

$$ I_{\alpha }(\phi ) = \lim _{n \rightarrow \infty } (2 n)^{\alpha } \Gamma (\alpha +1) \,\tilde{I}_{\alpha }(\phi ) $$

(40)

where

$$ \tilde{I}_{\alpha }(\phi ) = \int _{0}^{\infty } dx \, \frac{J_{\alpha }(x/n) x^{1+\alpha }}{ 1 + {\eta _{1}}\xi ^2 x^2 + \xi ^4 x^4 }. $$

(41)

The integral $\tilde{I}_{\alpha }(\phi )$ is evaluated by means of the Hankel-Nicholson formula (29) ($\nu =\alpha ,$ $\mu =0$ for ${\eta _{1}}\ne 2, \, \mu =1$ for ${\eta _{1}}= 2$) which leads to

$$ \tilde{I}_{\alpha }(\phi ) = \left\{ \begin{array}{cc} \frac{z_{+}^{\alpha } \, K_{\alpha }(z_{+}/n) - z_{-}^{\alpha } \,K_{\alpha }(z_{-}/n)}{\sqrt{{\eta _{1}}^2-4}} &{} {\eta _{1}}\ne 2\\ \frac{K_{\alpha -1}(1/n)}{2n} = \frac{K_{1-\alpha }(1/n)}{2n} &{} {\eta _{1}}= 2. \end{array} \right.$$

(42)

To evaluate $\lim _{n \rightarrow \infty }\tilde{I}_{\alpha }(\phi )$ for $1> p >0$ we use the series expansion (Schwartz 2008) of the K-Bessel function

$$ K_{p}(x) =\frac{1}{2} \left[ \Gamma (p) \left( \frac{2}{x}\right) ^{p} \left( 1+ O(x^2) \right) + \Gamma (-p) \left( \frac{x}{2}\right) ^{p} \right. \left. \quad \left( 1+ O(x^2) \right) \right].$$

For ${\eta _{1}}=2$, $p = 1 - \alpha$, and $x=1/n$ the dominant contribution at $n \rightarrow \infty$ comes from the $O(x^{-p})$ term of the first series on the right hand side, which gives $K_{1-\alpha }(1/n) \approx \frac{1}{2} \Gamma (1 - \alpha ) (2n)^{1 - \alpha }$.

For ${\eta _{1}}\ne 2$ the $O(x^{-p})$ term of the first series cancels out due to the difference between the two Bessel functions, whereas the $O(x^{2-p})$ terms vanish at the limit $n \rightarrow \infty$. A finite contribution comes from the $O(x^p)$ term of the second series on the right hand side, i.e., $z_{+}^{\alpha } \, K_{\alpha }(z_{+}/n) - z_{-}^{\alpha } \,K_{\alpha }(z_{-}/n) \sim \frac{\Gamma (-\alpha )}{2(2n)^{\alpha }}\left( z_{+}^{2\alpha } - z_{-}^{2\alpha } \right)$.

Thus, based on the above asymptotic analysis of the K-Bessel function, (40), (41), and (42) we obtain the following equation (where $\Delta = \sqrt{{\eta _{1}}^2-4}$):

$$ I_{\alpha }(\phi ) =\left\{ \begin{array}{cc} \frac{\Gamma (1-\alpha )\Gamma (1+\alpha ) \, \left[ \left( {\eta _{1}}+ \Delta \right) ^{\alpha } - \left( {\eta _{1}}- \Delta \right) ^{\alpha }\right] }{2^{\alpha +1}\alpha \Delta } &{} {\eta _{1}}\ne 2\\ \frac{\Gamma (1-\alpha ) \, \Gamma (1+\alpha ) }{2} &{} {\eta _{1}}= 2. \end{array} \right. $$

(43)

Finally, based on (43), the definition (38) and (36), (24a) is proved. $\square$

Appendix D: Proof of Proposition 4

Proof

Based on the spectral density (13) it follows that the denominator in (23) is given by

$$\int _{{\mathbb{R}}^d} d{\mathbf{k}}\,k^{2\alpha } {\widetilde{C_{\rm{x}\rm{x}}}}^{BL}(k;{\varvec{\theta }}) = \frac{{\mathcal{S}}_{d} {k_{c}}^{d+2\alpha }}{\eta _{0} \xi ^d} \times \left( \frac{1}{d+2\alpha } + \frac{{\eta _{1}}{k_{c}}^2 \xi ^2}{d+2\alpha +2} + \frac{ {k_{c}}^4 \xi ^4}{d+2\alpha +4} \right). $$

(44)

Let us define the function $\phi (k) = k^{2\alpha } \,\widetilde{C_{\mathrm{xx}}}^{BL}(k;{\varvec{\theta }})$. The numerator in (23) is then given by $\sup _{{\mathbf{k}}\in {\mathbb{R}}^d} \, \phi (k) = \phi (\kappa ^{*})$ where $\kappa ^{*} = \underset{{k} \in [0, {k_{c}}]}{\arg \max } \phi (k)$.

1.1 Non-negative ${\eta _{1}}$

For ${\eta _{1}}\ge 0$, $\phi (k)$ is a monotonically increasing function of $k$; thus, $\kappa ^{*}= {k_{c}}$ and $\phi (\kappa ^{*}) = \frac{{k_{c}}^{2\alpha }}{\eta _{0} \xi ^d}\, \left( 1 + {\eta _{1}}\, {k_{c}}^2 \, \xi ^2 + {k_{c}}^4 \xi ^4 \right)$. In light of (44), this leads to (25a).

1.2 Negative ${\eta _{1}}$

For ${\eta _{1}}<0$, $\phi (k)$ develops local extrema at the wavenumbers that solve the equation $d\phi (k)/dk=0$, i.e., at the $\kappa _{\pm }$ given by (26c)³.

• 1. Complex $\kappa _{\pm }$

For ${\eta _{1}}^2 < 4 \alpha (\alpha +2) / (\alpha +1)^2$, the $\kappa _{\pm }$ are complex numbers and thus $\phi (k)$ does not develop local extrema for $k\in {\mathbb{R}}$. Hence, $\kappa ^{*} = {k_{c}}$ and $\lambda _{c}^{(\alpha )}$ is given by (25a).

• 2. Real $\kappa _{\pm }$

For $4 > {\eta _{1}}^2 \ge 4 \alpha (\alpha +2) / (\alpha +1)^2$, the $\kappa _{\pm }$ are real numbers, one corresponding to the position of a local minimum and the other to a local maximum.

• 1. If $\alpha >0$, $d\phi (k)/dk\propto 2\alpha \, (k\xi )^{2\alpha -1}$ for $k\ll 1$, and thus $\phi (k)$ increases monotonically. The maximum of $\phi (k)$ thus occurs at $\kappa _{-} < \kappa _{+}$.

• 2. If $\alpha =0$, then $d\phi (k)/dk\propto -2|{\eta _{1}}|\, (k\xi )$ and thus $\phi (k)$ decreases monotonically for $k\ll 1$. In this case also $\phi (k)$ becomes maximum at $\kappa _{-}=0,$ whereas the minimum occurs at $\kappa _{+} = \sqrt{|{\eta _{1}}|/2} / ({k_{c}}\xi )$. Again, we distinguish between two cases depending on the relation between $\kappa _{-}$ and ${k_{c}}.$

  • (a) If $\kappa _{-} > {k_{c}}$, then $\kappa ^{*} = {k_{c}}$ and $\lambda _{c}^{(\alpha )}$ is given by (25a). For $\alpha =0$ it holds that $\kappa _{-}=0$ and thus $\kappa _{-} < {k_{c}}$.

  • (b) If ${k_{c}}> \kappa _{-}$, we further distinguish the following cases

    • i. If k _c  < κ ₊, then κ* = κ ₋ and λ ^(α) _c is given by (25b).

    • ii. If k _c > κ ₊, then κ^*= arg max(φ(κ ₋),φ(k _c)) and λ ^(α) _c is given by (25c).

$\square$