# Bayesian inversion with $\alpha$-stable priors

Jarkko Suuronen  School of Engineering Science, Lappeenranta-Lahti University
of Technology, FI-53850, Finland<EMAIL_ADDRESS>,  Tomas Soto  School
of Engineering Science, Lappeenranta-Lahti University of Technology, FI-53850,
Finland<EMAIL_ADDRESS>,  Neil K. Chada  Department of Actuarial
Mathematics and Statistics, Heriot Watt University, Edinburgh, EH14 4AS, UK
<EMAIL_ADDRESS>and  Lassi Roininen  School of Engineering Science,
Lappeenranta-Lahti University of Technology, FI-53850, Finland
<EMAIL_ADDRESS>
###### Abstract.

We propose using Lévy $\alpha$-stable distributions to construct priors for
Bayesian inverse problems. The construction is based on Markov fields with
stable-distributed increments. Special cases include the Cauchy and Gaussian
distributions, with stability indices $\alpha=1$, and $\alpha=2$,
respectively. Our target is to show that these priors provide a rich class of
priors for modeling rough features. The main technical issue is that the
$\alpha$-stable probability density functions lack closed-form expressions,
and this limits their applicability. For practical purposes, we need to
approximate probability density functions through numerical integration or
series expansions. For Bayesian inversion, the currently available
approximation methods are either too time-consuming or do not function within
the range of stability and radius arguments. To address the issue, we propose
a new hybrid approximation method for symmetric univariate and bivariate
$\alpha$-stable distributions that is both fast to evaluate and accurate
enough from a practical viewpoint. In the numerical implementation of
$\alpha$-stable random field priors, we use the constructed approximation
method. We show how the constructed priors can be used to solve specific
Bayesian inverse problems, such as the deconvolution problem and the inversion
of a function governed by an elliptic partial differential equation. We also
demonstrate hierarchical $\alpha$-stable priors in the one-dimensional
deconvolution problem. For all numerical examples, we use maximum a posteriori
estimation. To that end, we exploit the limited-memory BFGS and its bounded
variant for the estimator.

###### Key words and phrases:

Inverse problems, $\alpha$-stable processes, non-Gaussian priors, error bound,
deconvolution

###### 1991 Mathematics Subject Classification:

94A12, 86A22, 60G35, 62M99.

##  1\. Introduction

Inverse problems is the mathematical theory and practical interpretation of
noise-perturbed indirect observations. Bayesian statistical inversion is the
effort to formulate real-world inverse problems as Bayesian statistical
estimation problems [51, 26]. Bayesian inverse problems can be found in
medical and subsurface imaging, industrial applications, and near-space remote
sensing. The objective, for example in industrial tomography, is to detect
different materials, which may have isotropic, anisotropic, or inhomogeneous
features. This means that we typically aim to reconstruct a hidden substance
from indirect noise-perturbed measurements. Inhomogeneities include, for
example, material interfaces and rough features, and these are the main topics
of this paper.

Inverse problems are often formulated using a noise-perturbed measurement
equation

(1.1)
$\mathbf{y}=\mathcal{G}(u)+\bm{\eta},\quad\bm{\eta}\sim\mathcal{N}(\mathbf{0},\bm{C}),$

where $\mathbf{y}\in^{M}$ are noisy finite-dimensional measurements,
$\mathcal{G}$ is a linear or non-linear mapping from some function space to M,
$u:D\rightarrow$ is the unknown with $D\subseteq^{d},\,d=1,2,3$, and
$\bm{\eta}$ is noise, which we assume to be Gaussian. Our aim is to estimate
$u$ given one realization of $\mathbf{y}$.

Inverse problem methods can be roughly divided into deterministic and
statistical methods. We model $\mathbf{y},u,\bm{\eta}$ as random objects in a
statistical framework [26]. For practical computations, we discretize the
unknown $u$, and denote it by $\mathbf{u}$. Then, within the Bayesian inverse
problem framework, the solution can be represented through probability
distributions via Bayes’ theorem, that is, the posterior distribution

$\pi(\mathbf{u}\mid\mathbf{y})=\frac{\
\pi(\mathbf{y}\mid\mathbf{u})\pi(\mathbf{u})}{\pi(\mathbf{y})}\propto\pi(\mathbf{y}\mid\mathbf{u})\pi(\mathbf{u}),$

where $\pi(\mathbf{y}\mid\mathbf{u})$ is the likelihood, and $\pi(\mathbf{u})$
is the prior distribution of the unknown. We omit the normalization constant
$\pi(\mathbf{y})$, and instead use the unnormalized posterior distribution
from now on.

The choice of the prior $\pi(\mathbf{u})$ is practically the only tuneable
object in inversion. The traditional choices in inverse problems are Gaussian
and total variation priors for smoothing and edge-preserving inversion,
respectively [26, 46, 45, 30]. In this paper, we build upon the research line
starting from the observation that total variation priors do not provide
invariant estimators under mesh refinement [30]. Besov priors on a wavelet
basis were proposed as a solution to this problem [29]. Here, we extend the
study from Cauchy [50] priors to $\alpha$-stable priors, of which the Cauchy
priors are special cases with $\alpha=1$, and Gaussian priors are similarly
special cases with $\alpha=2$. In order to leverage $\alpha$-stable laws for
Bayesian inverse problems, we need approximations of $\alpha$-stable
probability densities evaluated very fast with reasonable precision [44, 43].
Our particular interest is to implement and use discretized $\alpha$-stable
random fields in Bayesian continuous-parameter estimation.

We note that traditionally, stable distributions have been employed in
financial applications like the modeling of asset time series [7]. They have
also been used in biomedical engineering [1], remote sensing [35], network
traffic statistical analysis [22], and digital signal processing [40], to name
a few. We extend the application to inverse problems here.

###  1.1. Contributions

Our main contribution is to implement computationally feasible numerical
approximations of symmetric $\alpha$-stable priors for Bayesian inverse
problems, which requires evaluating the univariate or multivariate probability
density functions of $\alpha$-stable random variables. The symmetric
$\alpha$-stable probability density functions do not have elementary function
expressions, except for the two special cases of Gaussian and Cauchy
distributions. Thus, the evaluation of the probability densities requires the
incorporation of an appropriate approximation method.

Unfortunately, none of the existing $\alpha$-stable density function
approximation methods is alone optimal for our needs. Being based on either
series expansions, integral expressions or Fourier transforms [43, 47, 17, 44,
3, 34, 11, 52], they are either computationally too heavy to evaluate within
Bayesian inversion, are not applicable for arbitrary values of $r$ and
stability indices $\alpha$, or do not provide a consistent seamless
approximation in the sense of the partial derivatives.

For this reason, we introduce a fast hybrid method to approximate the
$\alpha$-stable laws that leverage both bicubic spline interpolations at
precomputed probability density grids and asymptotic series approximations.
The novelty of the hybrid method is multi-part. First, the presented method is
orders of magnitude faster to evaluate than the prevalent methods that are
based on numerical integration. The method approximates the log densities of
the $\alpha$-stable distributions within 100-400 nanoseconds on a typical
Intel Xeon-based workstation. For the second, the method enables the
evaluation of the partial derivatives of the log densities with respect to
both $\alpha$ and $r$ with similar performance, which is crucial for the
estimators in the Bayesian inference. Finally, the method provides a
consistent and discontinuity-free approximation of the log densities on a wide
range of stability indices $\alpha$ and for any argument $r$ as opposed to the
series expansion methods.

We establish error bounds for the method through comprehensive computer-
assisted analysis of the partial derivatives of the $\alpha$-stable
probability density functions [33] and argue that the presented approach is
accurate enough for Bayesian inversion. Furthermore, we demonstrate various
$\alpha$-stable priors on a range of Bayesian inverse problems. To our
knowledge, general $\alpha$-stable priors have not been approximated and
implemented numerically for Bayesian inverse problems previously. The
numerical examples include deconvolution problems in one-dimensional and two-
dimensional grids. Additionally, we illustrate nonlinear Bayesian inversion
governed by an elliptic partial differential equation through $\alpha$-stable
priors. In the numerical examples, we resort to maximum a posteriori (MAP)
estimators:

(1.2)
$\mathbf{u}_{\rm{MAP}}:=\arg\max_{\mathbf{u}}\pi(\mathbf{u}\mid\mathbf{y}).$

We showcase the applicability of the presented methodology in a hierarchical
$\alpha$-stable prior, of which the stability index $\alpha$ or the scale
$\sigma$ are processes of their own. The presented hybrid method proved to be
very useful in such a hierarchical scenario because of the variable stability
indices. In essence, we extend previous works on parametric deep Gaussian
processes [54] to simple two-layer $\alpha$-stable processes. The
$\alpha$-stable priors have the ability to model both discontinuities and
smoothness simultaneously in a function of interest, which is a desired
property in Bayesian continuous-parameter estimation. The results are
promising for further development of $\alpha$-stable random field priors. For
instance, the hierarchical $\alpha$-stable priors may be very effective tools
to model processes that have simultaneously properties of a stationary
Gaussian process and a Cauchy random walk.

###  1.2. Outline

This paper is organized as follows: We provide the necessary background
material for the paper as an introduction to $\alpha$-stable priors in Section
2. This will lead onto Section 3, where we briefly review the existing methods
and our hybrid method for approximating $\alpha$-stable probability density
functions, and provide error bounds related to our method. Numerical
experiments with the $\alpha$-stable priors are provided in Section 4, where
we test our priors on the example problems. A summary of our findings and
future work is provided in Section 5.

##  2\. Models

In this section, we review and discuss the necessary prior forms based on
$\alpha$-stable distributions. We also present some basic properties and then
present the multivariate setting and how they can be defined.

###  2.1. Stable distributions

A random variable $W$ corresponding to a symmetric _stable_ distribution, also
known as _$\alpha$ -stable_ and _Lévy $\alpha$-stable_ distribution, can be
characterized in terms of a _stability index_ $\alpha\in(0,2]$ (sometimes also
called the _tail index_ or the _characteristic exponent_), and a scale
parameter $\sigma>0$, in the sense that its characteristic function is given
by

(2.1) $\mathbb{E}\left[\exp(\mathrm{i}\theta
W)\right]=\exp\left(-(\sigma|\theta|)^{\alpha}\right),\quad\theta\in,$

in which case we write

$W\sim\mathcal{S}_{\alpha}(\sigma).$

The parameter $\alpha$ is called the stability index because if $W_{1}$ and
$W_{2}$ are two independent copies of $W$ and $A$, $B>0$, then

(2.2) $AW_{1}+BW_{2}\stackrel{{\scriptstyle d}}{{=}}CW$

with

$C^{\alpha}=A^{\alpha}+B^{\alpha}.$

Hence, the symmetric $\alpha$-stable distributions are a family of continuous
probability distributions that are infinitely divisible and closed under
convolution. The monograph [47] is the standard reference for stable
distributions, including the wide class of non-symmetric stable distributions,
which we do not consider here.

It is clear from (2.1) that for $\alpha=2$, $W$ is normally distributed (with
zero mean and variance $2\sigma^{2}$). Likewise, for $\alpha=1$, $W$ has a
Cauchy distribution (with zero median and scale parameter being $\sigma$).
Besides these two special cases of $\alpha$-stable laws, there are no known
closed-form expressions based on elementary functions for the density
functions of symmetric stable distributions (the other special cases where
closed-form expressions are known to consist of non-symmetric distributions,
such as the Lévy distribution).

It holds that $\mathbb{E}[|W|^{\alpha}]=\infty$ for $\alpha<2$, implying that
an $\alpha$-stable distribution has infinite variance and that its mean is not
well-defined for $\alpha\leq 1$. It does, however, hold that
$\mathbb{E}[|W|^{\lambda}]<\infty$ for all $\lambda\in(0,\alpha)$.

Multivariate stable distributions can be defined in a similar but more
complicated manner, with a _spectral measure_ $\Lambda$ in place of the scale
parameter $\sigma$; see [47, Chapter 2]. For our purposes, it suffices to
recall the definition of _spherically contoured stable distributions_ : a
random vector $\bi{W}=\mathrel{\mathop{:}}(W_{1},W_{2},\cdots,W_{d})$ on d is
said to have a spherically contoured stable distribution if its characteristic
function is given by

$\mathbb{E}\left[\exp\left(\mathrm{i}\sum_{j=1}^{d}\theta_{j}W_{j}\right)\right]=\exp\left(-(\sigma|\theta|)^{\alpha}\right),\quad\theta\in^{d},$

where $\alpha\in(0,2]$ is again a stability index, $\sigma>0$ is a scale
parameter and $|\cdot|$ stands for the standard $\ell^{2}$-based Euclidean
norm on d. We refer to [43] for a treatment of such distributions.

###  2.2. $\alpha$-stable priors

A stochastic process $(W_{t})_{t\geq 0}$ is a symmetric $\alpha$-stable
process if $\sum_{j=1}^{n}a_{j}W_{t_{j}}$ is a symmetric $\alpha$-stable
random variable for all finite $\\{t_{1},\cdots,t_{n}\\}\subset[0,\infty)$ and
$\\{a_{1},\cdots,a_{n}\\}\subset$. We refer to [47, Chapter 3] for the
existence and construction of a wide variety of such processes. If the
previous definition is satisfied for the multivariate case
$\\{t_{1},\cdots,t_{n}\\}\subset^{K}$, the process is called an
$\alpha$-stable field.

In particular, we aim to apply discretized priors corresponding to a _Lévy
$\alpha$-stable motion_, which we take to mean an $\alpha$-stable process
$(W_{t})_{t\geq 0}$ with some given initial distribution $W_{0}\sim\mu$ and
independent increments that satisfies

$W_{t}-W_{s}\sim\mathcal{S}_{\alpha}(|t-s|^{1/\alpha}),\quad\textrm{for
all~{}}s,\,t\in[0,\infty),\,t\neq s.$

For $\alpha<2$, the Lévy $\alpha$-stable motion generally does not have
continuous sample paths. However, according to [48, Theorem 11.1], there
exists a version of this process with cádlág paths satisfying

(2.3) $\mathbb{P}(W_{t}=W_{t-})=1,\quad\textrm{for all~{}}t>0.$

We more generally refer to [48] for an overview of the analytical properties
of Lévy $\alpha$-stable motions and related processes, including a description
of their infinitesimal generators.

A Lévy $\alpha$-stable motion with initial distribution $\mu$ can be
discretized as follows. For $\triangle\in(0,1)$, define the Markov chain
$(u^{\triangle}_{k})_{k\geq 0}$ by $u^{\triangle}_{0}\sim\mu$, and
$u^{\triangle}_{k+1}-u^{\triangle}_{k}\sim\mathcal{S}_{\alpha}(\triangle^{1/\alpha})$
independently for all $k\geq 0$. Then, by writing $(W^{\triangle}_{t})_{t\geq
0}$ for the appropriately-scaled, piecewise constant cádlág process stemming
from the Markov chain $(u^{\triangle}_{k})_{k\geq 0}$, i.e.

$W^{\triangle}_{t}\mathrel{\mathop{:}}=u^{\triangle}_{\lfloor
t/\triangle\rfloor}.$

Using the basic properties of stable distributions along with (2.3), it is
easy to verify that $\lim_{\triangle\to 0^{+}}W^{\triangle}=W$ in the sense of
finite-dimensional distributions, i.e.

$\lim_{\triangle\to
0^{+}}\mathbb{E}\bigl{[}h\bigl{(}W^{\triangle}_{t_{1}},\cdots,W^{\triangle}_{t_{n}}\bigr{)}\bigr{]}=\mathbb{E}\bigl{[}h\bigl{(}W_{t_{1}},\cdots,W_{t_{n}}\bigr{)}\bigr{]},$

for all finite $\\{t_{1},\cdots,t_{n}\\}\subset[0,\infty)$ and bounded and
continuous functions $h\colon^{n}\to$.

An alternative way to construct such a discretization, localized to a finite
interval, is to partition the interval by $N$ equispaced points and define the
unnormalized density function of
$\mathbf{u}\mathrel{\mathop{:}}=(u_{i})_{i=1}^{N}$ on these points as

(2.4)
$\displaystyle\pi(\mathbf{u})\propto\mu(u_{1})\prod_{i=2}^{N}f(u_{i}-u_{i-1};\alpha,\sigma),$

where $\mu$ is the initial distribution of the above-mentioned process
$(W_{t})_{t\geq 0}$ and $f(\,\cdot\,;\alpha,\sigma)$ stands for the stable
density function with stability index $\alpha$ and appropriately-chosen scale
parameter $\sigma$.

The only two-dimensional $\alpha$-stable field we consider in this paper is a
simple generalization of the quasi-isotropic Cauchy first-order difference
prior [50], defined analogously to (2.2). That is, the probability density
function of an $\alpha$-stable random field $\mathbf{u}$ discretized through
finite differences on a two-dimensional rectangular domain $\Omega\subset^{2}$
is proportional to

(2.5)
$\pi(\mathbf{u})\propto\pi_{\partial\Omega}(\mathbf{u}_{\partial\Omega})\prod_{i,j\notin\partial\Omega}f_{B}(u_{i,j}-u_{i,j-1},u_{i,i}-u_{i-1,j};\alpha;\sigma),$

where $\partial\Omega$ denotes the set of the left and bottom indices on the
grid, and $f_{B}(\cdot,\cdot;\alpha,\sigma)$ the symmetric bivariate
$\alpha$-stable probability density function. The probability density function
$\pi_{\partial\Omega}$ is applied on the grid points at the left and bottom
boundary of the grid to make the resulting distribution of $\mathbf{u}$
proper.

###  2.3. Hierarchical $\alpha$-stable priors

Hierarchical priors are dominantly used within Gaussian priors [18, 54, 5].
With these priors, we can model discontinuities and other features with
varying scale or smoothness at the target function. Unfortunately, the
computational complexity of the canonical Gaussian priors is cubic with
respect to the number of training points, unless a special formulation of the
process is employed, like a stochastic partial differential equation [32]. The
hierarchical priors might require several layers on top of each other to
perform well, while having too many layers may not offer any additional
expression capability [20] but rather overfit in the data.

We intend to build and demonstrate simple, two-layer Markovian hierarchical
$\alpha$-stable priors that could be useful without the computational or
implementation complexity of hierarchical Gaussian processes. We model the
scale or the stability of a discretized $\alpha$-stable process as another
$\alpha$-stable process. This is possible due to the first-order difference
prior’s simple Markovian construction, which effectively allows expressing the
normalization constant of the joint distribution of the discretized process
and its parameters in a closed form. A hierarchical $\alpha$-stable difference
process $\mathbf{u}$ with scale $\sigma=G(c)$ and stability $\alpha=H(s)$
based on other discretized $\alpha$-stable difference processes could be
constructed as follows:

(2.6)
$\displaystyle\pi(\mathbf{u},\mathbf{c},\mathbf{s}\mid\mathbf{y})\propto\pi(\mathbf{y}\mid\mathbf{u})\pi(\mathbf{u}\mid\mathbf{c},\mathbf{s})\pi(\mathbf{c})\pi(\mathbf{s})=$
(2.7)
$\displaystyle\pi(\mathbf{y}\mid\mathbf{u})f\left(u_{1};H(s_{1}),G(c_{1})\right)f(c_{1};\alpha_{c},\sigma_{c})f(s_{1};\alpha_{s},\sigma_{s})\times$
(2.8)
$\displaystyle\prod_{i=2}^{N}f\left(u_{i}-u_{i-1};H(s_{i}),G(c_{i})\right)f(c_{i}-c_{i-1};\alpha_{c},\sigma_{c})f(s_{i}-s_{i-1};\alpha_{s},\sigma_{s}),$

where $H$ and $G$ are nonlinear functions with
$\textrm{Range}(G)\subseteq^{+}$, and $\textrm{Range}(H)\subseteq(0,2]$.
$f(\cdot;\alpha,\sigma)$ denotes the probability density function of a
symmetric univariate $\alpha$-stable random variable with stability $\alpha$,
scale $\sigma$. The conditional distribution
$\pi(\mathbf{u}\mid\mathbf{c},\mathbf{s})$ integrates to a constant,
regardless of $\mathbf{c}$ and $\mathbf{s}$, as do the priors
$\pi(\mathbf{c})$ and $\pi(\mathbf{s})$, because
$\sigma_{s},\sigma_{c},\alpha_{s}$ and $\alpha_{c}$ are fixed. The overall
joint prior distribution is thus proper.

The convergence properties to the continuous limit of the hierarchical
$\alpha$-stable processes in the continuous-time limit are unknown – a sum of
two Lévy $\alpha$-stable random variables with different stability indices
does not obey an $\alpha$-stable distribution. However, continuous-time
processes with a local stability index varying with the state of the process,
commonly called _stable-like processes_ , are well-studied in the literature;
see e.g. Chapter 7 in [27] and Theorem 5.2 in [28].

Further studies are thus needed regarding the matter, but as we demonstrate in
the numerical experiments, the hierarchical $\alpha$-stable process
constructions are promising. Unfortunately, the presented hierarchical priors
cannot be applied to the $\alpha$-stable difference priors when the spatial
dimension is greater than one. That is because the normalization constants of
the priors are intractable due to their construction upon the distributions of
increments between nearest neighbors. However, a Matérn-like stochastic
partial differential equation prior could be optionally employed instead of
the difference priors [50], which would effectively allow incorporating deep
$\alpha$-stable processes thanks to the tractable normalization constants.

##  3\. Approximation of $\alpha$-stable probability density functions

Before presenting our hybrid method of approximating symmetric spherically-
contoured $\alpha$-stable density functions, we perform a brief literature
overview of the existing approximation methods as a motivation for our
contributions. Later, we provide various relative error bounds for the
probability density approximations given by our method. For simplicity, we
denote with $r$ both the argument of the univariate probability density
functions and the Euclidean distance of the arguments of the multivariate
$\alpha$-stable probability density functions to the origin, that is,
$r:=\|\mathbf{r}\|$.

Unless otherwise specified, the approximations are applied for $\sigma=1$.
Recall that for general symmetric $\alpha$-stable laws, the probability
density functions for the other scale parameters are given by
$f(r;\alpha,\sigma)=\frac{1}{\sigma^{d}}f(\frac{r}{\sigma};\alpha,1)$, where
$d$ stands for the dimensionality of the distribution [43].

###  3.1. Prevalent approximation methods for symmetric univariate and
bivariate $\alpha$-stable laws

A canonical method to approximate the $\alpha$-stable probability density
functions is to evaluate the inverse Fourier transform of the characteristic
function. For the symmetric univariate $\alpha$-stable distributions given by
(2.1) with $\sigma=1$, the probability density function can be expressed as
[42]

(3.1)
$f(r;\alpha)=\frac{1}{\pi}\int_{0}^{\infty}\cos(rt)\exp(-t^{\alpha})\mathrm{d}t.$

In theory, numerical integration allows evaluating the density of any
$\alpha$-stable distribution at an arbitrary point $r$ with specified
precision. The integral may be impractical to evaluate for large $r$ and small
$\alpha$ due to the severe oscillations [42], so oscillatory integral
techniques have been proposed to address the issue [3]. Furthermore, the
discrete Fourier transform method [49] exploits the low computational
complexity of the fast Fourier transform, but requires interpolation to
evaluate the density outside the grid [36, 34].

Additionally, there is an alternative integral representation formula [44],
that we call Nolan’s method for short, for the univariate $\alpha$-stable
density function when $\alpha\neq 1$. For simplicity, if we consider only
symmetric $\alpha$-stable distributions and $\alpha\neq 1$, the probability
density function for an $\alpha$-stable random variable with $\mu=0$ and
$\sigma=1$ given by the method is [44]

(3.2)
$f(r;\alpha)=\frac{\alpha|r|^{\frac{1}{\alpha-1}}}{\pi|\alpha-1|}\int_{0}^{\pi/2}Q(t,\alpha)\exp\left(-|r|^{\frac{\alpha}{\alpha-1}}Q(t,\alpha)\right)\mathrm{d}t,$

where $Q(t,\alpha)=\left(\frac{\cos(t)}{\sin(\alpha
t)}\right)^{\frac{\alpha}{\alpha-1}}\frac{\cos(\alpha t-t)}{\cos(t)}$. In
contrast to the inverse Fourier method, the integrand in Nolan’s method is
compactly supported and non-oscillatory. Unfortunately, when
$|\alpha-1|<0.02$, the integrand becomes extremely peaky and narrow, and then
the method is impractical unless arbitrary precision arithmetic is used for
the integration [44]. Even otherwise Nolan’s method can be tricky to evaluate
since the domain of integration should be evaluated in parts near the peak of
the integrand. The peak is located at $t_{p}$, which satisfies the equation
$Q(t_{p},\alpha)|r|^{\frac{\alpha}{\alpha-1}}=1$ [44]. Thus, the univariate
symmetric $\alpha$-stable laws can be accurately evaluated with either Nolan’s
method or the inverse Fourier transform method depending on the values of
$\alpha$ and $r$. The approximation methods based on the Fast Fourier
transform [49, 6, 34, 36] are also worth mentioning. They are simple to
implement and relatively fast to evaluate but must be used in conjunction with
interpolation to approximate the density at a point that is not part of the
FFT grid. It has been reported that the FFT-based approximation is accurate
only for large $\alpha$ [10, 34].

The approximations based on the integral representations of $\alpha$-stable
laws are complemented by series expansions. A well-known series expansion for
the univariate $\alpha$-stable density is of form [8]

(3.3)
$f(r;\alpha)\sim-\frac{1}{\pi}\sum_{k=1}^{\infty}\frac{(-1)^{k}\Gamma(k\alpha+1)\sin(\frac{k\pi\alpha}{2})}{k!}r^{-k\alpha-1},$

which is an asymptotic expansion for $\alpha\in(1,2)$ for $r\to\infty$, and
converges pointwise to the true density for $\alpha\in(0,1)$. There is a
similar series expansion, also outlined in [8], which is an asymptotic series
for $\alpha\in(0,1)$ and a converging series for $\alpha\in(1,2)$ at $r\to
0^{+}$. Furthermore, there are methods that provide a converging series
approximation for the symmetric univariate probability densities for
$\alpha\in(0,2)$ by combining two separate power series expansions [11], or
approximate the inverse Fourier transform of the characteristic function by
domain partitioning and implementing different series expansions within them
[17].

The methodology for approximating spherically contoured multivariate
$\alpha$-stable distributions is similar to the univariate one. There are
several integral expressions for their probability density functions (see e.g.
[43]), such as

(3.4) $\displaystyle f_{M}(r;\alpha_{M})$
$\displaystyle=\frac{2^{1-d/2}}{\Gamma(d/2)}\int_{0}^{\infty}(rs)^{\frac{d}{2}}J_{d/2-1}(rs)\exp{\left(-s^{\alpha}\right)}\mathrm{d}s$

where $J_{\nu}$ is the Bessel function of the first kind. Analogously to the
univariate case, multivariate spherically contoured $\alpha$-stable laws have
an absolutely converging series expansion for $r>0$ and $\alpha\in(0,1)$,
which is an asymptotic expansion for $\alpha\in(1,2)$ and $r\to\infty$ [43]:

(3.5)
$\eqalign{f_{M}(r;\alpha)\sim\frac{-1}{2\pi^{d/2+1}}\sum_{k=1}^{\infty}\frac{(-1)^{k}\Gamma(\frac{k\alpha+2}{2})\Gamma(\frac{k\alpha+d}{2})\sin(\frac{k\alpha\pi}{2})}{k!}\Bigl{(}\frac{r}{2}\Bigr{)}^{-(k\alpha+1)}.}$

Similarly, there is an absolutely converging series expansion for $r>0$ and
$\alpha\in(1,2)$, which is an asymptotic expansion for $\alpha\in(0,1)$ and
$r\to 0^{+}$ [43].

None of the existing approximation methods as such is suitable to be used in
Bayesian inversion. The presented approximation techniques based on numerical
integration are not applicable for the tails of the $\alpha$-stable
distributions because the relative accuracy of the integration cannot be
ensured due to limited floating point precision. They and the advanced series
expansions [17] may be too time-consuming to perform within Bayesian
continuous-parameter estimation since the probability density functions must
be evaluated up to several hundreds of thousands of times even in modest-
dimensional settings. Some of the methods are also applicable only on fixed
$\alpha$ [34], and not on a continuous range of them. Certain approximation
methods fail to approximate the log density of the $\alpha$-stable
distributions with a close to uniform relative accuracy for all $r$ because
they are based on combining two power series expansions [11], which may affect
negatively the estimators. Finally, the typical asymptotic series expansions
are fast to evaluate when the number of used terms is low, but not accurate
enough or even applicable for all $r$, let alone $\alpha$.

###  3.2. Hybrid method for approximating $\alpha$-stable laws

$0$$0.9$00.511.52$\alpha$$29$$30$$r$ Figure 1. Regions of the hybrid
interpolation method for approximating symmetric $\alpha$-stable laws.
Turquoise: the first bicubic interpolation grid. Violet: the second bicubic
interpolation grid. Red: the transition region of the spline and the
asymptotic series. Yellow: the asymptotic series expansion for
$r\rightarrow\infty$. Gray: the implemented approximation method is not
employed.

To address the issues of the existing approximation methods in Bayesian
inversion, our hybrid approximation method is divided into parts, and
different techniques are employed within them. When $r$ is small, we
approximate the $\alpha$-stable laws with two-variable bicubic splines that
are fitted on grids of precomputed $\alpha$-stable log densities with varying
radius $r$ and stability $\alpha$. To our knowledge, a similar rationale has
been presented by [34], but our method works also for $\alpha<1$. Furthermore,
our method approximates the probability densities on a continuous range of the
stability indices $\alpha$ with consistent accuracy and regularity, and for a
continuous range of $\alpha$. The feature is essential for hierarchical
$\alpha$-stable prior constructions. We employ the Julia library
Interpolations.jl to evaluate the bicubic splines. Figure 1 depicts the
overall approximation method. We present the details of the parts of the
hybrid method below.

####  3.2.1. The first bicubic spline approximation

The first bicubic spline grid of precomputed log densities is applied when
$r\in[0,0.9],\,\alpha\in[0.5,1.9]$. This region is illustrated in turquoise in
Figure 1. We divide the domain $[0,1.0]\times[0.5,1.9]$ uniformly to the
intervals of $h_{r}=0.01$ and $h_{\alpha}=5\cdot 10^{-4}$, and evaluate
numerically the densities using the integral methods. For the one-dimensional
$\alpha$-stable density we use the Fourier integral formula from (3.1) for
$|\alpha-1|<0.2$, and otherwise we employ Nolan’s method of (3.2). Instead of
introducing improvised heuristics for the integration and domain partitioning,
we count on the Fourier integral for the aforementioned stability values that
are particularly tricky in Nolan’s method. The numerically evaluated densities
of both methods agree each with a least 12 decimals for $|\alpha-1|\geq 0.2$,
so using only the Fourier approach would suffice.

We use the very same spline grid for the bivariate symmetric $\alpha$-stable
laws, and we employ the integral expression with the Bessel function in (3.4).
In both cases, we set the derivative with respect to $r$ at $r=0$ to zero in
the spline to ensure the method approximates the partial derivative with
respect to $r$ properly at $r=0$. At $r=0.9$, we set the partial derivatives
with respect to $r$ to agree with those given by the second spline, which
guarantees the $C^{1}$ continuity of the spline approximation for the log
density at $r=0.9$, which is the upper limit of the first spline.

####  3.2.2. The second bicubic spline approximation

The second bicubic spline grid is constructed in
$r\in[0,30],\,\alpha\in[0.5,1.9]$ with the intervals of $h_{r}=0.01$ and
$h_{\alpha}=5\cdot 10^{-4}$ for the nodes of the grid. The domain is
illustrated in violet in Figure 1. We apply the same approach as in the first
grid For the one-dimensional $\alpha$-stable laws, we use Nolan’s method for
approximating the log densities in the grid points $|\alpha-1|>0.2$, and
otherwise we use the inverse Fourier method. Likewise, we use the integral
expression incorporating a Bessel function for the bivariate $\alpha$-stable
laws to construct the log density grid.

We apply natural boundary conditions with respect to $r$ and $\alpha$ in the
spline, that is, the second derivatives with respect to them are set to zero
on the boundaries of the spline. However, we limit the usage of the spline
only for $r\in(0.9,29.6),\,\alpha\in[0.5,1.9]$. Again, we want the overall
approximation and its partial derivatives to be continuous. That is why we
intentionally let this second spline grid, which is also the largest of them
all, overlap with the domains of the first and the third spline grids, and
modify their boundary conditions to match the partial derivatives with respect
to $r$ of the second spline.

####  3.2.3. Series expansions for approximating the tails

As we pointed out in the overview, approximating the $\alpha$-stable laws for
$r\geq 0$ is ineffective with numerical integration-based approaches, so we
use the asymptotic series expansions for the tails of the distributions. For
the univariate $\alpha$-stable log densities, we use (3.3), and for the
bivariate tails, we use (3.5). We use three terms in both series expansions,
which we found to provide feasible accurate tail density approximations. The
domain of the series expansions is depicted in yellow in Figure 1.

The lower limit of $r=30$ for the series expansions was selected because it
has been reported to provide a practical rounded bound for the asymptotic tail
expansions of the univariate stable densities up to $\alpha=1.999$. However,
the number of required terms in the series can be several tens or over a
hundred for very high accuracy [52]. Naturally, such a high number of terms
would result in a loss of computational efficiency. We keep the same lower
limit for the tail expansions of the symmetric bivariate $\alpha$-stable
distributions since we have numerically verified the difference of the log
densities given by the tail approximation and the numerical integration-based
approach are of the same magnitude as in the univariate case.

####  3.2.4. The transition region

The switch from the spline grids to the asymptotic tails series is made
seamlessly with the help of the third bicubic spline. The spline is employed
for $r\in[29.6,30],\,\alpha\in[0.5,1.9]$, and we call this region the
transition region of the approximation. The purpose of the transition region
is to ensure the log-density approximation and its partial derivatives are
continuous everywhere in its domain. The domain of the transition region
spline is illustrated in red in Figure 1. We use the same grid node intervals
of $h_{r}=0.01$ and $h_{\alpha}=5\cdot 10^{-4}$ in this domain as we do in the
other two grids, and the methodology is identical for both in the univariate
and bivariate $\alpha$-stable laws. The spline is constructed as follows. Let
us denote $r_{a}=29.6$, $r_{b}=30.0$, and $\triangle=r_{b}-r_{a}$. Let the
log-density approximation given by the asymptotic series expansion from (3.3)
((3.5) in the bivariate case) at the spline grid node $r_{a},\alpha_{j}$ be
$f^{s}_{a,j}$ and its derivative with respect to $r$ be $D^{s}_{a,j}$.
Additionally, let $f^{d}_{b,j}$ stand for the log-density approximation given
by the numerical integration at $r_{b},\alpha_{j}$, and $D^{d}_{b,j}$ denote
its derivative with respect to $r$. We set the value in the transition region
spline grid point value at $i,j$ to follow an auxiliary cubic Hermite
interpolation as follows:

(3.6) $\displaystyle f_{i,j}=$ $\displaystyle\frac{(3\triangle
q_{i}^{2}-2q_{i}^{3})}{\triangle^{3}}f^{s}_{b,j}+\frac{(\triangle^{3}-3\triangle
q_{i}^{2}+2q_{i}^{3})}{\triangle^{3}}f^{d}_{a,j}+\frac{q_{i}^{2}(q_{i}-\triangle)}{\triangle^{2}}D^{d}_{b,j}$
(3.7)
$\displaystyle+\frac{q_{i}(q_{i}-\triangle)^{2}}{\triangle^{2}}D^{s}_{a,j},$

where $q_{i}=r_{i}-r_{a}$.

Equation (3.6) is applied for each stability $\alpha_{j}$ within the
transition grid separately. The Hermite interpolation is only applied during
the construction of the transition spline because the evaluation of the
constructed grid is performed by the bicubic spline library in Julia.
Introducing the Hermite interpolated data points as an additional step at the
transition region helps to avoid abrupt changes in the derivatives of the log
densities near the boundary of the transition region and the tail
approximation.

The derivatives of the spline with respect to $r$ on the boundary $r=29.6$ are
set to follow the values given by the second spline. The same derivative at
$r=30$ is set to follow the derivative given by the asymptotic tail expansion,
respectively.

###  3.3. Error bounds of the hybrid method

We obtain the following error bound for the symmetric $\alpha$-stable log
densities in the domain of the spline grids:

$\sup_{r\leq 30,\;0.5\leq\alpha\leq
1.9}\,|\log\,f_{T}(r;\alpha)-\log\,f_{A}(r;\alpha)|\leq\begin{cases}0.00038&\quad\textrm{(univariate
case);}\\\ 0.22&\quad\textrm{(bivariate case),}\end{cases}$

where $f_{T}\mathrel{\mathop{:}}=f(r;\alpha)$ stands for the (true) density of
the symmetric $\alpha$-stable distribution with $\sigma=1$, and
$f_{A}\mathrel{\mathop{:}}=f_{A}(r;\alpha)$ denotes density given by the
bicubic spline interpolation. The error estimate for the bivariate log density
is orders of magnitudes higher than for the univariate one due to the
significantly larger suprema for the partial derivatives within the domain of
the splines, particularly with small $\alpha$. If the lower bound of $\alpha$
of the approximation domain was increased to 0.7, the bivariate log-density
error estimate would decrease to 0.013, accordingly. The accuracy of the
approximations is enough for our needs in the Bayesian inversion.

For the relative error bounds of the tails, we have the following estimates.
Denoting by $\mathcal{S}_{3}(r;\alpha)$ the sum in (3.3) (resp. (3.5)) with
$3$ in place of $\infty$, and by $f_{T}(r;\alpha)$ the true density, we have

$\sup_{r>30,\;0.5\leq\alpha\leq
1.9}\,|\log\,f_{T}(r;\alpha)-\log\,\mathcal{S}_{3}(r;\alpha)|\leq\begin{cases}0.00097&\quad\textrm{(univariate
case);}\\\ 0.0017&\quad\textrm{(bivariate case).}\end{cases}$

For clarity, we briefly explain the rationale of the procedures behind the
derived bounds. First, we assume the integration error of the probability
densities to be zero within the spline grid points. For the second, we ignore
the transition region spline from the error estimates, because it can be left
out of the hybrid method at the expense of having a less regular overall
approximation. Then we make use of the properties of bicubic splines as
follows [24].

######  Theorem 3.1.

Let the true $\alpha$-stable density and the bicubic spline approximated
density be denoted as above. The error resulting from the bicubic spline
approximation can then be approximated by [12, 24]

(3.8) $\displaystyle\|\log f_{T}-\log f_{A}\|_{\infty}\leq$
$\displaystyle\frac{5}{384}\|(\log
f_{T})^{(4,0)}\|_{\infty}h_{r}^{4}+\frac{81}{64}\|(\log
f_{T})^{(2,2)}\|_{\infty}h_{r}^{2}h_{\alpha}^{2}$ (3.9)
$\displaystyle+\frac{5}{384}\|(\log f_{T})^{(0,4)}\|_{\infty}h_{\alpha}^{4},$

where $h_{r}$ and $h_{\alpha}$ stand for the intervals of the $\log$ density
interpolation grid cells in the directions of the radius and the stability
index, respectively, and the superscripts (i,j) stand for partial derivatives
of the form $\frac{\partial^{i+j}}{\partial r^{i}\partial\alpha^{j}}$.

Due to the lack of any sort of a closed-form expression for $f_{T}$,
estimating the partial derivatives of $\log f_{T}$ appearing in the suprema in
(3.8) involves estimating the partial derivatives of $f_{T}$ from above, and
$f_{T}$ itself from below. Both types of estimates are tricky to do in a
precise manner. When estimating the partial derivatives of $f_{T}$, we will
use several strategies that are variably efficient for different regions of
$(r;\alpha)$, and then use the first-order variants of these estimates in
conjunction with a precomputed grid and the fundamental theorem of calculus
for the lower bounds of $f_{T}$. The full details of these estimates are
presented in the Supplementary Material.

First, we can obtain crude uniform bounds (with respect to $r$) for each
$(f_{T})^{(i,j)}$ by e.g. differentiating (3.1) under the integral sign and
eliminating the resulting oscillatory term simply using the triangle
inequality. We can somewhat refine this pointwise for “moderate” $r>1$ by
using standard oscillatory integral techniques (which basically amount to
partial integration against a sufficiently quickly vanishing function,
resulting in an upper bound of order $r^{-1}$). This is largely the best we
can do for moderate values of $r$, where neither of the asymptotic expansions
(see (3.3) and the subsequent discussion) come close to approximating $f_{T}$
well with only a couple summands – the region of said “moderate” values will,
of course, depend on $\alpha$.

For larger values of $r$, we may exploit the expansion (3.3) pointwise with
e.g. $2$–$3$ summands and an explicit (albeit complicated) expression for the
remainder term, due to Bergström [8]. Some of the integrals we encounter here
are highly intractable in the mathematical sense of the word, but non-
oscillating and well-behaved enough for efficient numerical estimation,
yielding an upper bound that decreases to order roughly comparable to that of
$f_{T}$ for $r\to\infty$, and where the asymptotic constants stay sufficiently
tame for $\alpha\in[0.5,1.9]$. This all holds, mutatis mutandis, for $r\to
0^{+}$ as well.

As a result, we obtain several different kinds of upper bounds for
$|(f_{T})^{(i,j)}(r;\alpha)|$, pointwise with respect to $(r,\alpha)$. With
some minor additional work, we may loosen the estimates very slightly so that
they will be uniform for the spline grid cells $r\in[r_{j},r_{j}+h_{r}]$ and
$\alpha\in[\alpha_{i},\alpha_{i}+h_{\\{}\alpha]$, where $i$ and $j$ are
indexed over the numbers of the spline grid points.With the upper estimates
for the partial derivatives of $f_{T}$, we then estimate $f_{T}$ from below
with reasonable accuracy. Namely, noting that $f_{T}(r;\alpha)$ is for fixed
$\alpha$ always a decreasing function of $r$, we may precompute
$f_{T}(r_{j};\alpha_{i})$ at the nodes of the discretized grid, and thus use
the fundamental theorem of calculus to obtain the lower bound for
$f_{T}(r_{j};\alpha_{i})$ within the spline grid cells.

The discussion above also applies to the bivariate case, with the additional
difficulty of the Bessel function of the first kind $J_{0}$ in the
representation (3.4). In the Supplementary Material, we present analogous
asymptotic expansions with respect to $r\to 0^{+}$ and $r\to\infty$ with
quantitative remainder term estimates. In particular for $r\to\infty$, we
present a modification of Bergström’s [8] complex-analytical treatment, which
allows us to obtain estimates for the remainder term in the bivariate case
which are not immediate from the asymptotic expansions presented in [43].

###  3.4. Practical considerations on the hybrid method

We consider the computational efficiency of the presented hybrid scheme as
good. The evaluation of an $\alpha$-stable log density takes approximately 100
nanoseconds in the domain of the spline grids, and about 400 nanoseconds in
the asymptotic tail expansions on a workstation equipped with Intel Xeon CPU
E5-2698 v4 central processing unit. For instance, evaluation of the univariate
$\alpha$-stable density through Nolan’s method or the inverse Fourier
transform approach with a relative tolerance of $10^{-10}$ takes approximately
10-150 microseconds for $\alpha=1.6\,r\in[0,30]$ using Julia’s QuadGK.jl
library, which applies adaptive Gauss–Kronrod quadrature. However, for
$\alpha=1.6\,r=150$, both integration-based methods take over 2 milliseconds
to evaluate, and the computational efficiency degrades the further the greater
the radius argument $r$ is because the integrals are close to zero, and hence
slower to evaluate with the specified relative tolerance. On a typical
Bayesian inverse problem, the number of separate evaluations of probability
density functions can be tens or hundreds of thousands in a posterior
distribution. Furthermore, an approximation based on combining power series
expansions such as the method of [11] is likely not any faster than the hybrid
method that combines splines with well-known asymptotic tail expansions since
the method requires evaluation of two different series expansions
simultaneously. That is why a hybrid method is essential from a computational
efficiency viewpoint, and we argue that our method certainly meets those
requirements.

In the presented methodology, we do not consider the stability indices of
$\alpha<0.5$ or $\alpha>1.9$. In theory, values of $\alpha$ close to
$2-10^{-6}$ can be effectively approximated with the presented hybrid method,
because the series approximation and the integration-based spline
interpolation agree well at $r=30$ for that high $\alpha$. However, it must be
borne in mind that the closer the stability index is to 2, the worse the
relative error of the approximation will be because of the magnitudes of the
partial derivatives with respect to $\alpha$. As a remedy to that, an extra
spline grid could be incorporated in the approximation for $1.9\leq\alpha\leq
2-10^{-6}$, but the problem in the approach is to smoothly join the splines
through their boundary conditions with respect to $\alpha$, which would
further complicate the overall workflow. On the other hand, low stability
indices are not in our interest, and including them would require increasing
the number of nodes in the precomputed log-density grids to sustain the
accuracy of the approximation. Alternatively, yet another spline grid could be
employed for smaller $\alpha$ with a possibly greater range of $r$, which
would again make the method more complex. We remark that the inverse Fourier
transform-based method copes badly at approximating the univariate
$\alpha$-stable densities with a small stability index, whereas Nolan’s method
still sustains its effectiveness thanks to the compact integrand it adapts.

Finally, we elaborate on the selection of the spline parameters. The first and
the transition region splines ($r\in[0,0.9]$ and $r\in[29.6,30.0]$) are useful
also because the resulting systems of equations of the spline coefficients
involving the boundary conditions of the splines are smaller, and hence easier
to solve than directly incorporating them into the coefficients of the largest
spline grid. That is because the interpolation library Interpolations.jl
applies effectively the symmetries of the spline coefficient equations when
natural boundary conditions are used. The limit of $r=29.6$ for the transition
region splines was selected as a good compromise between the regularity of the
approximation and its relative accuracy. The limit of the first spline
($r=0.9$) can be argued to be adequate for the first spline because the
partial derivative with respect to $r$ is the largest with the small $\alpha$,
and the largest magnitudes are obtained for $\alpha<0.9$. We have not
optimized the limit further because all spline grids have the same node
intervals. However, it would be straightforward to decrease the spacing
$h_{r}$ in the first spline, if the lower stability indices will be of
interest in the future. As a proof of concept, we regard the current grid
intervals of $h_{r}=0.01$ and $h_{\alpha}=5\cdot 10^{-4}$ as a practical
compromise between accuracy and the required memory to store grids. The
current parameters result in a grid of approximately $3000\times 3000$ nodes
in the second interpolation spline. Although all the grids in the study need
only a few hundred megabytes of storage in total, we emphasize that the
asymptotic series expansions are virtually the largest source of the relative
error in the overall approximation despite the proven error bounds. Thus, we
consider decreasing the node intervals impractical for now.

##  4\. Numerical examples

We demonstrate the $\alpha$-stable priors in three numerical experiments. We
employ the priors first in a deconvolution, which is a well-known linear
inverse problem. Moreover, the same priors are used in estimating the
conductivity field of an elliptic partial differential equation in two spatial
dimensions. For the time being, we only use MAP estimators in the
reconstructions because full Bayesian inference with the presented random
field priors requires the usage of MCMC methods, which have been shown to
struggle with such heavy-tailed priors [50]. Since the assessment of the
reconstructions in inverse problems cannot be usually accomplished in a
unified manner, we do not intentionally tabulate any metrics of the
reconstructions, such as $L^{2}$ errors of the reconstructions, in the
manuscript. We demonstrate the MAP estimators of the $\alpha$-stable priors by
varying the stability index $\alpha$ and the scale $\sigma$ in the examples.
We did not optimize their ranges in the experiments. Rather, they were
selected so that their effect on the estimators would be meaningful to
illustrate. The Julia codes of the experiments can be found at
https://github.com/suurj/alpha-stable.

###  4.1. MAP estimation

Evaluation of the MAP estimators (1.2) in Bayesian continuous-parameter
estimation is usually performed with the help of a nonlinear conjugate
gradient algorithm, a quasi-Newton method, a matrix-free truncated Newton
method, or a combination of them [2, 25, 38, 53].

In our numerical examples, the maximization of the log posteriors is done
through the L-BFGS method in the deconvolution experiments. Moreover, we
resort to the bounded L-BFGS algorithm [55] at the inversion of the
conductivity field of a linear elliptic partial differential equation. As the
numerical implementations of the limited-memory BFGS algorithms, we use
Optim.jl [37] for the unconstrained L-BFGS, and a Julia wrapper LBFGSB.jl of
the original Fortran-implementation of L-BGFS-B [55]. Lastly, we want to
emphasize that the presented $\alpha$-stable random field priors often make
the posteriors multimodal [50], and finding global maxima from such
distributions is difficult. Moreover, different optimization algorithms may
converge to different local minima.

###  4.2. One-dimensional deconvolution

In the first numerical experiment, we demonstrate the properties of the first-
order $\alpha$-stable difference priors. We discretize the target function
$u$, which includes both discontinuities and piecewise realizations of a
Gaussian process with Matérn covariance, at 500 equispaced grid points on
$[-1,1]$. We convolve $\mathbf{u}$ with a normalized Gaussian kernel of

(4.1) $\psi(p)=25\exp\left(-50\|p\|^{2}\right)$

though a matrix approximation, which is denoted by $\mathbf{F}$. We construct
the noise-corrupted data $\mathbf{y}$ at 60 equispaced points within the
support of the target function by adding white Gaussian noise to the discrete
convolution as follows:

(4.2)
$\mathbf{y}=\mathbf{F}\mathbf{u}+\bm{\epsilon},\,\bm{\epsilon}\sim\mathcal{N}(\bm{0},0.02^{2}\mathbf{I}).$

Thus, the likelihood $\pi(\mathbf{y}\mid\mathbf{u})$ of the forward model is
Gaussian. We avoid committing an inverse crime by reconstructing $\mathbf{u}$
on a 120-dimensional equispaced grid at $[-1,1]$ and use the corresponding
matrix approximation for the convolution with the kernel (4.1) on the
reconstruction grid to evaluate the forward model and the likelihood. The
ground truth function $u$ and the measurement data $\mathbf{y}$ are depicted
in Figure 2(a).

We employ four different $\alpha$-stable priors in the one-dimensional
setting. Namely, we exploit the hierarchical $\alpha$-stable priors defined in
(2.6), so that

  1. (1)

fix the scale $\sigma$ and stability index $\alpha$ of the prior of the
increments of $\mathbf{u}$,

  2. (2)

consider the stability index of the prior of the increments of $\mathbf{u}$ as
a process that depends on another $\alpha$-process $\mathbf{s}$ and fix the
scale of the increments,

  3. (3)

consider the scale of the prior of the increments of $\mathbf{u}$ as a process
that depends on another $\alpha$-process $\mathbf{c}$ and fix the stability
index of the increments,

  4. (4)

and consider both the stability index and the scale of the prior of the
increments of $\mathbf{u}$ as processes that depend on another
$\alpha$-processes $\mathbf{s}$ and $\mathbf{c}$.

The MAP estimates with the priors are illustrated in Figures 2, 3, 4, and 5.
The selected discretization of the processes implies that the dimensions of
the posterior distributions are 120-, 240-, 240-, and 360-dimensional,
respectively.

The MAP estimates with the non-hierarchical $\alpha$-stable first-order
difference priors in Figure 2 demonstrate the effect of altering the stability
$\alpha$ or scale $\sigma$ of the distribution of the increments in the prior
of $\mathbf{u}$. As a rule of thumb, the smaller the stability $\alpha$ is,
the stronger the prior favors non-Gaussian increments, so they are usually
close to zero. The larger the scale $\sigma$ is, the greater the variability
is allowed within the increments. Stability indices of $1\leq\alpha\leq 2$,
are particularly useful for reconstructing the target function in this example
case. Those priors are able to favor the existence of Gaussian-like parts of
the ground truth function when needed. If the estimation was done using
stationary Gaussian priors, the MAP estimate would be either over-smoothed and
incapable to locate the discontinuity at the boxcar, or it could detect the
discontinuity at the expense of being very sensitive to noise.

Considering the stability of the prior of $\mathbf{u}$ as another first-order
$\alpha$-stable process, turned out to be less successful. We let the scale of
the process $\mathbf{u}$ to follow an $\alpha$-stable process with scale
$\sigma=0.01$ and set its untransformed stability process $\mathbf{s}$ to
follow an $\alpha$-stable process with the parameters tabulated in Figure 3.
To guarantee that $0.51\leq\alpha\leq 1.9$, we apply a transformation (2.6)

(4.3) $\alpha=H(s)=0.51+1.39S(s),$

where $S(x)=\frac{1}{1+\exp(-x)}$. In both Figures 3 and 5, the stability
processes are shown in their transformed values. The stability process seems
to be either close to constant ($\approx 1.25$) or decreasing towards the
right side of the domain in all the tabulated cases. However, there is some
variation in the stability process in the middle of the domain when the
untransformed process has the parameters $\alpha=0.8,\sigma=0.1$. The
phenomenon may suggest that the stability index being a process does not work
well as a prior. When the stability of the untransformed stability process is
$\alpha_{s}=1.4$ and its scale $\sigma_{s}=0.05$ (Equation (2.6)), the
reconstruction of $\mathbf{u}$ is smooth at first, but as the stability
decreases, the reconstruction becomes more discontinuous and non-Gaussian.

In the one-dimensional deconvolution experiment, the best results in terms of
the reconstruction agreement with the ground truth are obtained when the scale
of $\mathbf{u}$ process is considered a process instead of its stability. We
fix the stability of $\mathbf{u}$ to $\alpha=1.9$ and instead model the
untransformed discretized scale process $\mathbf{c}$ with another
$\alpha$-stable process with the parameters shown in Figure 4. The final scale
process (2.6) is given by a transformation

(4.4) $\sigma=G(c)=0.001+0.05S(c).$

The reconstructions where the untransformed scale process has the scale of
$0.05\leq\sigma_{c}\leq 0.1$ and stability between $0.8\leq\alpha_{c}\leq
1.9$, agree well with the ground truth and even with each other.

For the last, setting both the scale and the stability of $\mathbf{u}$ as
$\alpha$-stable processes seem to suffer from the same issue as the stability
process case. We set the scale of the untransformed stability index process to
$\sigma_{s}=0.05$ (Equation (2.6)), and the stability of the untransformed
scale parameter process to $\alpha_{c}=1.9$. Hence, the scale parameters
$\sigma$ in Figure 5 refer to the scale of the untransformed stability process
($\sigma_{s}$), and the the stability indices to the untransformed scale
process ($\alpha_{c}$). We transform the processes with the same sigmoid
functions as in the other two cases, using Equations (4.4) and (4.3). Either
one or both of the parameter processes remain close to constant throughout the
domain, and the MAP estimates for $\mathbf{u}$ are no better than in the
simpler hierarchical $\alpha$-stable priors. Whether the poor reconstructions
are caused by overfitting, poorly selected hyperparameters of the processes
$\bm{s}$ and $\bm{c}$, unidentifiability, or something else, shall be
investigated in further studies.

(a) Ground truth function (red) and noisy convolutions (blue).

Figure 2. Ground truth, measurements, and the reconstructions with fixed
stability index and scale parameter in the one-dimensional deconvolution
experiment with $\alpha$-stable difference prior. Red lines: ground truth.
Black lines: MAP estimate for the function.  Figure 3. Reconstructions when
considering the stability $\alpha$ as a process. Red lines: ground truth.
Black lines: MAP estimate for the function. Blue lines: MAP estimate of the
stability process.  Figure 4. Reconstructions when considering the scale
$\sigma$ as a process. Red lines: ground truth. Black lines: MAP estimate for
the function. Green lines: MAP estimate of the scale process.  Figure 5.
Reconstructions when considering both the stability $\alpha$ and the scale
$\sigma$ as processes. Red lines: ground truth. Black lines: MAP estimate for
the function. Blue lines: MAP estimate of the stability process. Green lines:
MAP estimate of the scale process. The left axes of the subfigures stand for
the stability indices, and the right ones denote the scale.

###  4.3. Two-dimensional deconvolution

We also conduct a deconvolution experiment in two dimensions with similar
settings as in the one-dimensional case. The ground truth function and the
reconstructions are plotted in Figure 6. We aim to reconstruct the blurred
test function with the help of spherically symmetric bivariate $\alpha$-stable
first-order difference priors (2.5). The ground truth function is supported on
$[-1,1]^{2}$. It is evaluated at a uniform equispaced grid of size $333\times
333$, after which is interpolated at $100^{2}$ points that are scattered
according to a low-discrepancy sequence within the domain of the target
function. The dataset $\mathbf{y}$ is then generated with the help of a matrix
approximation for the convolution with a Gaussian kernel of $\psi$:

$\psi(\mathbf{p})=\frac{150}{\pi}\exp\left(-150\|\mathbf{p}\|^{2}\right),$

and noncorrelated Gaussian noise with variance of $0.05^{2}$ is added to the
result, analogously to Equation (4.2). We use the interpolation and the non-
gridded evaluation points to ensure that there are no artifacts in the
reconstructions that could be explained by too-regular lattice discretizations
of $\mathbf{u}$ in the data generation step. An equispaced grid of $256\times
256$ nodes is used in the MAP estimators to prevent committing an inverse
crime, and we employ a matrix approximation for the convolution model also in
the reconstruction step.

The MAP estimates of the reconstructions are consistent with the one-
dimensional deconvolution experiment. Increasing the stability of the
$\alpha$-stable difference priors manifests in more Gaussian-like features in
the MAP estimates. The distribution $\alpha=0.51$ and $\sigma=0.01$ is
probably too spiky and heavy-tailed as the difference prior since the
reconstruction lacks any features resembling the ground truth objects. A
notable feature is the existence of diagonal discontinuities at certain MAP
estimates, like in the case $\alpha=0.8$ and $\sigma=0.1$ in the object that
consists of two overlapping spheres. Although the construction of the
$\alpha$-stable difference prior incorporates bivariate symmetrically
contoured $\alpha$-stable distributions, the prior is likely not fully
isotropic. In fact, even the isotropic and upwind total variation priors are
not perfectly isotropic, and a method has been proposed to alleviate the issue
[16]. Unfortunately, the technique cannot be applied to the presented
$\alpha$-stable priors, so the matter of alleviating the diagonally
anisotropic reconstructions must be considered separately.

(a) Ground truth function.

(b) Deconvolution of the function.

Figure 6. Ground truth function, its deconvolution, and the MAP estimate
reconstructions of the two-dimensional deconvolution problem with isotropic
$\alpha$-stable difference priors with varying scale $\sigma$ and stability
$\alpha$.

###  4.4. Inversion of an elliptic partial differential equation

As the third numerical experiment with the $\alpha$-stable priors, we consider
the nonlinear inverse problem of estimating a conductivity field $k\in
L^{\infty}(\Omega)$, with Lipschitz domain $\Omega\subset\mathbb{R}^{2}$ of an
elliptic partial differential equation (PDE):

(4.5) $\displaystyle-\nabla\cdot(k\nabla u$ $\displaystyle)=g,\quad
x\in\Omega$ (4.6) $\displaystyle u=0,\quad x\in\partial\Omega,$

with prescribed zero Dirichlet boundary conditions, where $u\in
H^{1}_{0}(\Omega)$ denotes the solution of the PDE, and $g\in
L^{\infty}(\Omega)$. The inversion is done using noisy observations
$\mathbf{y}$ of the solution of the PDE as the likelihood for $k$. We
discretize the PDE with the standard finite difference method. The noise is
assumed to be Gaussian with observations taking the form

(4.7)
$y_{i}=u_{i}+\epsilon_{i},\quad\epsilon_{i}\sim\mathcal{N}(0,0.001^{2}).$

Like in the other two experiments, we evaluate only the MAP estimator of the
problem. We use a bounded limited-memory BGFS algorithm (L-BFGS-B) to
calculate the MAP estimate. We employ the constraint $10^{-5}<k<10^{2}$ to
ensure the well-posedness of the elliptic PDE and to keep the condition number
of the matrix of the system of equations small. The gradient of the log
posterior with respect to the discretized $k$ is calculated through a discrete
adjoint method [19, 23]. In other words, we solve the adjoint equation to
obtain the adjoint $\mathbf{q}$ via the equation

(4.8)
$\left(\frac{\partial\mathbf{E}}{\partial\mathbf{u}}\right)^{T}\mathbf{q}=\left(\frac{\partial\pi\left(\mathbf{y}\mid\mathbf{u}(\mathbf{k})\right)}{\partial\mathbf{u}}\right)^{T},$

where $\mathbf{E}$ denotes the system of finite difference equations of the
discretized PDE, and $\pi\left(\mathbf{y}\mid\mathbf{u}(\mathbf{k})\right)$
the Gaussian likelihood function, which merely consists of solving the PDE
with given $\mathbf{k}$ and evaluating the fidelity of the obtained solution
with respect to $\mathbf{y}$. The gradient of the log posterior
$\pi(\mathbf{k}\mid\mathbf{y})$ with respect to the discretized conductivity
field $\mathbf{k}$ is then

(4.9)
$\frac{\partial\pi(\mathbf{k}\mid\mathbf{y})}{\partial\mathbf{k}}=-\mathbf{q}\frac{\partial\mathbf{E}}{\partial\mathbf{k}},$

since the likelihood depends on $\mathbf{k}$ only through $\mathbf{u}$.

We estimate the conductivity field with the same bivariate $\alpha$-stable
difference priors as we do in the two-dimensional deconvolution. We use a
reconstruction grid of $128\times 128$. To simulate the measurements, and to
avoid committing an inverse crime, we calculate the solution of the PDE using
a larger finite difference grid with a size of $223\times 223$ and interpolate
the solution at $25\times 25$ points that are positioned on the reconstruction
grid according to a low-discrepancy sequence, after we add the noise to the
dataset according to Equation (4.7). The source term function $g$ of (4.5),
the solution of the PDE, and the ground truth conductivity, as well as the
reconstructions, are plotted in Figure 7.

The shape of a double-sphere object in the conductivity field is captured the
best with the priors with smaller stability indices while increasing the
stability seems to favor smooth reconstructions. On the other hand, having a
large scale $\sigma$ may make the prior uninformative. Judging by the shape
and distribution of the values within the reconstruction, the prior with
$\alpha=0.8$ and scale $\sigma=0.1$ could be the best of the tested parameter
choices in this case.

(a) Source term function $g$.

(b) Solution of the PDE with the true $k$.

(c) True conductivity $k$.

Figure 7. Functions used in the nonlinear inversion of the conductivity
function $k$ of an elliptic PDE, and the MAP estimates with symmetric
$\alpha$-stable difference priors with varying choices of the scale $\sigma$
and stability $\alpha$.

##  5\. Conclusion

This work was motivated by the desire to implement approximations of the
$\alpha$-stable random field priors for Bayesian inverse problems. Because
both the Cauchy and Gaussian fields are special cases of the $\alpha$-stable
random fields, our objective was to extend the prior selection to general
$\alpha$-stable priors, which could prove useful in reconstructions where both
Gaussian and non-Gaussian features are present. As the $\alpha$-stable density
functions mostly lack the closed-form expressions, we introduced a
computationally feasible hybrid method for approximating the symmetric
univariate and bivariate $\alpha$-stable probability density functions. The
novelty of the presented method in comparison with the existing approximation
methods is its accuracy and, especially, its performance. The method allows
evaluation of the $\alpha$-stable probability log-density functions within a
stability index range of $\alpha\in[0.5,1.9]$ and radius argument range of
$r\in[0,\infty)$. Furthermore, we provided error bounds for the log-density
approximations through careful computer-assisted analysis. The range of the
stability index and the relative accuracy of the hybrid method could be
potentially improved by introducing more spline grids for large $\alpha$ and
small $\alpha$. Another further improvement would be to extend the domain of
the splines to higher $r$ and optimize the spacing of the spline grids for
greater accuracy.

In the numerical experiments, we employed finite-difference approximations of
the $\alpha$-stable first-order random motion priors at one-dimensional and
two-dimensional deconvolution, and we also addressed the estimation of a
function governed by an elliptic partial differential equation with the same
priors. The MAP estimation was implemented through the standard L-BFGS method
and its bounded variant. Our objective was to illustrate how the parameters,
such as the stability and scale of the $\alpha$-stable priors, can affect the
estimation of the unknown functions. The results are promising in the sense
that the presented priors are both computationally viable, manifest in useful
MAP estimates, and are yet novel compared to the existing random field priors
like Gaussian, Cauchy, Besov, and total variation priors. We also introduced
hierarchical $\alpha$-stable priors in our one-dimensional deconvolution
example. Although the results are already encouraging, we believe the full
potential of the hierarchical $\alpha$-stable priors is yet to be found. An
approximation method like ours is needed for hierarchical $\alpha$-stable
priors because the density functions must be evaluated on a continuous range
of stability indices.

As we introduced new $\alpha$-stable priors and provided examples through MAP-
estimates, we consider extending the estimators to full inference as well as
other $\alpha$-stable priors. For future work, we will consider Bayesian
neural networks [39] with $\alpha$-stable weights, which are possibly non-
symmetric in contrast to all the distributions in this study. We believe the
developed approximations to turn out useful in that case due to the recent
studies on Bayesian neural networks with Cauchy and Gaussian weights [31].
Alternatively, $\alpha$-stable random field approximations through the
stochastic partial differential equation approach could be beneficial [9, 32,
50]. Another consideration would be to test these priors on ensemble Kalman
methods [13, 14, 15], which have been used and tested with hierarchical Cauchy
processes.

## Acknowledgments

We thank Simo Särkkä and Sari Lasanen for useful and interesting discussions.
This work has been funded by the Research Council of Finland (project number
353095).

The authors declare no conflict of interest. All the authors participated in
the presented work with the four criteria as required per IOP Publishing, and
have read and accepted the manuscript before submission.

##  Appendix A Error bounds

We establish quantitative pointwise bounds for symmetric stable density
functions $f(r;\alpha)$ and their partial derivatives, based on series
expansions due to Bergström.

In this appendix, the symbol $\mathcal{R}$ stands for a generic remainder term
associated to a partial series expansion of a given density function. In
particular, the symbol $\mathcal{R}$ will generally have a _different_ meaning
from one line to another. Its precise meaning will always be clear from
context.

## Univariate bounds for $r\to 0$

Recall that for $r>0$, the density is given by the inverse Fourier transform
as follows:

(A.1)
$f(r;\alpha)=\frac{1}{\pi}\int_{0}^{\infty}\cos(rt)\mathrm{e}^{-t^{\alpha}}\mathrm{d}t=\frac{1}{\pi}\Re\Bigl{\\{}\int_{0}^{\infty}\mathrm{e}^{-\mathrm{i}rt-t^{\alpha}}\mathrm{d}t\Bigr{\\}}.$

Let us first recall the asymptotic expansion for $r\to 0^{+}$. We may apply
the Taylor series expansion for the cosine function in the first integral in
(A.1), yielding

$\displaystyle f(r;\alpha)$
$\displaystyle=\frac{1}{\pi}\sum_{k=0}^{n}\frac{(-1)^{k}r^{2k}}{(2k)!}\int_{0}^{\infty}t^{2k}\mathrm{e}^{-t^{\alpha}}\mathrm{d}t+\mathcal{R}(r;\alpha)$
$\displaystyle=\frac{1}{\pi\alpha}\sum_{k=0}^{n}\frac{(-1)^{k}\Gamma(\frac{2k+1}{\alpha})}{(2k)!}r^{2k}+\mathcal{R}(r;\alpha),$

for any $n\in\mathbb{N}_{0}$. Since
$\|\cos^{(\ell)}\|_{L^{\infty}(0,\infty)}\leq 1$ for all $\ell\in\mathbb{N}$,
we can apply Taylor’s theorem and arrive at the following estimate for the
remainder term:

$|\mathcal{R}(r;\alpha)|\leq\frac{r^{2n+2}}{\pi(2n+2)!}\int_{0}^{\infty}t^{2n+2}\mathrm{e}^{-t^{\alpha}}\mathrm{d}t=\frac{\Gamma(\frac{2n+3}{\alpha})}{\pi\alpha(2n+2)!}r^{2n+2}.$

Similarly, for $\ell\in\mathbb{N}$ we can differentiate the first integral in
(A.1) with respect to $r$ to get

(A.2) $\frac{\partial^{\ell}}{\partial
r^{\ell}}f(r;\alpha)=\frac{(-1)^{\lceil\frac{\ell}{2}\rceil}}{\pi\alpha}\sum_{k=0}^{n}\frac{(-1)^{k}\Gamma(\frac{2k+1+2\lceil\frac{\ell}{2}\rceil}{\alpha})}{(2k+o(\ell))!}r^{2k+o(\ell)}+\mathcal{R}(r;\alpha),$

where $o(\ell)=1$ if $\ell$ is odd and $o(\ell)=0$ otherwise, and

$|\mathcal{R}(r;\alpha)|\leq\frac{\Gamma(\frac{2n+3+2\lceil\frac{\ell}{2}\rceil}{\alpha})}{\pi\alpha(2n+2+o(\ell))!}r^{2n+2+o(\ell)}.$

For partial derivatives with respect to $\alpha$, we can use the decomposition
(A.2) for $\ell_{1}\in\mathbb{N}_{0}$, $\ell_{2}\in\mathbb{N}$ and
$\ell\mathrel{\mathop{:}}=\ell_{1}+\ell_{2}\in\mathbb{N}$ to get a
decomposition of the form

(A.3) $\frac{\partial^{\ell}}{\partial
r^{\ell_{1}}\partial\alpha^{\ell_{2}}}f(r;\alpha)=\frac{(-1)^{\lceil\frac{\ell_{1}}{2}\rceil}}{\pi}\sum_{k=0}^{n}\frac{(-1)^{k}}{(2k+o(\ell_{1}))!}\frac{\partial^{\ell_{2}}}{\partial\alpha^{\ell_{2}}}\Bigl{[}\frac{\Gamma(\frac{2k+1+2\lceil\frac{\ell_{1}}{2}\rceil}{\alpha})}{\alpha}\Bigr{]}r^{2k+o(\ell_{1})}+\mathcal{R}(r;\alpha).$

Here the partial derivatives in the summands can be computed explicitly in
terms of polygamma functions, and the remainder term can be estimated as

(A.4)
$|\mathcal{R}(r;\alpha)|\leq\frac{r^{2n+2+o(\ell_{1})}}{\pi\alpha^{\ell_{2}+1}(2n+2+o(\ell_{1}))!}\int_{0}^{\infty}|\log(t)|^{\ell_{2}}t^{\frac{2n+3+2\lceil\frac{\ell_{1}}{2}\rceil}{\alpha}-1}|p_{\ell_{2}}(t)|\mathrm{e}^{-t}\mathrm{d}t,$

where $p_{\ell_{2}}$ stands for the polynomial of degree $\ell_{2}$ given by

$\frac{\partial^{\ell_{2}}}{\partial\alpha^{\ell_{2}}}[\mathrm{e}^{-t^{\alpha}}]=\log(t)^{\ell_{2}}p_{\ell_{2}}(t^{\alpha})\mathrm{e}^{-t^{\alpha}},$

(and $p_{0}\equiv 1$). Generally this integral cannot be evaluated exactly,
but for given values of the parameters it can be estimated numerically rather
efficiently, since the integrand is neither oscillating nor overly peaked.

## Univariate bounds for $r\to\infty$

As in [8], the latter integral in (A.1) can be rotated from the positive real
axis to the line $\\{\tau\mathrm{e}^{\mathrm{i}\varphi}\,:\,\tau>0\\}$ for an
arbitrary $\varphi\in(-\frac{\pi}{\max(2\alpha,1)},0)$, resulting in

(A.5)
$f(r;\alpha)=\frac{1}{\pi}\Re\Bigl{\\{}\mathrm{e}^{\mathrm{i}\varphi}\int_{0}^{\infty}\mathrm{e}^{\mathrm{e}^{\mathrm{i}\beta_{2}}r\tau}\mathrm{e}^{\mathrm{e}^{\mathrm{i}\beta_{1}}\tau^{\alpha}}\mathrm{d}\tau\Bigr{\\}},$

where
$\beta_{1}\mathrel{\mathop{:}}=\beta_{1}(\alpha)\mathrel{\mathop{:}}=\pi+\alpha\varphi\in(\frac{\pi}{2},\pi)$
and
$\beta_{2}\mathrel{\mathop{:}}=\frac{3\pi}{2}+\varphi\in(\frac{\pi}{2},\frac{3\pi}{2})$.
Expanding the term
$\mathrm{e}^{\mathrm{e}^{\mathrm{i}\beta_{1}}\tau^{\alpha}}$ and performing
some elementary calculations for the summands yields

$\displaystyle f(r;\alpha)$
$\displaystyle=\sum_{k=1}^{n}\frac{(-1)^{k+1}\Gamma(k\alpha+1)\sin(\frac{k\alpha\pi}{2})}{\pi
k!}r^{-k\alpha-1}$ (A.6)
$\displaystyle\qquad\qquad+\underbrace{\frac{1}{\pi}\Re\Bigl{\\{}\mathrm{e}^{\mathrm{i}(\varphi+(n+1)\beta_{1})}\int_{0}^{\infty}\tau^{(n+1)\alpha}M_{n+1}(\mathrm{e}^{\mathrm{i}\beta_{1}}\tau^{\alpha})\mathrm{e}^{\mathrm{e}^{\mathrm{i}\beta_{2}}r\tau}\mathrm{d}\tau\Bigr{\\}}}_{=\mathrel{\mathop{:}}\mathcal{R}(r;\alpha)},$

for all $n\in\mathbb{N}$, where $M_{n+1}$ is (the analytic continuation of)
the function $z\mapsto(e^{z}-\sum_{k=0}^{n}\frac{z^{k}}{k!})/z^{n+1}$. Writing
$M_{0}(z)=\mathrm{e}^{z}$ for notational convenience, the above expansion also
holds for $n\in\\{-1,0\\}$, with the understanding that the sum is zero in
this case.

The functions $M_{k}$ above satisfy the bound $|M_{k}(z)|\leq\frac{1}{k!}$ for
$z$ with negative real part, as can be easily verified for $k=0$ and
consequently proved inductively using the recursive formula
$\frac{\mathrm{d}}{\mathrm{d}z}[z^{k}M_{k}(z)]=z^{k-1}M_{k-1}(z)$ for $k\geq
1$. Since the term $\mathrm{e}^{i\beta_{1}}\tau^{\alpha}$ in the integral in
(A.6) has a negative real part, the error term can thus be estimated by

$|\mathcal{R}(r;\alpha)|\leq\frac{1}{\pi(n+1)!}\int_{0}^{\infty}\tau^{(n+1)\alpha}\mathrm{e}^{\sin(\varphi)r\tau}\mathrm{d}\tau=\frac{\Gamma((n+1)\alpha+1)}{\pi(n+1)!|\sin(\varphi)|^{(n+1)\alpha+1}}r^{-(n+1)\alpha-1},$

and since the true value of $\mathcal{R}$ does not actually depend on the
auxiliary parameter $\varphi\in(-\frac{\pi}{\max(2\alpha,1)},0)$, the above
estimate can be improved to

(A.7)
$|\mathcal{R}(r;\alpha)|\leq\frac{\Gamma((n+1)\alpha+1)}{\pi(n+1)!\sin(\pi_{\alpha})^{(n+1)\alpha+1}}r^{-(n+1)\alpha-1},$

where $\pi_{\alpha}\mathrel{\mathop{:}}=\frac{\pi}{2\max(\alpha,1)}$.

Differentiating (A.6) termwise with respect to $r$, we have

(A.8) $\frac{\partial^{\ell}}{\partial
r^{\ell}}f(r;\alpha)=\frac{(-1)^{\ell}}{\pi}\sum_{k=1}^{n}\frac{(-1)^{k+1}}{k!}\frac{\Gamma(k\alpha+\ell+1)\sin(\frac{k\alpha\pi}{2})}{r^{k\alpha+\ell+1}}+\mathcal{R}(r;\alpha),$

with

(A.9)
$|\mathcal{R}(r;\alpha)|\leq\frac{\Gamma((n+1)\alpha+\ell+1)}{\pi(n+1)!\sin(\pi_{\alpha})^{(n+1)\alpha+\ell+1}}r^{-(n+1)\alpha-\ell-1}.$

Finally, for pure and mixed partial derivatives with respect to $\alpha$, we
may use the previous expansion as a stepping stone to obtain

(A.10) $\frac{\partial^{\ell}}{\partial
r^{\ell_{1}}\partial\alpha^{\ell_{2}}}f(r;\alpha)=\frac{(-1)^{\ell_{1}}}{\pi}\sum_{k=1}^{n}\frac{(-1)^{k+1}}{k!}\frac{\partial^{\ell_{2}}}{\partial\alpha^{\ell_{2}}}\Bigl{[}\frac{\Gamma(k\alpha+\ell_{1}+1)\sin(\frac{k\alpha\pi}{2})}{r^{k\alpha+\ell_{1}+1}}\Bigr{]}+\mathcal{R}(r;\alpha),$

with

(A.11)
$|\mathcal{R}(r;\alpha)|\leq\frac{1}{\pi}\int_{0}^{\infty}\tau^{\ell_{1}}\big{|}\frac{\partial^{\ell_{2}}}{\partial\alpha^{\ell_{2}}}\bigl{[}\mathrm{e}^{\mathrm{i}(n+1)\beta_{1}}\tau^{(n+1)\alpha}M_{n+1}\bigl{(}\mathrm{e}^{\mathrm{i}\beta_{1}}\tau^{\alpha}\bigr{)}\bigr{]}\big{|}\mathrm{e}^{\sin(\varphi)r\tau}\mathrm{d}\tau.$

In (A.10), the derivatives in the summands can again be computed in terms of
the polygamma functions if necessary. In (A.11), the partial derivatives can
be estimated in terms of the functions $M_{k}$ introduced above by iterating
the recursive formula $M_{k}^{\prime}=M_{k}-kM_{k+1}$. The integral in (A.11)
will be of the order
$\mathcal{O}(r^{-(n+1)\alpha-\ell_{1}-1}\log(r)^{\ell_{2}})$ for large values
of $r$, and one can a posteriori take $\varphi\to-\pi_{\alpha}$ since again
the true value of $\mathcal{R}(r;\alpha)$ does not depend on $\varphi$.

## Bivariate bounds for $r\to 0$

A random two-dimensional vector $\bi{X}$ obeying a symmetric and _spherically
contoured_ bivariate stable distribution with stability index $\alpha\in(0,2)$
can be described in terms of its characteristic function:

$\mathbb{E}[\exp(\mathrm{i}\bi{t}^{\intercal}\bi{X})]=\exp(-|\bi{X}|^{\alpha}),\qquad\bi{t}\in^{2},$

where $|\cdot|$ refers to the standard $\ell^{2}$-based Euclidean norm. We
refer to e.g. [43] for a comprehensive account on such distributions.

The density function $f_{\bi{X}}$ of $\bi{X}$ can be expressed at
$\bi{x}\in^{2}$ by the inverse Fourier transform of the characteristic
function above, resulting in

$f_{\bi{X}}(\bi{x};\alpha)=\frac{1}{2\pi}\int_{0}^{\infty}J_{0}(|\bi{x}|t)t\mathrm{e}^{-t^{\alpha}}\mathrm{d}t;$

see e.g. [21]. Here $J_{\nu}$ with $\nu=0$ stands for the Bessel function of
the first kind, which can for $\nu\in\mathbb{N}_{0}$ and $z\in\mathbb{C}$ be
expressed as

(A.12)
$J_{\nu}(z)=\sum_{k=0}^{\infty}\frac{(-1)^{k}}{k!(k+\nu)!}\Bigl{(}\frac{z}{2}\Bigr{)}^{2k+\nu}.$

In this section we are interested in estimating partial derivatives of the
density function $f_{\bi{X}}$ in terms of $|\bi{x}|$ and $\alpha$. For this
purpose, we write

(A.13)
$f(r;\alpha)\mathrel{\mathop{:}}=\frac{1}{2\pi}\int_{0}^{\infty}J_{0}(rt)t\mathrm{e}^{-t^{\alpha}}\mathrm{d}t,\quad
r\geq 0,$

for the radial density function of $\bi{X}$.

By well-known integral representation formulas [4, Equation (4.9.11)], we have
the uniform bounds $\|J_{\nu}\|_{L^{\infty}()}\leq 1$ on the real line for all
$\nu\in\mathbb{N}_{0}$, and by standard recurrence relations [4, Equation
(4.6.6)] we thus have $\|J^{(\ell)}_{\nu}\|_{L^{\infty}()}\leq 1$ for all
$\nu\in\mathbb{N}_{0}$ and derivatives $J^{(\ell)}_{\nu}$ of $J_{\nu}$. Hence
we may expand the function $J_{0}$ in (A.13) to obtain

$f(r;\alpha)=\frac{1}{2\pi\alpha}\sum_{k=0}^{n}\frac{(-1)^{k}\Gamma(\frac{2k+2}{\alpha})}{(2^{k}k!)^{2}}r^{2k}+\mathcal{R}(r;\alpha),$

with

$|\mathcal{R}(r;\alpha)|\leq\frac{\Gamma(\frac{2n+4}{\alpha})}{2\pi\alpha(2n+2)!}r^{2n+2}.$

Similarly, differentiating (A.13) with respect to $r$, and (A.12) for $\nu=0$
with respect to $z$, yields

$\displaystyle\qquad\frac{\partial^{\ell}}{\partial r^{\ell}}f(r;\alpha)$
$\displaystyle=\frac{(-1)^{\lceil\frac{\ell}{2}\rceil}}{2\pi\alpha}\sum_{k=0}^{n}\frac{(-1)^{k}\bigl{\\{}\prod_{i=1}^{\ell}(2k+o(\ell)+i)\bigr{\\}}\Gamma(\frac{2k+2+2\lceil\frac{\ell}{2}\rceil}{\alpha})}{(2^{k+\lceil\frac{\ell}{2}\rceil}(k+\lceil\frac{\ell}{2}\rceil)!)^{2}}r^{2k+o(\ell)}+\mathcal{R}(r;\alpha),$

for all $\ell\in\mathbb{N}$, where

$|\mathcal{R}(r;\alpha)|\leq\frac{\Gamma(\frac{2n+4+2\lceil\frac{\ell}{2}\rceil}{\alpha})}{2\pi\alpha(2n+2+o(\ell))!}r^{2n+2+o(\ell)}.$

Subsequently

$\displaystyle\frac{\partial^{\ell}}{\partial
r^{\ell_{1}}\partial\alpha^{\ell_{2}}}f(r;\alpha)$
$\displaystyle=\frac{(-1)^{\lceil\frac{\ell_{1}}{2}\rceil}}{2\pi}\sum_{k=0}^{n}\frac{(-1)^{k}\bigl{\\{}\prod_{i=1}^{\ell_{1}}(2k+o(\ell_{1})+i)\bigr{\\}}}{(2^{k+\lceil\frac{\ell_{1}}{2}\rceil}(k+\lceil\frac{\ell_{1}}{2}\rceil)!)^{2}}\frac{\partial^{\ell_{2}}}{\partial\alpha^{\ell_{2}}}\Bigl{[}\frac{\Gamma(\frac{2k+2+2\lceil\frac{\ell_{1}}{2}\rceil}{\alpha})}{\alpha}\Bigr{]}r^{2k+o(\ell_{1})}$
$\displaystyle\qquad\qquad+\mathcal{R}(r;\alpha),$

for all $\ell_{1}\in\mathbb{N}_{0}$ and $\ell_{2}\in\mathbb{N}$, where

$|\mathcal{R}(r;\alpha)|\leq\frac{r^{2n+2+o(\ell_{1})}}{2\pi\alpha^{\ell_{2}+1}(2n+2+o(\ell_{1}))!}\int_{0}^{\infty}|\log(t)|^{\ell_{2}}t^{\frac{2n+4+2\lceil\frac{\ell_{1}}{2}\rceil}{\alpha}-1}|p_{\ell_{2}}(t)|\mathrm{e}^{-t}\mathrm{d}t.$

Here $p_{\ell_{2}}$ stands for the polynomials introduced in (A.4) above.
Again, the integrals appearing on the right-hand side well-behaved enough for
numerical approximation; see the discussion following (A.4).

## Bivariate bounds for $r\to\infty$

Nolan [43] records an asymptotic expansion for the density function of the
_amplitude distribution_ of $\bi{X}$ as $r\to\infty$, which for the radial
density function (A.13) translates as

(A.14)
$f(r;\alpha)=\frac{1}{\pi^{2}}\sum_{k=1}^{n}\frac{(-1)^{k+1}2^{k\alpha}\Gamma(\frac{k\alpha+2}{2})^{2}\sin(\frac{k\alpha\pi}{2})}{k!}r^{-k\alpha-2}+\mathcal{O}\bigl{(}r^{-(n+1)\alpha-2}\bigr{)},$

for $\alpha\in(0,2)$, $r>0$ and $n\in\mathbb{N}$. This can be obtained by
expressing $\bi{X}$ as a sub-Gaussian vector with respect to a certain totally
skewed univariate stable distribution, which admits a similar asymptotic
series expansion as described in [8]. Below we will obtain (A.14) in an
alternative way that allows us to quantitatively control the error term.

We may rewrite (A.13) as

(A.15)
$f(r;\alpha)=\frac{1}{2\pi}\Re\Bigl{\\{}\int_{0}^{\infty}H_{0}(rt)t\mathrm{e}^{-t^{\alpha}}\mathrm{d}t\Bigr{\\}},$

where $H_{0}\colon\mathbb{C}\setminus(-\infty,0]\to\mathbb{C}$ stands for the
so-called Hankel function of the first kind [4, Section 4.7], defined in terms
of the Bessel functions of the first kind ($J_{0}$) and second kind ($Y_{0}$)
by

$H_{0}(z)=J_{0}(z)+\mathrm{i}Y_{0}(z).$

The functions $H_{\nu}$, $\nu\in\mathbb{N}_{0}$, can be defined similarly in
terms of $J_{\nu}$ and $Y_{\nu}$. By well-known recurrection relations and
connection formulas [41, Sections 10.6 and 10.4], each derivative
$H_{\nu}^{(\ell)}$ can again be expressed as a linear combination of functions
$H_{\nu^{\prime}}$ with $(\nu-\ell)_{+}\leq\nu^{\prime}\leq\nu+\ell$.

Following the contour integration procedure described in [8], we can then
rotate the integral in (A.15) from the positive real axis to the line
$\\{\mathrm{e}^{i\varphi}\tau\,:\,\tau>0\\}$ for an arbitrary
$\varphi\in(0,\frac{\pi}{\max(2\alpha),1})$ to get

(A.16)
$f(r;\alpha)=\frac{1}{2\pi}\Re\Bigl{\\{}\mathrm{e}^{2\mathrm{i}\varphi}\int_{0}^{\infty}H_{0}\bigl{(}\mathrm{e}^{\mathrm{i}\varphi}r\tau\bigr{)}\tau\mathrm{e}^{\mathrm{e}^{\mathrm{i}\beta_{1}}\tau^{\alpha}}\mathrm{d}\tau\Bigr{\\}},$

where
$\beta_{1}\mathrel{\mathop{:}}=\beta_{1}(\alpha)\mathrel{\mathop{:}}=\pi+\alpha\varphi$
is as in (A.5). More precisely, this contour integration and the associated
limiting procedure is permitted because

(A.17)
$|H_{\nu}(z)|\leq\begin{cases}c_{1,\nu}\bigl{(}|z|^{-\nu}+|\log(z)|\bigr{)}&\quad\text{for}\quad
0<|z|\leq 1\\\
c_{2,\nu}|z|^{-\frac{1}{2}}\mathrm{e}^{-\Im(z)}&\quad\text{for}\quad|z|\geq
1\end{cases}$

for all $\nu\in\mathbb{N}_{0}$, with certain constants $c_{i,\nu}$ that we
will not elaborate on; we refer to [4, Sections 4.5 and 4.8] for this and more
comprehensive asymptotic expansions for Hankel functions.

Because of (A.17), we calso expand the term
$\mathrm{e}^{\mathrm{e}^{\mathrm{i}\beta_{1}}\tau^{\alpha}}$ in (A.16) and
integrate termwise to get

$f(r;\alpha)=\frac{1}{2\pi}\sum_{k=0}^{n}\frac{1}{k!}\Re\Bigl{\\{}\mathrm{e}^{\mathrm{i}(2\varphi+k\beta_{1})}\int_{0}^{\infty}H_{0}\bigl{(}\mathrm{e}^{\mathrm{i}\varphi}\tau\bigr{)}\tau^{k\alpha+1}\mathrm{d}\tau\Bigr{\\}}r^{-k\alpha-2}+\mathcal{R}(r;\alpha).$

If we can show that the remainder term above is of the order
$\mathcal{O}(r^{-(n+1)\alpha-2})$ as $r\to\infty$, it follows automatically
from the uniqueness property of asymptotic expansions that the principal term
coincides with the one in (A.14). To this end we recall the functions $M_{k}$
from (A.6) and write

$\mathcal{R}(r;\alpha)=\frac{1}{2\pi}\Re\Bigl{\\{}\mathrm{e}^{\mathrm{i}(2\varphi+(n+1)\beta_{1})}\int_{0}^{\infty}H_{0}\bigl{(}\mathrm{e}^{\mathrm{i}\varphi}r\tau\bigr{)}\tau^{(n+1)\alpha+1}M_{n+1}\bigl{(}\mathrm{e}^{\mathrm{i}\beta_{1}}\tau^{\alpha}\bigr{)}\mathrm{d}\tau\Bigr{\\}},$

and since the term $\mathrm{e}^{\mathrm{i}\beta_{1}}\tau^{\alpha}$ above has a
negative real part by assumption, we may estimate

$|\mathcal{R}(r;\alpha)|\leq\frac{1}{2\pi(n+1)!}\Bigl{(}\int_{0}^{\infty}|H_{0}(e^{\mathrm{i}\varphi}\tau)|\tau^{(n+1)\alpha+1}\mathrm{d}\tau\Bigr{)}r^{-(n+1)\alpha-2}$

which is of the desired form. In light of (A.17), we may further take
$\varphi\to\frac{\pi}{2\max(\alpha,1)}=\pi_{\alpha}$ so that

(A.18)
$|\mathcal{R}(r;\alpha)|\leq\frac{1}{2\pi(n+1)!}\Bigl{(}\int_{0}^{\infty}|H_{0}(e^{\mathrm{i}\pi_{\alpha}}\tau)|\tau^{(n+1)\alpha+1}\mathrm{d}\tau\Bigr{)}r^{-(n+1)\alpha-2}.$

Similarly,

$\displaystyle\frac{\partial^{\ell}}{\partial r^{\ell}}f(r;\alpha)$
$\displaystyle=\frac{(-1)^{\ell}}{\pi^{2}}\sum_{k=1}^{n}\frac{(-1)^{k+1}2^{k\alpha}\bigl{\\{}\prod_{i=2}^{\ell+1}(k\alpha+i)\bigr{\\}}\Gamma(\frac{k\alpha+2}{2})^{2}\sin(\frac{k\alpha\pi}{2})}{k!}r^{-k\alpha-\ell-2}$
$\displaystyle\qquad\qquad+\mathcal{R}(r;\alpha)$

with

$|\mathcal{R}(r;\alpha)|\leq\frac{1}{2\pi(n+1)!}\Bigl{(}\int_{0}^{\infty}|H_{0}^{(\ell)}(e^{\mathrm{i}\pi_{\alpha}}\tau)|\tau^{(n+1)\alpha+\ell+1}\mathrm{d}\tau\Bigr{)}r^{-(n+1)\alpha-\ell-2},$

and further

$\displaystyle\frac{\partial^{\ell}}{\partial
r^{\ell_{1}}\partial\alpha^{\ell_{2}}}f(r;\alpha)=$
$\displaystyle\frac{(-1)^{\ell_{1}}}{\pi^{2}}\sum_{k=1}^{n}\frac{(-1)^{k+1}}{k!}\frac{\partial^{\ell_{2}}}{\partial\alpha^{\ell_{2}}}\Bigl{[}\frac{2^{k\alpha}\bigl{\\{}\prod_{i=2}^{\ell+1}(k\alpha+i)\bigr{\\}}\Gamma(\frac{k\alpha+2}{2})^{2}\sin(\frac{k\alpha\pi}{2})}{r^{k\alpha+\ell_{1}+2}}\Bigr{]}+\mathcal{R}(r;\alpha),$

with

(A.19)
$|\mathcal{R}(r;\alpha)|\leq\frac{1}{2\pi}\int_{0}^{\infty}\big{|}H_{0}^{(\ell_{1})}\bigl{(}\mathrm{e}^{\mathrm{i}\varphi}r\tau\bigr{)}\big{|}\tau^{\ell_{1}+1}\Big{|}\frac{\partial^{\ell_{2}}}{\partial\alpha^{\ell_{2}}}\bigl{[}\mathrm{e}^{\mathrm{i}(n+1)\beta_{1}}\tau^{(n+1)\alpha}M_{n+1}\bigl{(}\mathrm{e}^{\mathrm{i}\beta_{1}}\tau^{\alpha}\bigr{)}\bigr{]}\Big{|}\mathrm{d}\tau.$

The latter integrand can again be estimated in terms of the functions $M_{k}$;
see (A.11) and the relevant discussion. One will then end up with integrals
that are analytically untractable, but not outside the reach of numerical
estimation.

## Uniform and oscillatory bounds for $0\ll r\ll\infty$

Here we present rather crude uniform bounds, and some refinements based on
oscillatory integral techniques, of
$|\frac{\partial^{\ell_{1}+\ell_{2}}}{\partial
r^{\ell_{1}}\partial\alpha^{\ell_{2}}}f(r;\alpha)|$ for “moderate” values of
$r$, where none of the asymptotic expansions discussed earlier come close to
approximating the true function with e.g. $2$–$3$ summands. The precise range
of these “moderate” values of $r$ of course depends on $\alpha$ and the
$\ell_{i}$’s. We present these estimates only in the univariate case, since
the bivariate case can be handled with minor adjustments.

First, from (A.1), we have

(A.20) $\frac{\partial^{\ell_{1}+\ell_{2}}}{\partial
r^{\ell_{1}}\partial\alpha^{\ell_{2}}}f(r;\alpha)=\frac{1}{\pi}\int_{0}^{\infty}\cos^{(\ell_{1})}(rt)\log(t)^{\ell_{2}}t^{\ell_{1}}p_{\ell_{2}}(t^{\alpha})\mathrm{e}^{-t^{\alpha}}\mathrm{d}t,$

where $\cos^{(\ell_{1})}$ stands for the $\ell_{1}$’th derivative of the
cosine function, and $p_{\ell_{2}}$ stands for the polynomial introduced in
(A.4). Thus, a simple application of the triangle inequality yields

$\Bigl{|}\frac{\partial^{\ell_{1}+\ell_{2}}}{\partial
r^{\ell_{1}}\partial\alpha^{\ell_{2}}}f(r;\alpha)\Bigr{|}\leq\frac{1}{\pi\alpha^{\ell_{2}+1}}\int_{0}^{\infty}|\log(t)|^{\ell_{2}}t^{\frac{\ell_{1}+1}{\alpha}-1}|p_{\ell_{2}}(t)|\mathrm{e}^{-t}\mathrm{d}t.$

The latter integral is again usually intractable (although it can be expressed
in terms of the standard gamma function if $\ell_{2}=0$), but numerically
estimable.

In case $\ell_{2}=0$ and $\ell_{1}\mathrel{\mathop{:}}=\ell>0$, we can
slightly refine (A.20) using the oscillatory integral technique of partial
integration against a function that decays sufficiently fast as $t\to 0^{+}$
and $t\to\infty$:

(A.21) $\Bigl{|}\frac{\partial^{\ell}}{\partial
r^{\ell}}f(r;\alpha)\Bigr{|}=\frac{1}{\pi
r}\Bigl{|}\int_{0}^{\infty}\cos^{(\ell-1)}(rt)\frac{\mathrm{d}}{\mathrm{d}t}\bigl{[}t^{\ell}\mathrm{e}^{-t^{\alpha}}\bigr{]}\mathrm{d}t\Bigr{|}\leq\frac{1}{rt}\int_{0}^{\infty}\bigl{|}\frac{\mathrm{d}}{\mathrm{d}t}\bigl{[}t^{\ell}\mathrm{e}^{-t^{\alpha}}\bigr{]}\bigr{|}\mathrm{d}t.$

The derivative inside the last integral can be easily computed, and one
arrives at an upper bound for the integral that can be expressed in terms of
the gamma function.

Similarly, when $\ell_{1}=0$ and and $\ell_{2}\mathrel{\mathop{:}}=\ell>0$, we
have

$\Bigl{|}\frac{\partial^{\ell}}{\partial\alpha^{\ell}}f(r;\alpha)\Bigr{|}=\frac{1}{\pi
r}\Bigl{|}\int_{0}^{\infty}\sin(rt)\frac{\mathrm{d}^{\ell+1}}{\mathrm{d}t\,\mathrm{d}\alpha^{\ell}}\bigl{[}\mathrm{e}^{-t^{\alpha}}\bigr{]}\mathrm{d}t\Bigr{|}\leq\frac{1}{\pi
r}\int_{0}^{\infty}\bigl{|}\frac{\mathrm{d}^{\ell+1}}{\mathrm{d}t\,\mathrm{d}\alpha^{\ell}}\bigl{[}\mathrm{e}^{-t^{\alpha}}\bigr{]}\bigr{|}\mathrm{d}t.$

Again, the derivatives inside the last integral can be computed, and after a
suitable change of variables we arrive at an integral directly directly
proportional to $\alpha^{-\ell}$, where the coefficient of $\alpha^{-\ell}$
can be calculated numerically.

This approach could also in principle be used for mixed derivatives of
$f(r;\alpha)$, and in (A.21) the partial integration trick could be applied
$\ell$ times instead of once. In both cases however the computations become
very unwieldy, and thus we discard them.

## Grid approximations

Recall from the main paper (Section 3.1) that one of our primary goals is to
estimate

(A.22) $\sup_{0<r<30,\;0.5\leq\alpha\leq
1.9}\Big{|}\frac{\partial^{\ell_{1}+\ell_{2}}}{\partial
r^{\ell_{1}}\partial\alpha^{\ell_{2}}}\log f(r;\alpha)\Big{|}$

for $(\ell_{1},\ell_{2})\in\\{(4,0),(2,2),(0,4)\\}$. Recalling that e.g.

$\displaystyle\bigl{(}\log
f(r;\alpha)\bigr{)}^{(4,0)}=-6\frac{f^{(1,0)}(r;\alpha)^{4}}{f(r;\alpha)^{4}}+12\frac{f^{(2,0)}(r;\alpha)f^{(1,0)}(r;\alpha)^{2}}{f(r;\alpha)^{3}}$
$\displaystyle\qquad\qquad-4\frac{f^{(3,0)}(r;\alpha)f^{(1,0)}(r;\alpha)}{f(r;\alpha)^{2}}-3\frac{f^{(2,0)}(r;\alpha)^{2}}{f(r;\alpha)^{2}}+\frac{f^{(4,0)}(r;\alpha)}{f(r;\alpha)},$

and similarly for the other relevant partial derivatives, we may estimate the
suprema in (A.22) by applying the triangle inequality on sums of the above
sort, and by estimating

  1. (1)

the absolute values of the partial derivatives
$f^{(\ell_{1},\ell_{2})}(r;\alpha)$ from above, relying on the pointwise
estimates established in the previous sections, and

  2. (2)

$f(r;\alpha)$ itself from below.

In this section we will give an overview of a numerical procedure to estimate
these functions in a manner that is applicable for estimating (A.22).

We start with a fixed grid of the $(r,\alpha)$-space, with density given by a
parameter $\triangle>0$ such that

$\frac{1.9-0.5}{\triangle}=\mathrel{\mathop{:}}i^{*},\quad\text{and}\quad\frac{30}{\triangle}=\mathrel{\mathop{:}}j^{*},$

are positive integers. In our numerical simulations we have used
$\triangle\mathrel{\mathop{:}}=\frac{2}{10^{3}}$ in order to conserve
computational resources, but other choices are possible as well.

Write then

$Q_{i,j}\mathrel{\mathop{:}}=\bigl{[}(j-1)\triangle,j\triangle\bigr{]}\times\bigl{[}0.5+(i-1)\triangle,0.5+i\triangle\bigr{]}\quad\text{for
}i\in 1{:}i^{*}\text{ and }j\in 1{:}j^{*},$

so that

$\bigcup_{i,j}Q_{i,j}=[0,30]\times[0.5,1.9],$

and the interiors of the $Q_{i,j}$’s are pairwise disjoint. The idea is to
numerically estimate $f(r;\alpha)$ and its partial derivatives uniformly in
these squares, of which there is a finite amount, and the uniform estimates
will not be too far off from the pointwise estimates we obtained in the
earlier sections if the parameter $\triangle$ is sufficiently small.

To this direction, let us first discuss the estimates of the absolute values
of the partial derivatives $f^{(\ell_{1},\ell_{2})}(r;\alpha)$ for
$(r,\alpha)\in Q_{i,j}$. Many of the pointwise estimates we have discussed
above – for example (A.2), (A.3), (A.4), (A.6), (A.7), (A.8),(A.9), (A.10),
(A.11) and the bivariate versions of these estimates – consist of terms that
are monotonous with respect to both $r$ and $\alpha$ either directly, or after
some simple additional upper estimates, such as applying the triangle
inequality in the sums like the ones in (A.2) and (A.8), and considering the
cases $r<1$ and $r\geq 1$ separately for terms of the form
$r^{k\alpha+\ell+1}$.

Some additional care has to be taken for some of the more complicated terms
appearing in the series expansions for the partial derivatives involving the
variable $\alpha$, as well as the respective remainder terms. We explain this
with examples pertaining to the univariate case, but everything here applies
to the bivariate case as well with obvious modifications.

First, the derivatives with respect to $\alpha$ appearing in the summands in
(A.3) and (A.10) can with some effort be computed for $\ell_{2}\in\\{1,2\\}$,
resulting terms involving the gamma function itself and the so-called
polygamma functions of orders $0$ and $1$, which we denote here by
$\psi(0,\cdot)$ and $\psi(1,\cdot)$ respectively. The function $\psi(1,\cdot)$
is strictly positive and decreasing on the entire positive real axis (see e.g.
[4, Theorem 1.2.5]), which implies the functions $|\psi(0,t)|$ and
$|\Gamma(t)|$ are decreasing for $0<t\leq t_{0}$ and increasing for $t>t_{0}$,
where $t_{0}\approx 1.46$ is the positive zero of $\psi(0,\cdot)$, which has
to be taken into account when estimating the derivatives involving $\Gamma$.
For $\ell_{2}>2$, we use (A.3) and (A.10) only with $n=-1$ and $n=0$
respectively, avoiding the need to to estimate these derivatives at all.

Secondly, the integrals appearing in (A.4) and (A.11) depend on $\alpha$ (the
latter also on $r$) in less than immediately obvious ways. For (A.4), we
simply note that

(A.23)
$t^{\frac{2n+3+2\lceil\frac{\ell_{1}}{2}\rceil}{\alpha}}\leq\max\Bigl{(}t^{\frac{2n+3+2\lceil\frac{\ell_{1}}{2}\rceil}{\alpha_{-}}},t^{\frac{2n+3+2\lceil\frac{\ell_{1}}{2}\rceil}{\alpha_{+}}}\Bigr{)}\quad\forall
t>0,$

if $\alpha_{-}\leq\alpha\leq\alpha_{+}$, so it suffices to numerically
precompute

$\int_{0}^{\infty}|\log(t)|^{\ell_{2}}\max\Bigl{(}t^{\frac{2n+3+2\lceil\frac{\ell_{1}}{2}\rceil}{0.5+(i-1)\triangle}},t^{\frac{2n+3+2\lceil\frac{\ell_{1}}{2}\rceil}{0.5+i\triangle}}\Bigr{)}t^{-1}|p_{\ell_{2}}(t)|\mathrm{e}^{-t}\mathrm{d}t,$

for $\ell_{2}\in 1{:}4$, $\lceil\frac{\ell_{1}}{2}\rceil\in\\{0,1\\}$ and
$i\in 1{:}i^{*}$, and use each of these for in the respective square
$Q_{i,j}$. For numerical integration, we use the Julia library QuadGK.jl.

Concerning integrals of the form (A.11), we recall that parameter
$\beta_{1}=\pi+\alpha\varphi$ depends on $\alpha$, and note that it possible
to write

(A.24)
$\frac{\partial^{\ell_{2}}}{\partial\alpha^{\ell_{2}}}\bigl{[}\mathrm{e}^{\mathrm{i}(n+1)\beta_{1}}\tau^{(n+1)\alpha}M_{n+1}\bigl{(}\mathrm{e}^{\mathrm{i}\beta_{1}}\tau^{\alpha}\bigr{)}\bigr{]}=\tau^{(n+1)\alpha}\bigl{(}\log(\tau)+\mathrm{i}\varphi\bigr{)}^{\ell_{2}}\sum_{k=0}^{\ell_{2}}b_{n,k}\tau^{k\alpha}M^{(k)}\bigl{(}\mathrm{e}^{\mathrm{i}\beta_{1}}\tau^{\alpha}\bigr{)},$

where the coefficients $b_{n,k}$ depend on $\alpha$ but their absolute values
do not. For example, for $\ell_{2}=2$, the latter expression can be written as

$\displaystyle\bigl{(}\mathrm{e}^{\mathrm{i}\beta_{1}}\tau^{\alpha}\bigr{)}^{n+1}\bigl{(}\log(t)+\mathrm{i}\varphi\bigr{)}^{2}\Bigl{(}$
$\displaystyle(n+1)^{2}M\bigl{(}\mathrm{e}^{\mathrm{i}\beta_{1}}\tau^{\alpha}\bigr{)}-(2n+3)\mathrm{e}^{\mathrm{i}\alpha\varphi}\tau^{\alpha}M^{\prime}\bigl{(}\mathrm{e}^{\mathrm{i}\beta_{1}}\tau^{\alpha}\bigr{)}$
$\displaystyle\qquad+\bigl{(}\mathrm{e}^{\mathrm{i}\alpha\varphi}\tau^{\alpha}\bigr{)}^{2}M^{\prime\prime}\bigl{(}\mathrm{e}^{\mathrm{i}\beta_{1}}\tau^{\alpha}\bigr{)}\Bigr{)}.$

The crux of all this is that by applying the triangle inequality to (A.24) and
using simple uniform estimates for $M$ and its derivatives (which is possible
since $\mathrm{e}^{\mathrm{i}\beta_{1}}\tau^{\alpha}$ by construction always
has negative real part) , the integral appearing in (A.11) can be decomposed
as a linear combination of integrals of the form

$\int_{0}^{\infty}\big{|}\log(\tau)+\mathrm{i}\varphi\big{|}^{\ell_{2}}\tau^{k\alpha+\ell_{1}}\mathrm{e}^{\sin(\varphi)r\tau}\mathrm{d}t,\quad
n+1\leq k\leq n+\ell_{2}+1.$

Making the substitution $r\tau=\mathrel{\mathop{:}}t$, doing some elementary
massaging for the resulting integrand, taking $\varphi\to-\pi_{\alpha}$ and
using a discretization estimate similar to (A.23), we again end up with a
finite collection of integrals that can be precomputed.

The bivariate version of this estimate, (A.19), comes with the additional
ingredient of the Hankel functions, which we use the Julia library
SpecialFunctions.jl to compute.

All in all, we have described several different ways to bound

(A.25) $\sup_{(r,\alpha)\in Q_{i,j}}|f^{(\ell_{1},\ell_{2})}(r;\alpha)|,$

and for each square $Q_{i,j}$, we may take the smallest out of these bounds as
the ultimate bound (A.25).

It then remains to estimate

$\inf_{(r,\alpha)\in Q_{i,j}}f(r;\alpha),$

from below. For this purpose, we have first precomputed $f(r;\alpha)$ at the
corners of the $Q_{i,j}$’s using different numerical integration routines for
different ranges of the parameters; see Section 3.1 of the main paper, where
this is explained in the context of the spline approximation.

We then observe that $f(r;\alpha)$ is for fixed $\alpha$ a decreasing function
of $r>0$, both in the univariate and in the bivariate case. Thus, writing
$Q_{i,j}=\mathrel{\mathop{:}}[r_{-,j},r_{+,j}]\times[\alpha_{-,i},\alpha_{-,i}]$,
we may use the fundamental theorem of calculus to obtain

$\displaystyle\inf_{(r,\alpha)\in Q_{i,j}}f(r;\alpha)$
$\displaystyle=\inf_{\alpha\in[\alpha_{-,i},\alpha_{-,i}]}f(r_{+,j};\alpha)$
$\displaystyle\geq\min\bigl{(}f(r_{+,j};\alpha_{-,i}),f(r_{+,j};\alpha_{+,i})\bigr{)}-\frac{\triangle}{2}\sup_{(r,\alpha)\in
Q_{i,j}}|f^{(0,1)}(r;\alpha)|,$

where the latter can further be estimated using the bounds discussed in the
context of (A.25). It turns out that with a small enough $\triangle$, such as
$\triangle=\frac{2}{10^{3}}$, this yields a fair lower bound for a function
like $f(r;\alpha)$ which is otherwise very difficult to bound from below.

## Tail error bounds

Here we establish relative error bounds for the series decompositions (A.6)
and (A.14), used in the main paper with $n=3$ and $r>30$.

First, in the univariate case, we can rewrite (A.6) as

$\displaystyle f(r;\alpha)$
$\displaystyle=\sum_{k=1}^{3}\frac{(-1)^{k+1}\Gamma(k\alpha+1)\sin(\frac{k\alpha\pi}{2})}{\pi
k!}r^{-k\alpha-1}+\mathcal{R}_{4}(r;\alpha)$
$\displaystyle\qquad\qquad=\mathrel{\mathop{:}}\sum_{k=1}^{3}c^{\mathcal{S}}_{k}(\alpha)r^{-k\alpha-1}+\mathcal{R}_{4}(r;\alpha)=\mathrel{\mathop{:}}\mathcal{S}_{3}(r;\alpha)+\mathcal{R}_{4}(r;\alpha),$

with

$|\mathcal{R}_{4}(r;\alpha)|\leq\frac{\Gamma(4\alpha+1)}{24\pi\sin(\pi_{\alpha})^{4\alpha+1}}r^{-4\alpha-1}\mathrel{\mathop{:}}=c^{\mathcal{R}}_{4}(\alpha)r^{-4\alpha-1}.$

Thus, for $r\geq 30$, we have the following estimates:

$\displaystyle\Big{|}1-\frac{f(r;\alpha)}{\mathcal{S}_{3}(r;\alpha)}\Big{|}=\Big{|}\frac{\mathcal{R}_{4}(r;\alpha)}{\mathcal{S}_{3}(r;\alpha)}\Big{|}$
$\displaystyle\leq\frac{c^{\mathcal{R}}_{4}(\alpha)r^{-4\alpha-1}}{|c^{\mathcal{S}}_{1}(\alpha)|r^{-\alpha-1}-|c^{\mathcal{S}}_{2}(\alpha)|r^{-2\alpha-1}-|c^{\mathcal{S}}_{3}(\alpha)|r^{-3\alpha-1}}$
$\displaystyle=\frac{c^{\mathcal{R}}_{4}(\alpha)r^{-3\alpha}}{|c^{\mathcal{S}}_{1}(\alpha)|-|c^{\mathcal{S}}_{2}(\alpha)|r^{-\alpha}-|c^{\mathcal{S}}_{3}(\alpha)|r^{-2\alpha}},$

assuming these calculations are valid in the sense of the latter denumerator
being positive for $r=30$ – numerical considerations show that this is indeed
the case. The rightmost quantity is then a decreasing function of $r\geq 30$,
and we thus get a uniform error bound by investigating it with $r=30$:

$\sup_{r>30,\;0.5\leq\alpha\leq
1.9}\,\Big{|}1-\frac{f(r;\alpha)}{\mathcal{S}_{3}(r;\alpha)}\Big{|}\leq
0.00096.$

In the bivariate case, we may proceed similarly, this time with (A.14) for
$n=3$ and bounding the remainder term with (A.18). We can estimate the
integral in (A.18) using grid-based numerical upper bounds like in the
previous section. We thus find

$\sup_{r>30,\;0.5\leq\alpha\leq
1.9}\,\Big{|}1-\frac{f(r;\alpha)}{\mathcal{S}_{3}(r;\alpha)}\Big{|}\leq
0.0016.$

We may then use the logarithm function’s basic continuity properties near $1$
to infer

$\sup_{r>30,\;0.5\leq\alpha\leq
1.9}\,|\log\,f(r;\alpha)-\log\,\mathcal{S}_{3}(r;\alpha)|\leq\begin{cases}0.00097&\quad\textrm{(univariate
case);}\\\ 0.0017&\quad\textrm{(bivariate case).}\end{cases}$

## References

  * [1] A. Achim, A. Basarab, G. Tzagkarakis, P. Tsakalides and D. Kouamé  “Reconstruction of Ultrasound RF Echoes Modeled as Stable Random Variables”  In _IEEE Transactions On Computational Imaging_ 1, 2015, pp. 86–95
  * [2] A. K. Alekseev, I. M. Navon and J. L. Steward  “Comparison of advanced large-scale minimization algorithms for the solution of inverse ill-posed problems”  In _Optimization Methods and Software_ 24.1 Taylor & Francis, 2009, pp. 63–87
  * [3] S. Ament and M. O’Neil  “Accurate and efficient numerical calculation of stable densities via optimized quadrature and asymptotics”  In _Statistics And Computing_ 28, 2018
  * [4] George E. Andrews, Richard Askey and Ranjan Roy  “Special Functions”, Encyclopedia of Mathematics and its Applications  Cambridge University Press, 1999  DOI: 10.1017/CBO9781107325937
  * [5] A. Arjas, L. Roininen, M. J. Sillanpää and A. Hauptmann  “Blind Hierarchical Deconvolution”  In _2020 IEEE 30th International Workshop on Machine Learning for Signal Processing (MLSP)_ , 2020, pp. 1–6  DOI: 10.1109/MLSP49062.2020.9231822
  * [6] I Belov  “On the computation of the probability density function of $\alpha$-stable distributions”  In _Mathematical Modelling and Analysis_ , 2005
  * [7] Igoris Belovas, Leonidas Sakalauskas, Vadimas Starikovicius and Edward Sun  “Mixed-Stable Models: An Application to High-Frequency Financial Data”  In _Entropy_ 23, 2021, pp. 739  DOI: 10.3390/e23060739
  * [8] Harald Bergström  “On some expansions of stable distribution functions”  In _Arkiv för Matematik_ 2.4 Institut Mittag-Leffler, 1952, pp. 375 –378
  * [9] D. Bolin  “Spatial Matérn fields driven by non-Gaussian noise.”  In _Scandinavian Journal of Statistics_ 41.3, 2014, pp. 557–579
  * [10] Szymon Borak, Wolfgang Härdle and Rafał Weron  “Stable Distributions”  In _Statistical Tools for Finance and Insurance_ Berlin, Heidelberg: Springer Berlin Heidelberg, 2005, pp. 21–44  DOI: 10.1007/3-540-27395-6˙1
  * [11] Arias Calluari, K. Alonso-Marroquin and M. Harré  “Closed-form solutions for the Lévy-stable distribution”  In _Physical Review E_ , 2018
  * [12] R. Carlson and C. Hall  “Error bounds for bicubic spline interpolation”  In _Journal Of Approximation Theory_ 7, 1973, pp. 41–47
  * [13] N. K. Chada, M. Iglesias, L. Roininen and A. M. Stuart  “Parameterizations for ensemble Kalman Inversion.”  In _Inverse Problems_ 34.5, 2019
  * [14] N. K. Chada, C. Schillings and S. Weissmann  “On the incorporation of box-constraints for ensemble Kalman inversion.”  In _Foundations of Data Science_ 4.1, 2019
  * [15] N. K. Chada, A. M. Stuart and X. T. Tong  “Tikhinov regularization within ensemble Kalman inversion.”  In _SIAM Journal on Numerical Analysis_ 55.2, 2019
  * [16] Laurent Condat  “Discrete Total Variation: New Definition and Minimization”  In _SIAM Journal on Imaging Sciences_ 10, 2017, pp. 1258–1290
  * [17] J. Crisanto-Neto, M. Luz, P. E. Raposo and G Viswanathan  “An efficient series approximation for the Lévy $\alpha$ -stable symmetric distribution”  In _Physics Letters A_ 382, 2018
  * [18] A. C. Damianou and N. D. Lawrence  “Deep Gaussian processes”  In _AISTATS_ , 2013, pp. 207–2015
  * [19] Austen C. Duffy  “An Introduction to Gradient Computation by the Discrete Adjoint Method” Preprint at http://computationalmathematics.org, 2009
  * [20] M. M. Dunlop, M. A. Girolami, A. M. Stuart and A. Teckentrup  “How deep are deep Gaussian processes?”  In _Journal of Machine Learning Research_ 19.54, 2018, pp. 1–46
  * [21] Kuruoglu E. E. and Zerubia J.  “Modeling SAR images with a generalization of the rayleigh distribution”  In _IEEE transactions on image processing_ 13 Institute of ElectricalElectronics Engineers, New York, 2004, pp. 527–533  DOI: 10.1109/tip.2003.818017
  * [22] Xiaohu Ge, Guangxi Zhu and Yaoting Zhy  “On the testing for alpha-stable distributions of network traffic”  In _Computer Communications_ 27, 2004, pp. 447–457
  * [23] Dan Givoli  “A tutorial on the adjoint method for inverse problems”  In _Computer Methods in Applied Mechanics and Engineering_ 380, 2021, pp. 113810
  * [24] C. A. Hall and W. W. Meyer  “Optimal error bounds for cubic spline interpolation.”  In _Journal of Approximation Theory_ 16.2, 1976, pp. 105–122
  * [25] Mohd Asrul Hery Ibrahim and Wah Leong  “The Hybrid BFGS-CG Method in Solving Unconstrained Optimization Problems”  In _Abstract and Applied Analysis_ 2014, 2014, pp. 6  DOI: 10.1155/2014/507102
  * [26] J. Kaipio and E. Somersalo  “Statistical and computational inverse problems”  Springer-Verlag, 2004
  * [27] Vassili N. Kolokoltsov  “Markov Processes, Semigroups and Generators”  Berlin, New York: De Gruyter, 2011  DOI: doi:10.1515/9783110250114
  * [28] F. Kühn  “Lévy Matters VI: Lévy-Type Processes: Moments, Construction and Heat Kernel Estimates”, Lecture Notes in Mathematics  Springer International Publishing, 2017
  * [29] Matti Lassas, Eero Saksman and Samuli Siltanen  “Discretization-invariant Bayesian inversion and Besov space priors”  In _Inverse Problems & Imaging_ 3.1, 2009, pp. 87–122
  * [30] Matti Lassas and Samuli Siltanen  “Can one use total variation prior for edge-preserving Bayesian inversion?”  In _Inverse Problems_ 20.5 IOP Publishing, 2004, pp. 1537–1563  DOI: 10.1088/0266-5611/20/5/013
  * [31] Chen Li, Matthew Dunlop and Georg Stadler  “Bayesian neural network priors for edge-preserving inversion”  In _Inverse Problems and Imaging_ , 2022
  * [32] F. Lindgren, H. Rue and J. Lindström  “An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach.”  In _J. R. Stat. Soc. Ser. B Stat. Methodol_ 73, 2011, pp. 423–498
  * [33] Muneya Matsui and Akimichi Takemura  “Some Improvements in Numerical Evaluation of Symmetric Stable Density and Its Derivatives”  In _Communications in Statistics - Theory and Methods_ 35, 2004, pp. 149–172
  * [34] Christian Menn and Svetlozar T. Rachev  “Calibrated FFT-based density approximations for $\alpha$-stable distributions”  In _Computational Statistics & Data Analysis_ 50.8, 2006, pp. 1891–1904  DOI: https://doi.org/10.1016/j.csda.2005.03.004
  * [35] S. Mercan and M. Alam  “Anomaly detection in hyperspectral imagery using Stable Distribution”  In _Proceedings Of SPIE - The International Society For Optical Engineering_ 8049, 2011
  * [36] S. Mittnik, T. Doganoglu and D. Chenyao  “Computing the probability density function of the stable Paretian distribution”  In _Mathematical And Computer Modelling_ 29, 1999, pp. 235–240
  * [37] Patrick Kofod Mogensen and Asbjørn Nilsen Riseth  “Optim: A mathematical optimization package for Julia”  In _Journal of Open Source Software_ 3.24, 2018, pp. 615  DOI: 10.21105/joss.00615
  * [38] Ludovic Métivier, Romain Brossier, Stéphane Operto and Jean Virieux  “Full Waveform Inversion and the Truncated Newton Method”  In _SIAM Journal on Scientific Computing_ 35, 2013, pp. B401–B437  DOI: 10.1137/120877854
  * [39] R. M. Neal  “Bayesian Learning for Neural Networks”  Springer-Verlag, 1996
  * [40] C. Nikias and M. Shao  “Signal processing with alpha-stable distributions and applications”, 1995
  * [41] “NIST Digital Library of Mathematical Functions” F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, eds., http://dlmf.nist.gov/, Release 1.1.7 of 2022-10-15  URL: http://dlmf.nist.gov/
  * [42] J. Nolan  “An algorithm for evaluating stable densities in Zolotarev’s (M) parameterization”  In _Mathematical And Computer Modelling_ 29, 1999, pp. 229–233
  * [43] J. Nolan  “Multivariate elliptically contoured stable distributions: theory and estimation”  In _Computational Statistics_ 28, 2013, pp. 2067–2089
  * [44] J. Nolan  “Numerical calculation of stable densities and distribution functions”  In _Communications In Statistics. Stochastic Models_ 13, 1997, pp. 759–774
  * [45] C. E. Rasmussen and C. K. I. Williams  “Gaussian Processes for Machine Learning”  MIT Press, 2006
  * [46] L. Roininen, J. M. J. Huttunen and S. Lasanen  “Whittle-Matérn priors for Bayesian statistical inversion with applications in electrical impedance tomography”  In _Inverse problems and Imaging_ 8, 2014
  * [47] G. Samorodnitsky and M. S. Taqqu  “Stable non-Gaussian Random Processes”  Chapman & Hall, 1994
  * [48] K. Sato  “Lévy Processes and Infinitely Divisible Distributions”  Cambridge University Press, 1999
  * [49] M. Scherer, S. Rachev, S. Young and F. A Fabozzi  “FFT-based approximation of tempered stable and tempered infinitely divisible distributions” Preprint at https://statistik.econ.kit.edu/, 2009
  * [50] J. Suuronen, N. K. Chada and L. Roininen  “Cauchy Markov Random Field Priors for Bayesian Inversion”  In _Statistics and Computing_ 32, 2022, pp. 33
  * [51] A. Tarantola  “Inverse Problem Theory and Methods for Model Parameter Estimation”  Elsevier, 1987
  * [52] Mahdi Teimouri and Hamidreza Amindavar  “A Novel Approach to Calculate Stable Densities”  In _Proceedings of the World Congress on Engineering_ 1, 2008
  * [53] Nan Ye, Farbod Roosta-Khorasani and Tiangang Cui  “Optimization Methods for Inverse Problems”  In _2017 MATRIX Annals_ Springer International Publishing, 2019, pp. 121–140
  * [54] Zheng Zhao, Muhammad Emzir and Simo Särkkä  “Deep state-space Gaussian processes”  In _Statistics and Computing_ 31, 2021
  * [55] Ciyou Zhu, Richard H. Byrd, Peihuang Lu and Jorge Nocedal  “Algorithm 778: L-BFGS-B: Fortran Subroutines for Large-Scale Bound-Constrained Optimization”  New York, NY, USA: Association for Computing Machinery, 1997, pp. 550–560