Bayesian Inference with Constraints Through Projection and Randomized Linear Least Squares Problems

In many applications, the signal that is being reconstructed must satisfy some constraints, e.g., material densities being positive, pixels of an image lying between 0 and 1, or the energy of a signal being bounded. Therefore, adding such prior information to a Bayesian model improves the meaningfulness of the resulting posterior distribution.

A common method to model constraints in a prior distribution is truncation, i.e., choosing an unconstrained prior distribution and remove all probability outside of the constraint set. Another method is to reparametrize the signal, e.g., if the signal x is assumed to be positive, then it can be written as x = exp(z) and one can choose an unconstrained prior distribution for z. These two methods focus on the interior of the constraint set, but very often, signals of interest lie on the boundary of the constraint set.

As natural alternative, one can first sample from an unconstrained posterior distribution and then apply a suitable projection onto a convex constraint set C. The resulting implicit posterior will consist of the original density in the interior of the constraint set, and a “lower dimensional” distribution on the boundary, see the illustration below.

If the forward model is linear and both the likelihood and prior are Gaussian distributions, then the unconstrained sampling and projection can be done in a single step by repeatedly solving randomized constrained linear least square problems of the form

min_x∈C || A x - b_s ||₂² + α || x - c_s ||₂² ,

where b_s and c_s are samples from suitable random distributions.  The solutions of these optimization problems will be samples of a projected posterior.

Consider the problem of deblurring a noisy signal with prior knowledge that the signal lies in the interval [0,1]. The figure above shows the pointwise median and 99% credibility interval for such a problem using no constraints (left), only the lower bound (middle) and both upper and lower bounds (right). The projection framework significantly improves the quality of the median near the boundary.

This example was created by PhD student Jasper Everink.

The Horseshoe Prior for Edge-Preserving Reconstruction

If we want to reconstruct an image with sharp edges, we can use a Cauchy or Laplace distribution of the difference between neighbor pixels as the prior (this is related to total variation regularization). Unfortunately, these priors are computationally demanding. As as a computationally attractive alternative, we can use a so-called horseshoe prior that resembles the Cauchy and Laplace priors - see the figure below where the horseshoe prior's density is guaranteed to lie in the shaded band.

The main advantage of the horseshoe prior is that it imposes a conditionally Gaussian distribution on the differences, which allows more efficient computations because Gaussians can be handled with efficient least-squares methods. The difficulty is that the horseshoe prior has hyperpriors with heavy tails, and to alleviate this issue we propose an extended horseshoe prior that uses a scale-mixture representation of the heavy-tailed hyperpriors.

For the posterior, we exploit the extended horseshoe prior to compute conditional distributions for the associated parameters that can be sampled in closed form due to conjugacy. This allows the application of a Gibbs sampler.

Here we illustrate our method with a 1D deconvolution problem for two different noise levels, and we compare with posterior statistics obtained with a Laplace-difference prior proposed by Uribe, Bardsley, Dong, Hansen & Riis (2022).  Compared to the Laplace-difference prior, we see that the horseshoe prior gives a sharper posterior (dotted lines) and a lower posterior uncertainty (the shaded area shows the 95% credibility interval). The solid line shows the ground truth.

The horseshoe prior originated in the statistical community, and here we introduce it in the setting of Bayesian inverse problem - see the paper by Uribe, Dong & Hansen (2022).

Boundary Detection in X-Ray CT Applications

When an X-ray passes through an object, its intensity is attenuated depending on the material's attenuation coefficient. The accumulated attenuation is the integral of the attenuation coefficient along the ray. We consider the case of 2D parallel-beam systems for which the measured data is called a sinogram.

Any X-ray going through the image is characterized by its slope θ with the horizontal axis. In practice, a full-angle geometry, i.e. when θ covers the entire interval [0,π), is not always possible. This could be due to obstacles in the CT measurement setup or to avoid unhealthy levels of radiation exposure. In such cases, a limited-angle geometry is needed, i.e., θ belongs to the reduced interval [0,θ_max) with θ_max ≪ π. In such a limited-angle setup, the quality of reconstructed image is compromised when conventional reconstruction methods, such as filtered back projection (FBP), are applied.

In many X-ray CT applications, the boundaries of objects contain valuable information, e.g., to differentiate between benign vs malignant tumors. A common approach in detecting boundaries is via an image segmentation post-processing step, but the quality of the estimated boundaries is compromised in a limited-angel imaging case.

We provide a novel goal-oriented, Bayesian framework for the limited-angel CT problem, in which we reconstruct the boundaries directly and quantify their uncertainty, and without the need for reconstructing the image. This approach avoids the error propagation and reduces the dimensionality of the problem from finding a 2D image to a 1D boundary of a region. In addition, our method estimates and performs UQ of the roughness/smoothness of the boundaries that carry important information for, e.g., medical imaging applications.

In the above plots, we see the performance of our method in estimating the boundaries of an object in the limited-angle CT setup with different θ_max. The arrows indicate the angles of the X-rays and thus region of the limited-angle interval [0, θ_max ). The estimated boundaries are shown in dark blue color, and the light blue region indicates the uncertainty in this estimation. The smaller the θ_max the larger the uncertainty region - which is consistent with the results from microlocal analysis.

A paper by Afkham, Dong, and Hansen is submitted to SIAM/ASA J. Uncertainty Quantification. A related paper is submitted to IEEE Trans. Image Processing.

Defect Detection in X-Ray CT of Subsea Pipes

X-ray computed tomography (CT) is used to monitor the condition of subsea oil or gas pipes in operation, in order to detect defects that might cause leaks. In this example, we use data obtained in the test facilities at FORCE Technology. Computations are done with our CUQIpy software which draws upon the Core Imaging Library (CIL) for the CT models.

We formulate a Bayesian inverse problem with built-in defect detection. The goal is to detect defects and quantify their uncertainties. We express the CT problem as y = A (x+d) + e, where y is the measured X-ray absorption, e is data noise, and A is the linear forward model representing the physics and geometry of the measurements. We use a novel representation x + d of the unknown image to be reconstructed, expressed as a sum of two images; x contains the pipe structure and d contains potential defects.

In the Bayesian setting we formulate the joint posterior distribution p(x,d | y) ∝ p(y | x,d) p(x) p(d), where p(x,d | y) is the joint posterior distribution representing the solution to the CT problem, p(y | x,d) is the likelihood that represents the data misfit, while p(x) and p(d) are prior distributions representing any knowledge we have about the unknown images.

We impose priors that promote the structures we are looking for. For the pipe structure x, we use the structural Gaussian prior proposed by Christensen, Riis, Uribe & Jørgensen (submitted). The above figure shows a sample from this prior which promotes the known layered pipe structure. We expect small and few defects, and therefore we impose a prior that promotes both sparsity and correlation (a gamma Markov random field) in the defect image d (as described in a paper by Christensen, Riis, Pereyra & Jørgensen). We explore the conditional posterior distributions p(x | y,d) and p(d | y,x) using Gibbs sampling.

The figure above shows the means of the x-samples (left), the d-samples (middle), and their sum with annotation (right). These results indicate that our methodology has successfully separated the defects from the overall pipe structure. In the figure below we zoom-in on the defect reconstruction d for the annotated defects indicated in the figure above, for further investigation of these defects. The figure shows the mean of the posterior samples, as well as related UQ in the form of the standard deviation of the samples. Engineers can use these results to identify critical defects in the subsea pipes.

This is joint work with Assoc. Prof. Marcelo Pereyra from Heriot-Watt University. The example was created by PhD student Silja L. Christensen.

Handling Model Errors in Inverse Problems

Using an ideal and complex mathematical model will produce a good reconstruction in an inverse problem, but at a high computational cost. Instead, one can use a simplified approximate model that requires less computational work, but this often results in a poorer reconstruction. We use the Bayesian framework to study the model error between the ideal and approximate models and find ways to reduce its influence when solving a Bayesian inverse problem.

We consider the general model y = F(x,η) + e , where y is the data, x is the reconstruction, F is the exact forward operator with model parameters η, and e is the noise. Given an approximate forward operator f and an initial estimate η₀ of η, we obtain the model y = f(x,η₀) + δ + e   with   δ = F(x,η) − f(x,η₀) , where δ denotes the model error. The corresponding likelihood then takes the form p(y | x) = p_v|x(y − f(x,η₀) | x)     with     v = δ + e .

Since p(v | x) is difficult to characterize, we may assume that δ and x are independent. Samples of δ can then be generated by "pushing forward" samples of x from the prior for x through the relation δ = F(x,η) − f(x,η0). Unfortunately, such samples of δ are inferior because they lack information related to the likelihood p. We can obtain better samples of δ by pushing forward samples of x obtained via the posterior — and this process can be repeated in an iterative fashion.

We develop methods based on the above ideas and apply them to reconstruction problems in magnetic resonance electrical impedance tomography (MREIT) where the solution x represents the electric conductivity inside a circular domain.

The plots below show that our method improves on results obtained by either neglecting the model error or representing it approximately. Specifically, we show the exact solution x (top left) and estimates of x obtained without the model error (top right), with the model error obtained by pushing forward samples from the prior (bottom left), and with the model error obtained by pushing forward samples from the posterior, repeated three times (bottom right).

This is joint work with Prof. Bangti Jin from Chinese University of Hong Kong. The example was created by PhD student Puyuan Mi.

Follow this LINK to see pictures from our team building event in the spring of 2023.

Specialist in convex optimization and numerical optimization algorithms.

Specialist in noise modeling and computational methods for image processing.

Link to portrait.

Specialist in numerical analysis, iterative solvers, and computational methods for inverse problems.

Link to portrait.

Specialist in methods and software for image reconstruction and their application in X-ray CT.

Specialist in functional analysis and physical modeling with uncertainty quantification.

Specialist in coupled-physics modeling, PDE-constrained optimization, functional analysis and electrical impedance tomography.

Link to portrait.

Project coordinator.

Project secretary.

Project: goal-oriented uncertainty quantification.

Project: computational UQ for PDE problems.

Project: machine learning and UQ for inverse problems.

Project: efficient and flexible computational methods.

Project: non-Gaussian priors.

Now with LUT University, Finland; homepage.

Project: Bayesian inference for inverse problems with sampling conditioned on functionals.

Project: prior modeling in computational UQ for inverse problems.

Project: UQ for tomographic reconstruction.

Project: computational UQ for inverse problems with implicit priors.

Project: dimensionality challenges in UQ for inverse problems.

Project: inverse problems with Besov prior.

Project: UQ for source localization and passive imaging in random media

Project: model error reduction in inverse problems. 

Project: computational UQ of hybrid inverse problems.

Professor Johnathan M. Bardsley, Department of Mathematical Sciences, University of Montana - specialist in computational methods for inverse problems and uncertainty quantification.

Associate Professor Julianne Chung, Department of Mathematics, Emory University - specialist in in computational methods for inverse problems and uncertainty quantification in imaging applications.

Associate Professor Matthias Chung, Department of Mathematics, Emory University - specialist in computational inverse problems, data analytics & learning, uncertainty quantification, and numerical optimization.

Professor Bangti Jin, Department of Mathematics, Chinese University of Hong Kong - specialist in theory and algorithms for inverse and ill-posed problems.

DTU Compute, Section for Mathematics - specialist in wavelets, frames, and signal analysis.

Associate Professor Marcelo Pereyra, Maxwell Institute for Mathematical Science and School of Mathematical and Computer Sciences, Heriot-Watts University - specialist in statistical, analytical and machine-learning paradigms.

Professor Tanja Tarvainen, Department of Applied Physics, University of Eastern Finland - specialist in computational methods for inverse problems and uncertainty quantification in imaging applications

Research Professor Faouzi Triki, Laboratoire Jean Kuntzman, Grenoble-Alpes University - specialist in inverse problems and mathematical modeling in optics.

Research Scientist Olivier Zahm, INRIA-Grenoble - specialist in model order reduction for uncertainty quantification.