Cauchy difference priors for edge-preserving Bayesian inversion

2.1 Lévy α-stable, Cauchy and Gaussian random walk priors

Our main aim in this paper is doing numerics, however here we will introduce the basic, and more formal, concepts of the processes we are interested in, as we feel that this is important background information for the benefit of the reader. But more importantly this sets the direction of the future studies, where more technical results will be done in order to prove the discretization-invariance of the posterior distribution.

Let {𝒳⁢(t),t∈[0,T]} be a stochastic process for some 0<T<∞. Furthermore, let Sα⁢(τ,β,μ) be Lévy α-stable distribution with parameters α, τ, β and μ, for some 0<α≤2,-1≤β≤1. We call 𝒳⁢(t) a continuous Lévy α-stable process starting from zero, if 𝒳⁢(0)=0, 𝒳 has independent increments and

for any 0≤s<t≤T. For more details on Lévy α-stable processes, see [37]. In order to illustrate the choices of the parameters, in Figure 1, we have plotted certain α-stable probability density functions with varying α and β. The case α=2 is the Gaussian random walk, and it is a special case and high end of stable distributions. The Cauchy walk is obtained with α=1. Parameter β models the skewness of the probability density, a useful feature in many applications, but here we mostly concentrate on symmetric densities hence β=0.

Figure 1

Probability density functions with varying α and β parameters.

For example, α-stable processes can be constructed by using independently scattered measures. A random measure M on [0,T] maps (Lebesgue) measurable sets A⊂[0,T] to random variables M⁢(A) in such a way that M⁢(⋃kAk)=∑kM⁢(Ak) almost surely for disjoint measurable sets Ak⊂𝕀. A random measure M is called independently scattered if M⁢(A1),…,M⁢(An) are independent random variables whenever the measurable sets Ai⊂[0,T], i=1,…,n, n∈ℕ, are pairwise disjoint. We require that

with constant skewness β. In our case, the measure m⁢(⋅) on [0,T] is the Lebesgue’s measure |⋅| multiplied with a constant λ>0. Then α-stable Lévy motion starting from zero is

where χA is the characteristic function of the set A.

For the continuous limit of the Lévy α-stable random walk, we apply independently scattered measures. Let us denote the discrete random walk at t=j⁢h by Xj, where j∈ℤ+ and h>0 is the discretization step. We obtain a random walk approximation (see in [37, Section 3.3])

It is easy to see that such a random walk approximation converges to the continuous Lévy α-stable motion as h→0 in distribution in the Skorokhod space of functions that are right-continuous and have left limits.

Our special interest is the Cauchy random walk. Given equation (2.4) and choosing α=1, β=0, we can write the joint distribution of the Cauchy random walk as a probability density

where λ is the regularizing constant and C is a normalization constant. We can apply this probability density directly as an a priori probability density in Bayes’ formula in equation (2.3).

In Figure 2, we have plotted realizations of Cauchy and Gaussian random walks. The latter is the Lévy α-stable random walk with α=2. The Cauchy random walk realizations have either small perturbations or big jumps. Hence, intuitively it feels plausible that using Cauchy distributions for differences lead to edge-preserving inversion. Gaussian random walk realizations do not prefer jumps, but continuous paths, hence the tendency for smoothing in Bayesian inversion.

Figure 2

Realizations of Cauchy and Gaussian random walks.

(a)

Cauchy random walk.

(b)

Gaussian random walk.

2.2 Edge-preserving inversion

To our best knowledge, for difference priors there has not been any simple criteria which would indicate, whether the prior favors edge-preserving solutions instead of smooth ones. The intuitive idea is to construct a prior, which promotes discontinuities, i.e. jumps. Hence, via a negation, we may say that in edge-preserving inversion, the priors do not promote smoothness everywhere.

Here, for simplicity, we will first consider a one-dimensional case and take h=λ=1. We will generalize the construction for two-dimensional case in the next sections. Hence, let us look at three consecutive coordinates Xj-1,Xj,Xj+1 of a random vector X. We could then say that X prefers smoothness at point j if Xj is (most probably) in the middle of Xj-1 and Xj+1, and that there is edge at point j if Xj is (most probably) either at Xj-1 or at Xj+1. The above mentioned idea means that we are interested in the conditional prior distribution of Xj on the condition that Xj-1 and Xj+1 are known. Without loss of generality we can then assume that Xj-1=-a and Xj+1=a. The conditional density for Cauchy difference prior for Xj can then be written as a probability density

This is simply a product of two Cauchy probability density functions, one from each neighbor of Xj. In order to see the properties of these functions, consider the second derivative of the probability density function with respect to Xj at zero:

In Figure 3, we have plotted these probability density functions. With |a|<1 we have a unimodal function, hence regularization promotes smooth solution around Xj=0. With |a|>1 we have a bimodal function, hence regularization promotes solutions around Xj=±a with smooth variations around these two values. This means that this kind of prior promotes jumps or small local oscillations.

Figure 3

Upper panel: Cauchy probability density function for X2, given fixed X1=-a and X3=a. Bottom panel: The same case for Gaussian difference prior and total variation prior.

(a)

Unimodal Cauchy a<1.

(b)

Flat Cauchy with |a|=1.

(c)

Bimodal Cauchy a>1.

(d)

Gaussian.

(e)

Total variation.

We can make similar plots for Gaussian difference prior and total variation prior. The Gaussian difference prior is simply a Gaussian-shaped function with maximum at Xj=0. Hence, providing smoothing properties. The total variation priors is constant for [-a,a] and decays exponentially outside this region. Hence, total variation neither punishes discontinuities or promotes edge-preserving inversion as it promotes all values equally between these two points. This is quite close the case for the Cauchy prior with |a|=1, i.e. in a sense the shape of the density function is similar to the total variation prior density function. The effect on reconstruction for both α→0 and α=2-ϵ can be understood with similar arguments. For α→0, the modalities become spikes, while for α=2-ϵ is bimodal, and becomes trimodal and finally unimodal for very small ϵ. If α=2, this can be considered as a special case.

2.3 Two-dimensional difference priors

Let us consider a two-dimensional lattice (h⁢j,h′⁢j′), where (j,j′)∈𝕀2⊂ℤ2 and discretization steps to each coordinate directions h,h′>0. A Gaussian difference prior can be constructed via difference equations with i.i.d. increments

The term σ2>0 is the regularization parameter. The prior is then constructed through products of conditional normal distributions

where L1 and L2 are difference matrices to two different coordinate directions, and C1=σ2⁢hh′⁢I and C2=σ2⁢h′h⁢I, where I is an identity matrix. These constructions are rather well-known, see for example [23]. We emphasize that equation (2.6) is not explicitly solvable, as only the inverse covariance matrix exists. For uniqueness of the least squares solution, we would need to impose e.g. boundary conditions.

Let us consider, with a simple example, why we scale the variances by hh′ and h′h. It is enough to consider the first difference equation in equation (2.6), as the other one is obtained through same arguments. We consider the grid densification in two steps: (1) We consider denser discretization along j-direction, which we refer as the regularization coordinate direction, and (2) denser discretization along j′-direction. Let Xj,j′-Xj-1,j′∼𝒩⁢(0,σ2). If we densify the lattice along the regularization coordinate direction j, this would correspond to one-dimensional Gaussian random walk and hence we need to add the discretization step in the variance, hence Xj,j′-Xj-1,j′∼𝒩⁢(0,σ2⁢h).

The second coordinate direction is a bit trickier. Let us consider a simplified situation as depicted in Figure 4. Thus, we define pixelwise mean

The reason for scaling by 2 is because this is the standard pixelwise arithmetic mean, i.e. we want to have natural scaling of pixels. Now the increments

If we choose X2-X1∼𝒩⁢(0,σ2), then we must choose X21-X11∼X22-X12∼𝒩⁢(0,2⁢σ2). Hence, from this follows the scaling by 1h′ in equation (2.6).

As an alternative, it could be appealing to use a stochastic partial differential equation Δ12⁢𝒳=𝒲, where Δ is the Laplacian. However, equivalently in the sense of probability distributions, you may consider a system of equations written as

whose inverse covariance operator in the least squares sense is the Laplacian Δ. We prefer to use formulation in equation (2.7) instead of Δ12⁢𝒳=𝒲, as working with fractional order Laplacian Δ12 is quite tedious numerically.

Similarly for the Cauchy differences we may write with i.i.d. increments

The result is obtained through similar argumentation as in the Gaussian case. If we make the grid denser along the regularization direction, this corresponds again to one-dimensional Cauchy random walk, i.e. we simply scale with discretization step h.

In the second coordinate direction we take again the same division as in Figure 4. If X2-X1∼Cauchy⁢(λ), then X21-X11∼X22-X12∼Cauchy⁢(λ) (note independent linear combination of Cauchy distributed random variable).

Through similar arguments, we can actually see that the scaling of the general Lévy α-stable difference prior is with i.i.d. increments

This, naturally, implies similar constructions in the higher dimensions.

Figure 5

Realizations of two-dimensional Cauchy, Gaussian and TV priors exhibiting blocky and smooth structures.

(a)

Realization of 2D Cauchy prior.

(b)

Realization of 2D Gaussian prior.

(c)

Realization of 2D TV prior.

We note that the prior density is not exponentially bounded. Therefore, special care is needed in estimation whenever the direct theory A in (2.2) has nontrivial null space.

In Figure 5, we have plotted realizations of two-dimensional Gaussian and Cauchy priors, defined in equations (2.6) and (2.8), respectively. We have also plotted a realization of a total variation prior. We have used zero-boundary conditions, but we have cropped the image in order to demonstrate blocky or smooth structures of the realizations. A notable feature in the Cauchy prior realization is that it looks blocky, and there are distinct sharp edges. For drawing the realizations, we use single-component Metropolis–Hastings, which is explained in next section, i.e. we generate a chain of realizations of the Cauchy prior.

We note that Cauchy priors are not rotation-invariant, and this can be clearly especially from the realization. This can explained through as a packing problem, i.e. in order to have rotation-invariant Cauchy priors, the edges should be spherical. But we cannot efficiently pack balls, hence the box-shaped realizations. In this current construction, the box-shapes are restricted to the coordinate directions. This can be changed by making a coordinate transformation and hence use directed nabla operator, but this would give us more complicated formulas. It would be indeed an interesting idea to use the directed derivatives, and let the direction be modeled in the hyperprior. Similarly TV prior is also neither rotation-invariant in its standard formulation, but it can be made by rotation-invariant by a reformulation as in [15]. The anisotropic TV prior, in a sense, has coordinate-wise edgy features as the Cauchy prior.

If we consider the edge-preserving features of the two-dimensional Cauchy prior, we can study processes along either of the coordinate directions. These one-dimensional processes look practically similar to the one-dimensional Cauchy random walks. For comparison purposes, we have also plotted realization of a Gaussian difference prior and TV priors obtained in a similar manner. As we can see, the Gaussian distribution promotes smoothness, not jumps. The realization of a TV prior is not as blocky as the Cauchy prior, but resembles more the Gaussian prior.

4.1 One-dimensional deconvolution

In Figure 6, we have a synthetic deconvolution example with Cauchy and total variation priors on different grids. The original unknown is a piecewise constant function, and we observe convolved noisy measurements. The constant convolution kernel is cosine-shaped. In the measurements, we have added white noise with standard deviation σ=0.3, which corresponds to ∼10% relative ratio.

Figure 6

The effect the discretisation level on the CM-estimate with Cauchy and total variation priors: Top panel: noisy deconvolved observations (solid line). Middle panel: Cauchy prior CM-estimates on different lattices (solid lines), bottom panel – The same as middle panel, but with total variation priors. The dashed line in all figures corresponds to the true target and the gray area corresponds to 2⁢σ credibility intervals.

(a)

Measurements.

(b)

Cauchy N=131.

(c)

Cauchy N=261.

(d)

Cauchy N=521.

(e)

Total variation N=131.

(f)

Total variation N=261.

(g)

Total variation N=521.

We have made the reconstructions with a number of unknowns N=131,261,521. The Cauchy regularization parameter λ in equation (2.5) was chosen in such a way that reconstructions look edge-preserving. Some care needs to be taken in choosing λ, similarly to using any Gaussian priors, in order not to use either to too strong or loose priors – this demands some tuning, and in the future studies this parameter should be modeled as a hyperparameter similarly to e.g. [11], where Gaussian prior parameters where modeled in the hyperprior. On denser meshes, we simply then change regularization parameter by scaling with discretization step h accordingly. The reconstructions on different grids with Cauchy priors clearly provide edge-preserving features. We have also plotted the 2⁢σ credibility intervals, and these look also rather similar with different N except in the jumps, where edges do not see to converge but rather grow. We suspect that this is because the length of the chain is too short. The posterior is discretization-invariant in the sense of weak convergence of measures. A simple example of the posterior of the random variable X2 given Y1=Xj-1+ε and Y2=Xj+1+ε, when Y1=Y2=X0=0, shows that

Thus the posterior variance of Xj is bounded even though we have only observations for the neighboring points. Similarly we scale the TV prior regularization parameter. We note that the discretization-invariance behavior of estimators cannot be obtained with total variation priors, as can be seen from the smoothing behavior when we increase the number of unknowns. For further discussion of the total variation prior examples, see Lassas and Siltanen 2004 [25].

In the single-component Metropolis–Hastings, we use 1,000,000 long chain in order to guarantee convergence. We need long chains in the single-component algorithm, as we are dealing with multimodal distributions, and the single-component algorithms are known get stuck in one mode and it takes a long run to get to the correct mode. Burn-in period is chosen to be half of the chain length for both total variation and Cauchy cases. The MCMC chain could be optimized, but as in this paper our goal is to simply demonstrate the Cauchy prior, we postpone the MCMC optimization studies for subsequent papers.

4.2 Two-dimensional X-ray tomography

Now we will apply the developed method to X-ray tomography. We take the Shepp–Logan phantom to be the original unknown. The measurements correspond to two-dimensional fan-beam geometry (see Figure 7 for depicted the geometry of the measurements). The distance between the source and the detector to the rotation axis is 4 and 2, respectively. The width of the detector is 3 and the detector element contains 200 pixels. We make measurements at 41 different angles (from -10∘ to 190∘ with 5∘ steps). In addition, we add white measurement noise with approximately 0.1 % noise level.

Figure 7

Measurement geometry for the synthetic two-dimensional tomography example.

We compute CM estates using single-component Metropolis–Hastings as described in Section 3. The chain is initiated from the back-projection estimate (AT⁢m scaled properly). In order to get convergence of the MCMC chain, we use 5 million as the chain length from which 2.5 million samples were neglected as the burn-in period. In order to save memory, we store only every tenth sample for the computation of the ensemble means. Here, chain length refers to the number of single component updates, hence in the highest resolution in 2D simulations we have 800×800×5×106=3.2×1012 samples. However, we do not need to evaluate the full likelihood, as we can simply use the sparsity of the forward theory matrix, and the iid property of the difference priors – this allows realistic computation times. Computation time for a reconstruction with 200×200 resolution is approximately 24 hours and, for a 400×400 resolution, the time is approximately 96 hours. This indicates that the algorithm’s computational cost is (nearly) proportional to the number of unknowns, as expected. The computations are carried with a PC workstation with dual Xeon E5-2680 CPUs and 128 GB memory. The core parts of the algorithm are implemented using the MEX/C interface.

For comparison, we also compute MAP estimates using Cauchy prior which is computationally less expensive. The MAP estimate is computed by minimizing the loss function

where the additional penalty term is added to ensure non-negativity of the solution to improve the stability of the optimization. The minimization is carried out by employing Gauss-Newton optimization such that μ and ϵ are decreased during the iteration. We note that, as being non-convex optimization problems, the Gauss–Newton iteration has tendency can converge to a (suboptimal) local minima which also seems to be case here (see below).

Figure 8

Comparison of different estimates on lattice sized 400×400 pixels.

(a)

Unknown: Shepp–Logan phantom.

(b)

Filtered back-projection estimate.

(c)

MAP estimate with Cauchy prior.

(d)

CM estimate with Cauchy prior.

(e)

CM estimate with TV prior.

(f)

CM estimate with Gaussian prior.

In Figure 8, we have plotted the original unknown and reconstructions with five different methods. The filtered back-projection estimate recovers most of the features, but it has severe artifacts. The CM estimate with Cauchy prior captures the features well without significant artefacts. The CM estimate with TV prior captures also the features well, but it is not as sharp as the Cauchy prior estimates. Also, the Cauchy estimate seems to be more robust, when we decrease the number of measurements or increase the noise. The estimate with Gaussian prior recovers the main structures, but it has both bad artifacts and also smooth structures.

Figure 8 also shows the MAP estimate computed using the Cauchy prior. This estimate also captures the features well, while the CM estimate has less artifacts in the reconstructions. However, as noted above, Gauss-Newton iteration has tendency to converge to (suboptimal) local minimas and that seems to be also the case here: starting the iteration from the CM estimate gives an MAP estimate with lower cost L0 (i.e. higher posterior probability) and this estimate is visually indisquistable from the CM estimate. But of course, starting from the CM is not a practically very meaningful approach. However, the finding global (optimal) minimas of a non-convex cost function is commonly known as a difficult task and, since our main motivation is in MCMC estimates, finding optimal MAP estimates is left to future studies.

Figure 9

Conditional mean estimates of the X-ray tomography problem on three different lattices.

(a)

200×200.

(b)

400×400.

(c)

800×800.

Figure 10

Cross-sections of the conditional mean estimates from the middle of the computational domain.

(a)

Horizontal.

(b)

Horizontal zoom-in.

(c)

Vertical.

(d)

Vertical zoom-in.

In Figures 9 and 10, we have CM estimates on different size lattices. The discretization has been taken into account as described earlier in this section, and hence the reconstructions look essentially same.