  We propose to use L\'evy {\alpha}-stable distributions for constructing
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 modelling rough features. The main technical issue is that the
{\alpha}-stable probability density functions do not have closed-form
expressions in general, and this limits their applicability. For practical
purposes, we need to approximate probability density functions through
numerical integration or series expansions. Current available approximation
methods are either too time-consuming or do not function within the range of
stability and radius arguments needed in Bayesian inversion. To address the
issue, we propose a new hybrid approximation method for symmetric univariate
and bivariate {\alpha}-stable distributions, which is both fast to evaluate,
and accurate enough from a practical viewpoint. Then we use approximation
method in the numerical implementation of {\alpha}-stable random field priors.
We demonstrate the applicability of the constructed priors on selected Bayesian
inverse problems which include 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. We employ maximum-a-posterior-based estimation at all
the numerical examples. To that end, we exploit the limited-memory BFGS and its
bounded variant for the estimator.
