Here we illustrate our work on computational UQ for inverse problems with a few examples from ongoing projects at DTU Compute. These examples illustrate various methodologies for implementing, applying and using UQ.

Tikhonov Regularization with Nonnegativity Constraints

*Johnathan M. Bardsley, Univ. of Montana and Per Christian Hansen*

Nonnegativity constraints appear in many applications where the underlying physics dictates that the solution cannot be negative. This is the case, e.g., for absorption coefficients in computed tomography, image intensities in astronomical imaging, and wave velocities in seismic travel-time tomography. Hence, nonnegativity constraints are very important for producing meaningful reconstructions,

We considers computational aspects of UQ for Tikhonov regularization problems with nonnegativity constraints, formulated in a Bayesian setting:

min_{x} || **A** **x** – **b** ||_{2}^{2} + λ || **x** ||_{2}^{2} subject to * x* ≥

If we replace the nonnegativity constraint with a positivity constraint, * x* >

The above figure shows a simple Gauss distribution (left), a truncated Gauss distribution with zero probability for *x* = 0 (middle), and a nonnegative Gaussian distribution with nonzero probability for *x* = 0 (right).

This figure shows a histogram of a simple 2D nonnegativity constrained least squares problem (i.e., α = 0) whose exact solution is (*x*_{1},*x*_{2}*)* = (0,1), associated with a 2D nonnegative Gauss distribution. Notice the significant "mass" of the histogram for *x*_{1} = 0. Except for very simple cases the corresponding probability density function for the posterior does not have a convenient analytical expression. Instead we use a Markov chain Monte Carlo (MCMC) method to sample the constrained regualrized solution's probability distribution, and we refer to our algorithm as a *Nonnegative Hierarchical Gibbs Sampler*.

We applied our algorithm to positron emission tomography (PET), where a radioactive tracer element is injected into a body and exhibits radioactive decay, resulting in photon emission. The emitted photons that leave the body are recorded by a photon detector, and the resulting mathematical reconstruction problem is equivalent to a classical X-ray CT problem. The above figure shows our simulation results.

- Top left: the mean of the computed image samples from our nonnegative hierarchical Gibbs sampler.
- Bottom left: histogram of the regularization parameters produced by the sampler.
- Top right: the MAP estimate using the mean of the regularization parameters.
- The pixel-wise standard deviations from the sampler. As expected, image pixels with higher intensity also have higher variance.

For details see: J. M. Bardsley and P. C. Hansen, *MCMC algorithms for computational UQ of nonnegativity constrained linear inverse problems*, in review.

*Mirza Karamehmedović*

This example illustrates UQ applied to a stochastic model of a physical phenomenon. Specifically, we consider a stochastic differential equation that models a simple electric direct-current problem which can be considered a "1D electrical impedance tomography (EIT) problem."

The expected resistivities *ρ*_{a} and *ρ*_{b }of the inclusion and the ambient are assumed to be known. The length *a* of the inclusion is estimated in terms of a probability density, based on measurements of potential *V _{m}* at

*dV _{x} = I ρ dx + σ dW_{x} .*

Here the ambient resistivity *ρ* is the constant *ρ*_{a} or *ρ*_{b}, depending on the considered interval. The constant *σ* is the same in the two intervals for *x*. Moreover, *dW _{x}* is a white-noise term that introduces the stochasticity. Through the use of a forward Fokker-Planck equation, we show that the corresponding model sampling density is

*p*(*V _{m}* |

and then Bayes' theorem leads to the posterior

*p*(*a* | *V _{m}*) = const ⋅

Here, *V _{m}* is the measured data,

The above animation shows the evolution of the posterior *p*(*a* | *V _{m}*) as we obtain an increasing number of observations. In this example the prior

This is ongoing work at DTU Compute.

*Nicolai A. B. Riis, Yiqiu Dong, and Per Christian Hansen*

This example illustrates the use of UQ in a case where the forward model has uncertain model-parameters. Specifically we consider a Computed Tomography (CT) problem where the projection angles are not known precisely - but we have a statistical model for them.

Our goal is to compute the unknown attenuation coefficient image, represented by the vector ** x**, from measured and noisy data represented by the vector

** b** =

where the matrix ** R**(

* b* =

where the *model discrepancy term* ** d** acquires its push-forward distribution from

** d** =

Inspired by work by Kaipio & Somersalo we marginalize both ** d** and

*π*(** b**|

We now assume that both ** d**|

− log *π*(** b**|

where *L*^{T}** L** =

min_{x} || ** L** (

We developed a new iterative algorithm MD-TV that alternates between updating the solution ** x** and the pair (

The above figure compares reconstructions for two choices of the regularization parameter λ from simulation studies. From left to right:

- TV denotes a standard TV algorithm that ignores the uncertainty in the projection angles and uses the assumed angles
*θ*_{o}. These reconstructions are noisy due to the uncertain projection angles. - MV-TV denotes our new algorithm which produces improved reconstructions with less noise.
- TV-opt denotes TV reconstructions based on the true - but unknown - projection angles; this defines the unattainable optimal reconstruction from the noisy data.

For details see: N. A. B. Riis, Y. Dong, and P. C. Hansen, *Computed tomography reconstruction with uncertain view angles by iteratively updated model discrepancy*, submitted.