Bayesian Inference for Inverse Problems

Abstract

Inverse problems arise everywhere we have indirect measurement. Regularization and Bayesian inference methods are two main approaches to handle inverse problems. Bayesian inference approach is more general and has much more tools for developing efficient methods for difficult problems. In this chapter, first, an overview of the Bayesian parameter estimation is presented, then we see the extension for inverse problems. The main difficulty is the great dimension of unknown quantity and the appropriate choice of the prior law. The second main difficulty is the computational aspects. Different approximate Bayesian computations and in particular the variational Bayesian approximation (VBA) methods are explained in details.

1. Introduction

Inverse problems arise in many scientific and engineering applications. In fact, almost always we want to infer on quantities, variables, distributions and functions which are not directly observable. Inferring on a hidden variable f via the observation of another variable g is the main objective in many scientific area [1, 2, 3].

Classically, the Bayesian methods have been developed for direct observation of a quantity, its parametric modeling and the estimation of the parameters of the model. Description and development of the Bayesian inference for the case of inverse problems is the main objective of this chapter. The chapter is organized as follows: First a few inverse problems are mentioned, mainly in two categories: those described by Ordinary Differential Equations (ODE) and Partial Differential Equations (PDE)’s and those described with integral equations. Then, we will see that two main problems arise: parameter estimation and inversion. First the Bayesian parameter estimation is described and then the Bayesian inversion.

2. Inverse problems

To see easily the two categories of inverse problems, first a very simple example is given. Consider the electrical circuit of Figure 1.

Figure 1.

A simple electrical circuit example to show two different expressions of inverse problems modeling: Ordinary differential equation (ODE) or Integral equation (IE).

Using the notations used on the figure, we can easily obtain the following ODE:

$$g\left(t\right)+\mathit{\theta }\frac{\partial g\left(t\right)}{\partial t}=f\left(t\right)$$E1

Then, using the Fourier transform (FT), we obtain easily the following integral equation:

$$f\left(t\right)=\int f\left(\tau \right)h\left(t-\tau \right)\mathrm{d}\tau$$E2

These two simple equations describe the same linear inverse problem, where we can distinguish the following mathematical problems:

• Forward problem: Given the parameter $\mathit{\theta }$ of the system and the input, f (t) predict the output g(t).

• Parameter estimation: Given the input f (t) and the output g(t), estimate the parameter $\mathit{\theta }.$

• System identification: Given the input f (t) and the output g(t), estimate the impulse response (IR) of the system h(t).

• Inverse problems:

  • Simple: Given the characteristics of the system (either the parameter $\mathit{\theta }$ or equivalently the impulse response h(t)) and the output g(t) estimate the input f (t);

  • Blind: Given the output g(t) estimate both the system, parameter $\mathit{\theta }$ or the impulse response h(t), and input f (t).

For general vocabulary and examples, see [2, 4, 5].

2.1 Examples of linear inverse problems

Here, a few examples of classical inverse problems are listed.

2.1.1 Deconvolution

When the forward problem is a convolution operation:

$$g\left(t\right)=\int f\left(\tau \right)h\left(t-\tau \right)\mathrm{d}\tau ,$$E3

the inverse problem is called Deconvolution (Figure 2). Figure 3 shows an example of deconvolution problem which arise in radio astronomy.

Figure 2.

Signal deconvolution problem.

Figure 3.

Image deconvolution or restoration inverse problem in sattelite imaging.

2.1.2 Image restoration

In many imaging systems, such as visual cameras, microscopes, telescopes or Infra Red cameras, due to some limitations such as limited aperture or limited resolution, the forward problem can be approximated by a 2D convolution equation:

$$g\left(x,y\right)=\iint f\left(x^{\prime },y^{\prime }\right)h\left(x-x^{\prime },y-y^{\prime }\right)\mathrm{d}x\mathrm{d}y.$$E4

The corresponding inverse problem is called image deconvolution or more often image restoration. The example given in Figure 4, is the case of satellite imaging [6, 7].

Figure 4.

Image reconstruction in CT. On the left, the projections $g\left(r,\varphi \right)$ and on the right the object f (x, y).

2.1.3 Image reconstruction in X ray computed tomography (CT)

In X-ray CT, assuming parallel geometry, where a ray is characterized by its angle $\varphi$ and its distance r from the center of the object f (x, y) the relation between the data $g\left(r,\varphi \right),$ called projections at angle $\varphi$ and the function f (x, y), called object, is given by the Radon transform:

$$g\left(r,\varphi \right)=\mathit{\iint }f\left(x,y\right)\delta \left(r-xcos\varphi -ysin\varphi \right)\mathrm{d}x\mathrm{d}y.$$E5

The inverse problem here is called Image reconstruction. A simulated example is shown in Figure 4.

2.1.4 Fourier synthesis

In many imaging systems, when using FT, it is possible to model the inverse problem with the following forward FT relation:

$$g\left(u,v\right)=\iint f\left(x,y\right)exp\left[-j\left(\mathit{ux}+\mathit{vy}\right)\right]\mathrm{d}x\mathrm{d}y,$$E6

where the data, after an appropriate FT, can fill partially the Fourier domain g(u, v) of the unknown interested function f (x, y) [8, 9]. Figure 5 shows the case of X-ray CT.

Figure 5.

Fourier synthesis (FS) inverse problems arising in many imaging systems. Here is illustrated the FS problem in X-ray CT.

2.1.5 General linear inverse problems

All the examples of the linear inverse problems listed above, can be summarized in the following general form:

$$g\left(\mathit{s}\right)=\int f\left(\mathit{r}\right)h\left(\mathit{s},\mathit{r}\right)\mathrm{d}\mathit{r}$$E7

where s can be either $t,\left(x,y\right),\left(r,\varphi \right)$ or (u, v) and r, respectively $\tau ,\left(x^{\prime },y^{\prime }\right),\left(x,y\right)$ and (x, y).

3. Bayesian parameter estimation

To introduce, in a very simple way, the Bayes rule for parameter estimation, we consider the case where we have a set of data: $\mathit{g}=\left\{g_1,\cdots ,g_n\right\}$ where we assign them a probability law $p\left(g_i|\mathit{\theta }\right)$ with a set of unknown parameters $\mathit{\theta }.$ The question now is how to infer $\mathit{\theta }$ from those data. We can immediately use the Bayes rule:

$$p\left(\mathit{\theta }|\mathit{g}\right)=\frac{p\left(\mathit{g}|\mathit{\theta }\right)p\left(\mathit{\theta }\right)}{p\left(\mathit{g}\right)}\propto l\left(\mathit{\theta }\right)p\left(\mathit{\theta }\right)$$E8

where:

• $l\left(\mathit{\theta }\right)\triangleq p\left(g|\mathit{\theta }\right)={\mathrm{\Pi }}_{i}p\left(g_i|\mathit{\theta }\right)$ is called the likelihood, representing the uncertainty in the data knowing the parameters;

• $p\left(\mathit{\theta }\right)$ is called the prior or a priori, a probability law assigned to the parameters to represent the prior knowledge (to the observation data) we may have on those parameters;

• the denominator p(g)

$$p\left(\mathit{g}\right)=\int p\left(g|\mathit{\theta }\right)p\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }$$E9

is called the evidence.

So, the process of using the Bayes rule for parameter estimation can be summarized as follows:

• Write the expression of the likelihood $p\left(\mathit{g}|\mathit{\theta }\right)$

• Assign the prior $p\left(\mathit{\theta }\right)$ to translate all we know about $\mathit{\theta }$ before observing the data g

• Apply the Bayes rule to obtain the expression of the posterior law $p\left(\mathit{\theta }|\mathit{g}\right)$

• Use the posterior $p\left(\mathit{\theta }|\mathit{g}\right)$ to do any inference on $\mathit{\theta }.$ For example:

  • Compute its expected value, called Expected A Posteriori (EAP) or Posterior Mean (PM):

    $${\widehat{\mathit{\theta }}}_{\mathit{PM}}=\int \mathit{\theta p}\left(\mathit{\theta }|\mathit{g}\right)\mathrm{d}\mathit{\theta }$$E10

  • Compute the value of $\mathit{\theta }$ for which the $p\left(\mathit{\theta }|\mathit{g}\right)$ is maximum; Maximum A Posteriori (MAP):

    $${\widehat{\mathit{\theta }}}_{\mathit{MAP}}=\underset{\mathit{\theta }}{argmax}\left\{p\left(\mathit{\theta }|\mathit{g}\right)\right\}$$E11

  • Sampling and exploring [Monte Carlo methods]

$$\mathit{\theta }\sim p\left(\mathit{\theta }|\mathit{g}\right)$$

which gives the possibility to obtain any statistical information we want to know about $\mathit{\theta }.$ For example, if we generate N samples $\left\{{\mathit{\theta }}_1,\cdots ,{\mathit{\theta }}_N\right\},$ for large enough N, we have:

$$\mathrm{E}\left\{\mathit{\theta }\right\}\simeq \frac{1}{N}\sum _{n=1}^N{\mathit{\theta }}_n.$$E12

3.1 One parameter case

When $\mathit{\theta }$ is a scalar quantity, then, we can also do the following computations:

• Compute the value of ${\mathit{\theta }}_{\mathit{Med}}$ such that:

$$P\left(\mathit{\theta }>{\mathit{\theta }}_{\mathit{Med}}\right)=P\left(\mathit{\theta }<{\mathit{\theta }}_{\mathit{Med}}\right)$$E13

which is called the median value. Its computation needs integration:

$${\int }_{-\infty }^{{\mathit{\theta }}_{\mathit{Med}}}p\left(\mathit{\theta }|g\right)\mathrm{d}\mathit{\theta }={\int }_{{\mathit{\theta }}_{\mathit{Med}}}^{\infty }p\left(\mathit{\theta }|\mathit{g}\right)\mathrm{d}\mathit{\theta }$$E14

• Compute the value ${\mathit{\theta }}_{\alpha },$ called $\alpha$ quantile, for which

$$P\left(\mathit{\theta }>{\mathit{\theta }}_{\alpha }\right)=1-P\left(\mathit{\theta }<{\mathit{\theta }}_{\alpha }\right)={\int }_{{\widehat{\mathit{\theta }}}_{\alpha }}^{\infty }p\left(\mathit{\theta }|\mathit{g}\right)\mathrm{d}\mathit{\theta }=1-\alpha$$E15

• Region of high probabilities: [needs integration methods]

$$\left[{\widehat{\mathit{\theta }}}_1,{\widehat{\mathit{\theta }}}_2\right]:{\int }_{{\widehat{\mathit{\theta }}}_1}^{{\widehat{\mathit{\theta }}}_2}p\left(\mathit{\theta }|\mathit{g}\right)\mathrm{d}\mathit{\theta }=1-\alpha$$

Bayes rule and Bayesian estimation can be illustrated as follows:

Two main points are of great importance:

• How to assign the prior $p\left(\mathit{\theta }\right)$ in the second step; and

• How to do the computations in the last step.

This last problem becomes more serious with multi parameter case.

3.2 Multi-parameter case

If we have more than one parameter, then $\mathit{\theta }={\left[{\mathit{\theta }}_1,\cdots ,{\mathit{\theta }}_n\right]}^{\prime }.$ The Bayes rule still holds:

$$p\left(\mathit{\theta }|g\right)=\frac{p\left(\mathit{g}|\mathit{\theta }\right)p\left(\mathit{\theta }|\mathit{g}\right)}{p\left(\mathit{g}\right)}$$E16

Now, again, we can compute:

• The Expected A Posteriori (EAP):

$${\widehat{\mathit{\theta }}}_{\mathit{PM}}=\mathit{\int }\mathit{\theta }p\left(\mathit{\theta }|\mathit{g}\right)\mathrm{d}\mathit{\theta },$$E17

but this needs efficient integration methods.

• The Maximum A Posteriori (MAP):

$${\widehat{\mathit{\theta }}}_{\mathit{MAP}}=\underset{\mathit{\theta }}{argmax}\left\{p\left(\mathit{\theta }|\mathit{g}\right)\right\}$$E18

but this needs efficient optimization methods.

• Sampling and exploring [Monte Carlo methods]

$$\mathit{\theta }\sim p\left(\mathit{\theta }|\mathit{g}\right)$$

but this needs efficient sampling methods.

• We may also try to localize the region of the highest probability:

$$P\left(\mathit{\theta }\in \mathrm{\Theta }\right)={\mathit{\int }}_{\mathrm{\Theta }}p\left(\mathit{\theta }|g\right)\mathrm{d}\mathit{\theta }=1-\alpha$$E19

for a given small $\alpha ,$ but this problem may not have a unique solution.

4. Bayesian inference for inverse problems

As described before, in inverse problems, the unknown f is a function (of time, space, wavelength, …) and the observable quantity g is also another function which is related to f via an operator $g=\mathcal{H}\left(f\right)+ϵ.$ When discretized, they can be represented by the great dimensional vectors f, g and $\mathit{g}=\mathit{H}\left(\mathit{f}\right)+\mathit{ϵ},$ where $ϵ$ represents all the errors (measurement, model and discretization). When the operator is a linear one, we have:

$$\mathit{g}=\mathit{H}\mathit{f}+\mathit{ϵ},$$E20

where f is a vector of length n, H the known forward model matrix of size $m\times n,$ and g and $\mathit{ϵ}$ two vectors of size m.

The Bayes rule for this case is written as:

$$p\left(\mathit{f}|g,\mathcal{M}\right)=\frac{p\left(\mathit{g}|\mathit{f},\mathcal{M}\right)p\left(\mathit{f}|,\mathcal{M}\right)}{p\left(\mathit{g},\mathcal{M}\right)}$$E21

where we introduce $\mathcal{M}$ to represent the model, $p\left(\mathit{g}|\mathit{f},\mathcal{M}\right),$ called commonly the likelihood, is obtained using the forward model (20) and the assigned probability law of the noise $p\left(\mathit{ϵ}\right),p\left(\mathit{f}|,\mathcal{M}\right)$ is the assigned prior model and $p\left(\mathit{f}|\mathit{g},\mathcal{M}\right)$ the posterior probability law. Figure 6 shows in a schematic way the main ingredients of the Bayesian inference for inverse problems.

Figure 6.

Illustration of the Bayesian inference for inverse problems.

This even very simple linear model has been used in many areas: linear inverse problems, compressed sensing, curve fitting and linear regression, machine learning, etc.

In inverse problems such as deconvolution, image restoration, f represent the input or original image, g represents the blurred and noisy image and H is the convolution operator matrix. In image reconstruction in Computed Tomography (CT), f represents the distribution of some internal property of an object, for example the density of the material and g represents the radiography data and H is the radiographic projection operator (discretized Radon transform operator).

In Compress Sensing, g is the compressed data, f is the uncompressed image and H the compressing matrix. In machine learning, g are the data, H is a dictionary and f represents the sparse coefficients of the projections of the data on that dictionary.

5. Hyperparameter estimation

When applying the Bayes rule, the main terms which are the likelihood and prior depend on parameters, which cannot be fixed in practical situation. We may thus want to estimate them from the data. In the Bayesian approach, this can be done easily:

$$p\left(\mathit{f},{\mathit{\theta }}_1,{\mathit{\theta }}_2|\mathit{g}\right)=\frac{p\left(\mathit{g}|\mathit{f},{\mathit{\theta }}_1\right)p\left(\mathit{f},{\mathit{\theta }}_2\right)p\left({\mathit{\theta }}_1\right)p\left({\mathit{\theta }}_2\right)}{p\left(\mathit{g}\right)}$$E22

where $p\left({\mathit{\theta }}_1\right)$ and $p\left({\mathit{\theta }}_2\right)$ are the prior probability laws assigned to ${\mathit{\theta }}_1$ and ${\mathit{\theta }}_2$ and often $p\left(\mathit{\theta }\right)=p\left({\mathit{\theta }}_1\right)p\left({\mathit{\theta }}_2\right).$ We can then write more succinctly:

$$p\left(\mathit{f},\mathit{\theta }|\mathit{g},{\mathit{\theta }}_0\right)=\frac{p\left(\mathit{g}|\mathit{f},{\mathit{\theta }}_1\right)p\left(\mathit{f},{\mathit{\theta }}_2\right)p\left(\mathit{\theta }\right)}{p\left(\mathit{g}\right)}$$E23

The scheme of this situation is illustrated in Figure 7.

Figure 7.

Illustration of the Bayesian approach for inverse problems with unknown hyperparameters.

From here, we have different directions for doing estimation:

5.1 Joint maximum a posteriori (JMAP)

Rewriting the expression of the joint posterior law:

$$p\left(\mathit{f},\mathit{\theta }|\mathit{g}\right)=\frac{p\left(\mathit{g}|\mathit{f},{\mathit{\theta }}_1\right)p\left(\mathit{f}|{\mathit{\theta }}_2\right)p\left(\mathit{\theta }\right)}{p\left(\mathit{g}\right)}\propto p\left(\mathit{g}|\mathit{f},{\mathit{\theta }}_1\right)p\left(\mathit{f},{\mathit{\theta }}_2\right)p\left(\mathit{\theta }\right)$$E24

where $\propto$ means equal up to a constant factor which is $1/p\left(\mathit{g}\right).$ In this case, we can try to optimize it with respect to its two arguments:

$$\left(\widehat{\mathit{f}},\widehat{\mathit{\theta }}\right)=\underset{\left(\mathit{f},\mathit{\theta }\right)}{argmax}\left\{p\left(\mathit{f},\widehat{\mathit{\theta }}|\mathit{g}\right)\right\}$$E25

This can be done, for example, by alternate optimization:

$$\{\begin{array}{l}{\widehat{\mathit{f}}}^{\left(k+1\right)}=arg{max}_{\mathit{f}}\left\{p\left(\mathit{f},{\widehat{\mathit{\theta }}}^{\left(k\right)}|\mathit{g}\right)\right\}\\ {\widehat{\mathit{\theta }}}^{\left(k+1\right)}=arg{max}_{\mathit{\theta }}\left\{p\left({\mathit{f}}^{\left(k\right)},\widehat{\mathit{\theta }}|\mathit{g}\right)\right\}\end{array}$$E26

When the optimization algorithm is successful, we have the optimal values of $\widehat{\mathit{f}}$ and $\widehat{\mathit{\theta }}.$ This method can be summarized as follows:

5.2 Marginalization over $\mathit{\theta }$

The main idea here is to consider $\mathit{\theta }$ as a nuisance parameter. Thus, integrating it out, we get

$$p\left(\mathit{f}|\mathit{g}\right)=\int p\left(\mathit{f},\mathit{\theta }|\mathit{g}\right)\mathrm{d}\mathit{\theta }$$E27

which can be used to infer on f. Also, if we still want to get estimates of $\mathit{\theta },$ we can first obtain an estimate $\widehat{\mathit{f}}$ for f and then, if needed, to use it as it is illustrated in the following scheme:

5.3 Marginalization over f

The main idea here is first find a good estimate for the parameters $\mathit{\theta }$ and then use it for the inference on f. So, first obtain:

$$p\left(\mathit{\theta }|\mathit{g}\right)=\int p\left(\mathit{f},\mathit{\theta }|\mathit{g}\right)\mathrm{d}\mathit{f}$$E28

which can be used to first estimate $\mathit{\theta }$ and then use it. For example, the method which is related to the Second type Maximum likelihood, first estimate $\widehat{\mathit{\theta }}$ by

$$\widehat{\mathit{\theta }}=\underset{\mathit{\theta }}{argmax}\left\{p\left(\widehat{\mathit{\theta }}|\mathit{g}\right)\right\}$$E29

and then use it with $p\left(\mathit{f}|\widehat{\mathit{\theta }},\mathit{g}\right)$ to infer on $\widehat{\mathit{f}}.$ For a flat prior model, $p\left(\mathit{\theta }|g\right)\propto p\left(g|\mathit{\theta }\right)$ which is called the likelihood and the estimator

$$\widehat{\mathit{\theta }}=\underset{\mathit{\theta }}{argmax}\left\{p\left(\widehat{\mathit{\theta }}|g\right)\right\}=\underset{\mathit{\theta }}{argmax}\left\{p\left(g|\widehat{\mathit{\theta }}\right)\right\}$$E30

is called Maximum Likelihood (ML) and the whole approach is called ML of second type. This method can be summarized as follows:

The main difficulty in this approach is that, rarely we can have an analytical expression for the first marginalization. To overcome this difficulty, many algorithms have been proposed to compute f. One of them is called Expectation- Maximization (EM) and its generalization (GEM). The main idea of these algorithms are summarized in the following subsections:

5.3.1 EM and GEM algorithms

To summarize these methods, we use the vocabulary of the main authors of EM method, where f is considered as hidden variable, g as incomplete data, (g, f) as complete data, $lnp\left(\mathit{g}|\mathit{\theta }\right)$ incomplete data log-likelihood and $lnp\left(\mathit{g},\mathit{f}|\mathit{\theta }\right)$ as complete data log-likelihood. Then, the following iterative algorithms describe the EM and GEM algorithms.:

• EM Iterative algorithm:

$$\{\begin{array}{ll}\mathrm{\text{E‐step:}}& Q\left(\mathit{\theta },{\widehat{\mathit{\theta }}}^{\left(k\right)}\right)={\mathrm{E}}_{p\left(\mathit{f}|\mathit{g},{\widehat{\mathit{\theta }}}^{\left(k\right)}\right)}\left\{lnp\left(\mathit{g},\mathit{f}|\mathit{\theta }\right)\right\}\\ \\ \mathrm{\text{M‐step:}}& {\widehat{\mathit{\theta }}}^{\left(k\right)}=arg{max}_{\mathit{\theta }}\left\{Q\left(\mathit{\theta },{\widehat{\mathit{\theta }}}^{\left(k-1\right)}\right)\right\}\end{array}$$E31

• GEM (Bayesian) algorithm:

$$\{\begin{array}{ll}\mathrm{\text{E‐step:}}& Q\left(\mathit{\theta },{\widehat{\mathit{\theta }}}^{\left(k\right)}\right)={\mathrm{E}}_{p\left(\mathit{f}|\mathit{g},{\widehat{\mathit{\theta }}}^{\left(k\right)}\right)}\left\{lnp\left(\mathit{g},\mathit{f}|\mathit{\theta }\right)+lnp\left(\mathit{\theta }\right)\right\}\\ \\ \mathrm{\text{M‐step:}}& {\widehat{\mathit{\theta }}}^{\left(k\right)}=arg{max}_{\mathit{\theta }}\left\{Q\left(\mathit{\theta },{\widehat{\mathit{\theta }}}^{\left(k-1\right)}\right)\right\}\end{array}$$E32

These methods can be summarized in the following scheme:

$$\boxed{p\left(\mathit{f},\mathit{\theta }|\mathit{g}\right)}\to \boxed{\mathrm{EM},\mathrm{GEM}}\to \widehat{\mathit{\theta }}\to \boxed{p\left(\mathit{f}|\widehat{\mathit{\theta }},\mathit{g}\right)}\to \widehat{\mathit{f}}$$

5.4 Variational Bayesian approximation

VBA is a powerful approach to do approximate Bayesian computation. It starts by first obtaining the expression of the joint $p\left(f,\mathit{\theta }|g\right)$ and then by approximating it with a simpler probability law $q\left(\mathit{f},\mathit{\theta }|\mathit{g}\right)$ which can be handled much easily for the computations. VBA can be summarized in the following steps:

• Approximate $p\left(\mathit{f},\mathit{\theta }|\mathit{g}\right)$ by $q\left(\mathit{f},\mathit{\theta }|\mathit{g}\right)=q_1\left(\mathit{f}|\mathit{g}\right)q_2\left(\mathit{\theta }|\mathit{g}\right)$ and then continue computations.

• To do this approximation, we need a criterion to qualify the approximation. The standard criterion to measure the proximity of two probability laws p and q is the Kullback–Leibler (KL) criterion $\mathrm{KL}\left(q\left(\mathit{f},\mathit{\theta }|\mathit{g}\right):p\left(\mathit{f},\mathit{\theta }|\mathit{g}\right)\right).$

• It is easy to show that:

$$\begin{array}{l}\mathrm{KL}\left(q:p\right)=\iint qlnq/p=\iint q_1{q}_{2}ln\frac{q_1{q}_2}{p}\\ =\int q_{1}ln{q}_1+\int q_{2}ln{q}_2-\iint qlnp\\ =-H\left(q_1\right)-H\left(q_2\right)-<lnp{>}_q\end{array}$$E33

• Alternate optimization of KL(q₁q₂:p) with respect to q₁ and q₂ results to:

$$\{\begin{array}{cc}{q}_1\left(\mathit{f}\right)& \propto exp\left[{⟨lnp\left(\mathit{g},\mathit{f},\mathit{\theta };\mathcal{M}\right)⟩}_{q_2\left(\mathit{\theta }\right)}\right]\\ q_2\left(\mathit{\theta }\right)& \propto exp\left[{⟨lnp\left(\mathit{g},\mathit{f},\mathit{\theta };\mathcal{M}\right)⟩}_{q_1\left(\mathit{f}\right)}\right]\end{array}$$E34

As KL(q₁q₂:p) is convex as a function of q₁ and q₂, the algorithm converges (locally) to the optimum solution. At the end, we have the expressions of q₁(f) and $q_2\left(\mathit{\theta }\right)$ which can, then, be used to infer on f and $\mathit{\theta }.$ VBA is summarized in the following scheme:

$$\boxed{p\left(\mathit{f},\mathit{\theta }|\mathit{g}\right)}\to \boxed{\begin{array}{c}\text{Variation}\\ \text{Bayesian}\\ \text{Approximation}\end{array}}\begin{array}{l}\to q_1\left(\mathit{f}\right)\to \widehat{\mathit{f}}\\ \to q_2\left(\mathit{\theta }\right)\to \widehat{\mathit{\theta }}\end{array}$$

In real applications, we choose parametric probability law for q₁(f) $q_2\left(\mathit{\theta }\right),$ and so, the iterations will be done on the parameters. What is interesting is that, choosing appropriate parametric models for q₁(f) and $q_2\left(\mathit{\theta }\right)$ we obtain either JMAP and GEM as special cases.

• Case 1: Deterministic or degenerate expressions $\to$ Joint MAP

$$\{\begin{array}{ll}{\widehat{q}}_1\left(\mathit{f}|\stackrel{\sim }{\mathit{f}}\right)& =\delta \left(\mathit{f}-\stackrel{\sim }{\mathit{f}}\right)\\ {\widehat{q}}_2\left(\mathit{\theta }|\stackrel{\sim }{\mathit{\theta }}\right)& =\delta \left(\mathit{\theta }-\stackrel{\sim }{\mathit{\theta }}\right)\end{array}\to \{\begin{array}{l}\stackrel{\sim }{\mathit{f}}=arg{max}_{\mathit{f}}\left\{p\left(\mathit{f},\stackrel{\sim }{\mathit{\theta }}|\mathit{g};\mathcal{M}\right)\right\}\\ \stackrel{\sim }{\mathit{\theta }}=arg{max}_{\mathit{\theta }}\left\{p\left(\stackrel{\sim }{\mathit{f}},\mathit{\theta }|\mathit{g};\mathcal{M}\right)\right\}\end{array}$$E35

• Case 2: Degenerate expression for $\mathit{\theta }$ and marginal expression for $\mathit{f}\to \mathrm{EM}$

$$\{\begin{array}{l}{\widehat{q}}_1\left(\mathit{f}\right)\propto p\left(\mathit{f}|\mathit{\theta },\mathit{g}\right)\\ {\widehat{q}}_2\left(\mathit{\theta }|\stackrel{\sim }{\mathit{\theta }}\right)=\delta \left(\mathit{\theta }-\stackrel{\sim }{\mathit{\theta }}\right)\end{array}\to \{\begin{array}{ll}Q\left(\mathit{\theta },\stackrel{\sim }{\mathit{\theta }}\right)& ={⟨lnp\left(\mathit{f},\mathit{\theta }|\mathit{g};\mathcal{M}\right)⟩}_{q_1\left(\mathit{f}|\stackrel{\sim }{\mathit{\theta }}\right)}\\ \stackrel{\sim }{\mathit{\theta }}& =arg{max}_{\mathit{\theta }}\left\{Q\left(\mathit{\theta }|\stackrel{\sim }{\mathit{\theta }}\right)\right\}\end{array}$$E36

• Case 3: q₁ and q₂ are chosen proportional to the marginals $p\left(\mathit{f}|\stackrel{\sim }{\mathit{\theta }},\mathit{g};\mathcal{M}\right)$ and $p\left(\mathit{\theta }|\stackrel{\sim }{\mathit{f}},\mathit{g};\mathcal{M}\right).$ This is a very appropriate choice for inverse problems, in particular cases where we use the exponential families and conjugate priors.

$$\{\begin{array}{ll}{\widehat{q}}_1\left(\mathit{f}\right)& \propto p\left(\mathit{f}|\stackrel{\sim }{\mathit{\theta }},\mathit{g};\mathcal{M}\right)\\ {\widehat{q}}_2\left(\mathit{\theta }\right)& \propto p\left(\mathit{\theta }|\stackrel{\sim }{\mathit{f}},\mathit{g};\mathcal{M}\right)\end{array}\to \{\begin{array}{l}\text{Accounts for the uncertainties of}\\ \stackrel{\sim }{\mathit{\theta }}\text{for}\widehat{\mathit{f}}\text{and vise versa}.\end{array}$$E37

In the following schemes these three cases are illustrated for comparison.

• JMAP Alternate optimization Algorithm:

$$\begin{array}{l}{\mathit{\theta }}^{\left(0\right)}\to \stackrel{\sim }{\mathit{\theta }}\to \boxed{\stackrel{\sim }{\mathit{f}}=arg{max}_{\mathit{f}}\left\{p\left(\mathit{f},\stackrel{\sim }{\mathit{\theta }}|\mathit{g}\right)\right\}}\to \stackrel{\sim }{\mathit{f}}\to \widehat{\mathit{f}}\\ \uparrow \downarrow \\ \widehat{\mathit{\theta }}\leftarrow \stackrel{\sim }{\mathit{\theta }}\leftarrow \boxed{\stackrel{\sim }{\mathit{\theta }}=arg{max}_{\mathit{\theta }}\left\{p\left(\stackrel{\sim }{\mathit{f}},\mathit{\theta }|\mathit{g}\right)\right\}}\leftarrow \stackrel{\sim }{\mathit{f}}\end{array}$$

• EM:

$$\begin{array}{l}{\mathit{\theta }}^{\left(0\right)}\to \stackrel{\sim }{\mathit{\theta }}\to \boxed{q_1\left(\mathit{f}\right)=p\left(\mathit{f}|\stackrel{\sim }{\mathit{\theta }},\mathit{g}\right)}\to q_1\left(\mathit{f}\right)\to \widehat{\mathit{f}}\\ \uparrow \downarrow \\ \widehat{\mathit{\theta }}\leftarrow \stackrel{\sim }{\mathit{\theta }}\leftarrow \boxed{\begin{array}{c}Q\left(\mathit{\theta },\stackrel{\sim }{\mathit{\theta }}\right)=⟨\ln p(\mathit{f},\mathit{\theta }\mid \mathit{g})⟩{}_{q_1\left(\mathit{f}\right)}\\ \stackrel{\sim }{\mathit{\theta }}=arg{max}_{\mathit{\theta }}\left\{Q\left(\mathit{\theta },\stackrel{\sim }{\mathit{\theta }}\right)\right\}\end{array}}\leftarrow q_1\left(\mathit{f}\right)\end{array}$$

• VBA:

$$\begin{array}{c}{\mathit{\theta }}^{\left(0\right)}\to q_2\left(\mathit{\theta }\right)\to \boxed{q_1\left(\mathit{f}\right)\propto exp\left[{⟨lnp\left(\mathit{f},\mathit{\theta }|\mathit{g}\right)⟩}_{q_2\left(\mathit{\theta }\right)}\right]}\to q_1\left(\mathit{f}\right)\to \widehat{\mathit{f}}\\ \uparrow \downarrow \\ \widehat{\mathit{\theta }}\leftarrow q_2\left(\mathit{\theta }\right)\leftarrow \boxed{q_2\mathit{\theta }\propto exp\left[{⟨lnp\left(\mathit{f},\mathit{\theta }|\mathit{g}\right)⟩}_{q_1\left(\mathit{f}\right)}\right]}\leftarrow q_1\left(\mathit{f}\right)\end{array}$$

5.5 Hierarchical priors

One last extension is the case where f, itself depends on another hidden variable z. So that we have:

$$p\left(\mathit{f},\mathit{z},\mathit{\theta }|\mathit{g}\right)\propto p\left(\mathit{g}|\mathit{f},{\mathit{\theta }}_1\right)p\left(\mathit{f}|\mathit{z},{\mathit{\theta }}_2\right)p\left(\mathit{z}|{\mathit{\theta }}_3\right)p\left(\mathit{\theta }\right),$$E38

where $\mathit{\theta }=\left({\mathit{\theta }}_1,{\mathit{\theta }}_2,{\mathit{\theta }}_3\right).$ This situation is shown in Figure 8. Again, here, we may only be interested to f or (f, z) or to all the three variables (f, z, $\mathit{\theta }$). Here too, we can either use methods of JMAP, marginalization or VBA to infer on these unknowns.

Figure 8.

Illustration of advanced Bayesian approach with hierarchical prior modeling with hidden variables.

6. Linear forward models and Gaussian case

Linear models are of importance. Gaussian prior laws are the most common and the easiest ones to handle. Also, many non-linear problems can be approximated by equivalent linear ones. Linear models with Gaussian prior laws are the easiest and powerful tools for a great number of scientific problems. In this section, an overview and some main properties are given.

6.1 Simple supervised case

Let consider the linear forward model we considered in previous section

$$\mathit{g}=\mathit{H}\mathit{f}+\mathit{ϵ},$$E39

and assign Gaussian laws to $\mathit{ϵ}$ and f which leads to:

$$\{\begin{array}{l}p\left(\mathit{g}|\mathit{f}\right)=\mathcal{N}\left(\mathit{g}|\mathit{H}\mathit{f},v_{ϵ}\mathit{I}\right)\propto exp\left[-\frac{1}{2v_{ϵ}}{‖\mathit{g}-\mathit{H}\mathit{f}‖}^2\right]\\ p\left(\mathit{f}\right)=\mathcal{N}\left(\mathit{f}|0,v_{\mathit{f}}\mathit{I}\right)\propto exp\left[-\frac{1}{2v_{\mathit{f}}}{‖\mathit{f}‖}^2\right]\end{array}$$E40

Using these expressions, we get:

$$\{\begin{array}{l}p\left(\mathit{f}|\mathit{g}\right)\propto exp\left[-\frac{1}{2v_{ϵ}}{‖\mathit{g}-\mathit{H}\mathit{f}‖}^2-\frac{1}{2v_{\mathit{f}}}{‖\mathit{f}‖}^2\right]\\ \propto exp\left[-\frac{1}{2v_{ϵ}}J\left(\mathit{f}\right)\right]\text{with}J\left(\mathit{f}\right)={‖\mathit{g}-\mathit{H}\mathit{f}‖}^2+\lambda {‖\mathit{f}‖}^2,\lambda =\frac{v_{ϵ}}{v_{\mathit{f}}}\end{array}$$E41

which can be summarized as:

$$p\left(\mathit{f}|\mathit{g}\right)=\mathcal{N}\left(\mathit{f}|\widehat{\mathit{f}},\widehat{\Sigma }\right)\text{with}\widehat{\mathit{f}}={\left[{\mathit{H}}^{\prime }\mathit{H}+\mathit{\lambda I}\right]}^{-1}{\mathit{H}}^{\prime }\mathit{g}\text{and}\widehat{\Sigma }=v_{ϵ}{\left[{\mathit{H}}^{\prime }\mathit{H}+\mathit{\lambda I}\right]}^{-1},$$E42

where $\lambda =\frac{v_{ϵ}}{v_f}.$ This case is summarized in Figure 9.

Figure 9.

Supervised linear Gaussian case.

This is the simplest case where we know exactly the expression of the posterior law and all the computations can be done explicitly. However, for great dimensional problems, where the vectors f and g are very great dimensional, we may even not be able to keep in memory the matrix H and surely not be able to compute the inverse of the matrix $\left[{\mathit{H}}^{\prime }\mathit{H}+\mathit{\lambda }\mathit{I}\right].$ In Section 9 on Bayesian computation, We will see how to do these computations.

6.2 Unsupervised case or hyperparameter estimation

In the previous section, we considered the linear models with Gaussian priors with known parameters $v_{ϵ}$ and v_f. In many practical situations these parameters are not known, and we want to estimate them too. For this, we can assign them too prior laws. As the variances are positive quantities and using the concept of conjugate priors, we can assign then Inverse Gamma priors:

$$\{\begin{array}{l}p\left(v_{ϵ}\right)=\mathcal{IG}\left(v_f|{\alpha }_{{ϵ}_0},{\beta }_{{ϵ}_0}\right)\\ p\left(v_f\right)=\mathcal{IG}\left(v_f|{\alpha }_{f_0},{\beta }_{f_0}\right)\end{array}$$E43

and using the likelihood $p\left(\mathit{g}|\mathit{f},v_{ϵ}\right)=\mathcal{N}\left(\mathit{g}|\mathit{H}\mathit{f},v_{ϵ}\mathit{I}\right)$ and the prior $p\left(\mathit{f}|v_{\mathit{f}}\right)=\mathcal{N}\left(\mathit{f}|0,v_{\mathit{f}}\mathit{I}\right),$ we can easily obtain the expressions of the following conditional posterior laws:

$$\{\begin{array}{l}p\left(\mathit{f}|\mathit{g},{\widehat{v}}_{ϵ},{\widehat{v}}_{\mathit{f}}\right)=\mathcal{N}\left(\mathit{f}|\widehat{\mathit{f}},\widehat{\Sigma }\right)\text{with}:\\ \widehat{\mathit{f}}={\left[{\mathit{H}}^t\mathit{H}+\widehat{\lambda }\mathit{I}\right]}^{-1}{\mathit{H}}^t\mathit{g}\\ \widehat{\Sigma }={\widehat{v}}_{ϵ}{\left[{\mathit{H}}^t\mathit{H}+\widehat{\lambda }\mathit{I}\right]}^{-1},\widehat{\lambda }={\widehat{v}}_{ϵ}/{\widehat{v}}_{\mathit{f}}\end{array}$$E44

and

$$\{\begin{array}{l}p\left(v_{ϵ}|\mathit{g},\mathit{f}\right)=\mathcal{IG}\left(v_{ϵ}|{\stackrel{\sim }{\alpha }}_{ϵ},\stackrel{\sim }{\beta }ϵ\right)\\ p\left(v_{\mathit{f}}|\mathit{g},\mathit{f}\right)=\mathcal{IG}\left(v_{\mathit{f}}|{\stackrel{\sim }{\alpha }}_{\mathit{f}},{\stackrel{\sim }{\beta }}_{\mathit{f}}\right)\end{array}$$E45

where all the details and in particular the expressions for ${\stackrel{\sim }{\alpha }}_{ϵ},{\stackrel{\sim }{\beta }}_{ϵ},{\stackrel{\sim }{\alpha }}_f,{\stackrel{\sim }{\beta }}_f$ can be found in [10].

As we can see, the expressions of $\widehat{\mathit{f}}$ and $\widehat{\mathbf{\Sigma }}$ are the same as in the previous case, except that the values of ${\widehat{v}}_{ϵ},$${\widehat{v}}_{\mathit{f}}$ and $\widehat{\lambda }$ have to be updated. They are obtained from the conditionals $p\left(v_{ϵ}|\mathit{g},\mathit{f}\right)$ and $p\left(v_{\mathit{f}}|\mathit{g},\mathit{f}\right)$ which depend on f. This shows that we can propose an iterative algorithm in two steps: Determine the expression of $p\left(\mathit{f}|\mathit{g},{\widehat{v}}_{ϵ},{\widehat{v}}_{\mathit{f}}\right)$ and using the values of in the previous iteration, we can propose an estimate for f, and then, using $p\left(v_{ϵ}|\mathit{g},\mathit{f}\right)$ and $p\left(v_{\mathit{f}}|\mathit{g},\mathit{f}\right),$ we can give estimates for ${\widehat{v}}_{ϵ}$ and ${\widehat{v}}_{\mathit{f}}$ which can again be used in the first step. It is interesting to know that all the three approaches of JMAP, GEM and VBA for this cas follow exactly this same iterative algorithm. The only differences will be in the update values of ${\stackrel{\sim }{\alpha }}_{ϵ},{\stackrel{\sim }{\beta }}_{ϵ},{\stackrel{\sim }{\alpha }}_{\mathit{f}},{\stackrel{\sim }{\beta }}_{\mathit{f}}$ and the choice of the estimators (MAP or PM) of ${\widehat{v}}_{ϵ}$ and ${\widehat{v}}_f.$

This case is also summarized in Figure 10.

Figure 10.

Bayesian inference scheme in linear systems and Gaussian priors. The posterior is also Gaussian, and all the computations can be done analytically.

7. Non-Gaussian priors

Very often, assuming that the noise is Gaussian is valid in many applications, but a Gaussian prior may not be adequate. Thus, the case of Non-Gaussian priors is of great importance. A very well known example is the case of Generalized Gaussian:

$$p\left(\mathit{f}\right)\propto exp\left[-\gamma \sum _j{\left|{\mathit{f}}_j\right|}^{\beta }\right].$$E46

The case of $\beta =2$ is the Gaussian case, $\beta >2$ gives the Super-Gaussian and $\beta <1$ is called Sub-Gaussian. Its particular case $\beta =1$ results to what is called Double Exponential (DE) prior law:

$$p\left(\mathit{f}\right)\propto exp\left[-\gamma \sum _j|{\mathit{f}}_j|\right]\propto exp\left[-\gamma {‖\mathit{f}‖}_1\right]$$E47

which, when using with a Gaussian likelihood, results to:

$$\{\begin{array}{ll}p\left(\mathit{f}|\mathit{g}\right)& \propto exp\left[-\frac{1}{2v_{ϵ}}J\left(\mathit{f}\right)\right]\text{with}\\ J\left(\mathit{f}\right)& =\frac{1}{2}{‖\mathit{g}-\mathit{H}\mathit{f}‖}^2+\lambda {‖\mathit{f}‖}_1,\lambda =\gamma v_{ϵ}.\end{array}$$E48

From this, we can see that the computation of MAP solution needs an appropriate optimization algorithm and the computations of the Posterior Mean (PM) or Posterior Covariance (PCov) or any other expectations become more difficult. However, as we will see later, VBA can be used to do approximate computations.

Another example is the Total Variation (TV) regularization method [11, 12, 13, 14] which can be interpreted as choosing the prior

$$p\left(\mathit{f}\right)\propto exp\left[-\gamma \sum _j|{\mathit{f}}_j-{\mathit{f}}_{j-1}|\right]\propto exp\left[-\gamma {‖\mathit{D}\mathit{f}‖}_1\right]$$E49

where D is the first order difference matrix.

This prior with a Gaussian model for noise results to:

$$\{\begin{array}{ll}p\left(\mathit{f}|\mathit{g}\right)& \propto exp\left[-\frac{1}{2\mathit{vϵ}}J\left(\mathit{f}\right)\right]\text{with}\\ J\left(\mathit{f}\right)=\frac{1}{2}{‖\mathit{g}-\mathit{H}\mathit{f}‖}^2+\lambda {‖\mathit{D}\mathit{f}‖}_1,\lambda =\gamma v_{ϵ}.\end{array}$$E50

One last example is using the Cauchy or more generally the Student-t distribution as the prior:

$$p\left(\mathit{f}\right)\propto exp\left[-\gamma \sum _{j}ln{\left(1+{\mathit{f}}_j^2\right)}^{\nu /2}\right]$$E51

which results to:

$$\{\begin{array}{ll}p\left(\mathit{f}|\mathit{g}\right)& \propto exp\left[-\frac{1}{2v_{ϵ}}J\left(\mathit{f}\right)\right]\text{with}\\ J\left(\mathit{f}\right)& =\frac{1}{2}{‖\mathit{g}-\mathit{H}\mathit{f}‖}^2+\lambda \sum _{j}ln{\left(1+{\mathit{f}}_j^2\right)}^{\nu /2},\lambda =\gamma v_{ϵ}.\end{array}$$E52

These three examples are of great importance. They have been used in the framework of MAP estimation and thus the optimization of the criteria J(f) for many linear inverse problems. However, the computation of other Bayesian estimators and uncertainty quantification (UQ) need again specific approximate solutions.

8. Hierarchical prior models

Even if simple Gaussian and non-Gaussian priors used in previous sections are of great importance and use in many applications, still they have, in many cases, limitations. For example, when we know that the signals have impulsive shapes or discontinuous or are piecewise continuous. The same limitations when we know, for example, that the images are composed of homogeneous regions with specified contours, or even, that the object under the test is composed of a limited number of homogeneous materials. Hierarchical models push farther these limitations of simple prior models. In the following, we consider three families of such hierarchical models: Sparsity aware models, Scaled Mixture models and Gauss-Markov-Potts models [10, 15, 16, 17].

8.1 Sparsity awarded hierarchical models

An easy way to consider the hierarchical sparsity awarded priors is to introduce a hidden variable, z and so consider the following Forward and prior models:

$$\{\begin{array}{l}\mathit{g}=\mathit{H}\mathit{f}+\mathit{ϵ},\\ \mathit{f}=\mathit{Dz}+\mathit{\zeta },\mathit{z}\text{sparse modeled}\mathrm{by}\text{Double}\mathrm{Exp}\left(\mathrm{DE}\right)\end{array}$$E53

with

$$\{\begin{array}{l}p\left(\mathit{g}|\mathit{f}\right)=\mathcal{N}\left(\mathit{g}|\mathit{H}\mathit{f},v_{ϵ}\mathit{I}\right)\\ p\left(\mathit{f}|\mathit{z}\right)=\mathcal{N}\left(\mathit{f}|\mathit{Dz},v_{\xi }\mathit{I}\right)\to \\ p\left(\mathit{z}\right)=\mathcal{DE}\left(\mathit{f}|\gamma \right)\propto exp\left[-\gamma {‖\mathit{z}‖}_1\right]\end{array}$$E54

Then, we have to find the expression of the joint posterior law $p\left(\mathit{f},\mathit{z}|\mathit{g}\right):$

$$\{\begin{array}{l}p\left(\mathit{f},\mathit{z}|\mathit{g}\right)\propto exp\left[-J\left(\mathit{f},\mathit{z}\right)\right]\text{with}\\ J\left(\mathit{f},\mathit{z}\right)=\frac{1}{2v_{ϵ}}{‖\mathit{g}-\mathit{H}\mathit{f}‖}_2^2+\frac{1}{2v_{\xi }}{‖\mathit{f}-\mathit{Dz}‖}_2^2+\gamma {‖\mathit{z}‖}_1\end{array}$$E55

from which we can infer on f and z [10, 16, 18, 19, 20, 21, 22].

For the unsupervised case, we can add the appropriate priors:

$$\{\begin{array}{l}p\left(\gamma \right)=\mathcal{IG}\left(\gamma |{\alpha }_{{\gamma }_z},{\beta }_{{\gamma }_z}\right)\\ p\left(v_{ϵ}\right)=\mathcal{IG}\left(v_{ϵ}|{\alpha }_{{ϵ}_0},{\beta }_{{ϵ}_0}\right)\\ p\left(v_{\xi }\right)=\mathcal{IG}\left(v_{\xi }|{\alpha }_{{\xi }_z},{\beta }_{{\xi }_z}\right)\end{array}$$E56

and thus obtain:

$$\{\begin{array}{l}\begin{array}{l}p\left(\mathit{f},\mathit{z},\gamma ,v_{ϵ},v_{\xi }|\mathit{g}\right)\propto exp\left[-J\left(\mathit{f},\mathit{z},\gamma ,v_{ϵ},v_{\xi }\right)\right]\text{with}\\ J\left(\mathit{f},\mathit{z},v_{ϵ},v_{\xi },\gamma \right)=\frac{1}{2v_{ϵ}}{‖\mathit{g}-\mathit{H}\mathit{f}‖}_2^2+\frac{1}{2v_{\xi }}{‖\mathit{f}-\mathit{Dz}‖}_2^2+\gamma {‖\mathit{z}‖}_1+\end{array}\\ \left({\alpha }_{{\gamma }_{\mathit{z}}}+n/2\right)ln\gamma +{\beta }_{{\gamma }_{\mathit{z}}}/\gamma +\\ \left({\alpha }_{{ϵ}_0}+m/2\right)ln{v}_{ϵ}+{\beta }_{{ϵ}_0}/v_{ϵ}+\\ \left({\alpha }_{{\xi }_{\mathit{z}}}+n/2\right)ln{v}_{\xi }+{\beta }_{{\xi }_{\mathit{z}}}/v_{\xi }\end{array}$$E57

It is interesting to note that the alternate optimization of this criterion gives the ADMM like algorithms [23, 24, 25] with the main advantage that here we have direct updates of the hyperparameters.

8.2 Scaled mixture models

Scaled Gaussian Mixture (SGM) models have been used in many applications to model rare events by their heavier tails with respect to Gaussian. They are also used in sparse signals modeling. A general SGM is defined as follows:

$$\mathcal{S}\left(\mathit{f}\right)=\int \mathcal{N}\left(\mathit{f}|0,v\right)p_m\left(v|\mathit{\theta }\right)\mathrm{d}v$$E58

where the variance of the Gaussian model $\mathcal{N}\left(\mathit{f}|0,v\right)$ is assumed to follow the mixing probability law $p_m\left(v|\mathit{\theta }\right).$ Between many possibilities for this mixing pdf is Inverse-Gamma which results to Student-t:

$$\mathcal{S}\left(\mathit{f}|\nu \right)=\int \mathcal{N}\left(\mathit{f}|0,v\right)\mathcal{IG}\left(v|\nu ,\nu \right)\mathrm{d}v$$E59

which have been extended to more general case:

$$\mathcal{S}\left(\mathit{f}|\alpha ,\beta \right)=\int \mathcal{N}\left(\mathit{f}|0,v\right)\mathcal{IG}\left(v|\alpha ,\beta \right)\mathrm{d}v$$E60

This pdf models have been used with success in many developments in Bayesian approach for inverse problems by:

$$p\left(\mathit{f}|\alpha ,\beta \right)=\underset{j}{\Pi }\int \mathcal{N}\left({\mathit{f}}_j|0,v_j\right)\mathcal{IG}\left(v_j|\alpha ,\beta \right)\mathrm{d}{v}_j$$E61

or

$$p\left(\mathit{f}|\alpha ,\beta \right)=\int \mathcal{N}\left(\mathit{f}|0,\mathit{v\Sigma }\right)\mathcal{IG}\left(v|\alpha ,\beta \right)\mathrm{d}v$$E62

Scaled Gaussian mixture models have been used extensively for modeling sparse signals. However, it happens very often that the signals or images are not sparse directly, but their gradients are, or more generally in a transformed domain such as Fourier or Wavelet domains. We have used these models extensively in hierarchical way:

$$\{\begin{array}{l}\mathit{g}=\mathit{H}\mathit{f}+\mathit{ϵ},\\ \mathit{f}=\mathit{Dz}+\zeta ,\mathit{z}\text{sparse Student}\to \{\begin{array}{l}p\left({\mathit{z}}_j|v_{{\mathit{z}}_j}\right)=\mathcal{N}\left({\mathit{z}}_j|0,v_{{\mathit{z}}_j}\right),\\ p\left(v_{{\mathit{z}}_j}\right)=\mathcal{IG}\left(v_{{\mathit{z}}_j}|{\alpha }_{{\mathit{z}}_0},{\beta }_{{\mathit{z}}_0}\right)\end{array}\end{array}$$E63

where D represents any linear transformations and ${\mathit{D}}^{–1}$ applied of f transforms it to a sparse vector z.

The whole relations of the likelihood and priors are summarized in below:

$$\{\begin{array}{l}p\left(\mathit{g}|\mathit{f}\right)=\mathcal{N}\left(\mathit{g}|\mathit{H}\mathit{f},v_{ϵ}\mathit{I}\right)\\ p\left(\mathit{f}|\mathit{z}\right)=\mathcal{N}\left(\mathit{f}|\mathit{Dz},v_{\xi }\mathit{I}\right)\\ p\left(\mathit{z}|{\upsilon }_{\mathit{z}}\right)=\mathcal{N}\left(\mathit{z}|0,V_{\mathit{z}}\right)\\ p\left({\upsilon }_{\mathit{z}}\right)=\underset{j}{\Pi }\mathcal{IG}\left(v_{{\mathit{z}}_j}|{\alpha }_{{\mathit{z}}_0},{\beta }_{{\mathit{z}}_0}\right)\\ p\left(v_{ϵ}\right)=\mathcal{IG}\left(v_{ϵ}|{\alpha }_{{ϵ}_0},{\beta }_{{ϵ}_0}\right)\\ p\left(v_{\xi }\right)=\mathcal{IG}\left(v_{\xi }|{\alpha }_{{\xi }_{\mathit{z}}},{\beta }_{{\xi }_{\mathit{z}}}\right)\end{array}$$E64

and the corresponding joint posterior of all the unknowns writes:

$$\{\begin{array}{l}\begin{array}{l}p\left(\mathit{f},\mathit{z},{\mathit{\upsilon }}_{\mathit{z}},v_{ϵ},v_{\xi }|\mathit{g}\right)\propto exp\left[-J\left(\mathit{f},\mathit{z},{\mathit{\upsilon }}_{\mathit{z}},v_{ϵ},v_{\xi }\right)\right]\\ J\left(\mathit{f},\mathit{z},{\mathit{\upsilon }}_{\mathit{z}},v_{ϵ},v_{\xi }\right)=\frac{1}{2v_{ϵ}}{‖\mathit{g}-\mathit{H}\mathit{f}‖}_2^2+\frac{1}{2v_{\xi }}{‖\mathit{f}-\mathit{Dz}‖}_2^2+{‖{{\mathit{V}}_{\mathit{z}}}^{\frac{-1}{2}}\mathit{z}‖}_2^2+\end{array}\\ {\sum }_j\left({\alpha }_{{\mathit{z}}_0}+1\right)ln{v}_{{\mathit{z}}_j}+{\beta }_{{\mathit{z}}_0}/v_{{\mathit{z}}_j}+\\ \left({\alpha }_{{ϵ}_0}+m/2\right)ln{v}_{ϵ}+{\beta }_{{ϵ}_0}/v_{ϵ}+\left({\alpha }_{{\xi }_{\mathit{z}}}+n/2\right)ln{v}_{\xi }+{\beta }_{{\xi }_{\mathit{z}}}/v_{\xi }\end{array}$$E65

Looking at this expression, we see that we have:

• Quadratic optimization with respect to f and z;

• Direct analytical expressions for the updates of the hyperparameters $v_{ϵ}$ and $v_{\xi };$

• Possibility to compute posterior mean and quantify uncertainties analytically via VBA.

A final case we consider is the case of Non-stationary noise and sparsity enforcing prior in the same framework.

$$\{\begin{array}{ll}\mathit{g}=\mathit{H}\mathit{f}+\mathit{ϵ},& \mathit{ϵ}\text{non stationary}\to \left\{\begin{array}{l}p\left({ϵ}_i|v_{{ϵ}_i}\right)=\mathcal{N}\left({ϵ}_i|0,v_{{ϵ}_i}\right),\\ p\left(v_{{ϵ}_i}\right)=\mathcal{IG}\left(v_{{ϵ}_i}|{\alpha }_{{ϵ}_0},{\beta }_{{ϵ}_0}\right)\end{array}\right.\\ \mathit{f}=\mathit{Dz}+\mathit{\zeta },& \mathit{z}\text{sparse Student}\to \left\{\begin{array}{l}p\left(z_j|v_{z_j}\right)=\mathcal{N}\left(z_j|0,v_{z_j}\right),\\ p\left(v_{z_j}\right)=\mathcal{IG}\left(v_{z_j}|{\alpha }_{z_0},{\beta }_{z_0}\right)\end{array}\right.\end{array}$$E66

Again here, all the expressions of likelihood and priors can be summarized as follows:

$$\{\begin{array}{l}p\left(\mathit{g}|\mathit{f}\right)=\mathcal{N}\left(\mathit{g}|\mathit{H}\mathit{f},V_{ϵ}\right)\\ p\left(\mathit{f}|\mathit{z}\right)=\mathcal{N}\left(\mathit{f}|\mathit{D}\mathit{z},v_{\xi }I\right)\\ p\left(\mathit{z}|{\mathit{\upsilon }}_{\mathit{z}}\right)=\mathcal{N}\left(\mathit{z}|0,V_{\mathit{z}}\right)\\ p\left({\mathit{\upsilon }}_{\mathit{z}}\right)=\underset{j}{\Pi }\mathcal{IG}\left(v_{{\mathit{z}}_j}|{\alpha }_{{\mathit{z}}_0},{\beta }_{{\mathit{z}}_0}\right)\\ p\left({\mathit{\upsilon }}_{ϵ}\right)=\underset{i}{\Pi }\mathcal{IG}\left(v_{{ϵ}_i}|{\alpha }_{{ϵ}_0},{\beta }_{{ϵ}_0}\right)\\ p\left(v_{\xi }\right)=\mathcal{IG}\left(v_{\xi }|{\alpha }_{{\xi }_{\mathit{z}}},{\beta }_{{\xi }_{\mathit{z}}}\right)\end{array}$$E67

and the joint posterior of all the unknowns become:

$$\{\begin{array}{l}\begin{array}{l}p\left(\mathit{f},\mathit{z},{\mathit{\upsilon }}_{\mathit{z}},{\mathit{\upsilon }}_{ϵ},v_{\xi }|\mathit{g}\right)\propto exp\left[-J\left(\mathit{f},\mathit{z},{\mathit{\upsilon }}_{\mathit{z}},{\mathit{\upsilon }}_{ϵ},v_{\xi }\right)\right]\\ J\left(\mathit{f},\mathit{z},{\mathit{\upsilon }}_{\mathit{z}},v_{ϵ},v_{\xi }\right)={‖{\mathit{V}}_{ϵ}^{\frac{-1}{2}}\left(\mathit{g}-\mathit{H}\mathit{f}\right)‖}_2^2+\frac{1}{2v_{\xi }}{‖\mathit{f}-\mathit{D}\mathit{z}‖}_2^2+{‖{{\mathit{V}}_{\mathit{z}}}^{\frac{-1}{2}}\mathit{z}‖}_2^2+\end{array}\\ {\sum }_j\left({\alpha }_{{\mathit{z}}_0}+1\right)ln{v}_{{\mathit{z}}_j}+{\beta }_{{\mathit{z}}_0}/v_{{\mathit{z}}_j}+\\ {\sum }_j\left({\alpha }_{{ϵ}_0}+1\right)ln{v}_{{ϵ}_i}+{\beta }_{{ϵ}_0}/v_{{ϵ}_i}+\\ \left({\alpha }_{{\xi }_{\mathit{z}}}+n/2\right)ln{v}_{\xi }+{\beta }_{{\xi }_{\mathit{z}}}/v_{\xi }\end{array}$$E68

The following scheme shows graphically this case.

8.3 A four level hierarchical model

To account separately for the measurement and forward modeling error, a more detailed and four level hierarchical model has been proposed:

$$\{\begin{array}{l}\mathit{g}={\mathit{g}}_0+\mathit{ϵ},\text{measurement error}\\ {\mathit{g}}_0=\mathit{H}\mathit{f}+\xi ,\text{modeling error}\\ \mathit{f}=\mathit{Dz}+\mathit{\zeta },\text{Prior model}\end{array}$$E69

which accounts for two terms of errors (variable splitting) and Sparsity enforcing in Transformed domain prior: $\mathit{f}=\mathit{Dz}+\mathit{\zeta }$ with z sparse, modeled itself by Normal-IG.

This model is presented graphically here.

In this model, there are three error terms: $\mathit{ϵ}$ the observation error, $\mathit{\zeta }$ the forward modeling error and $\mathit{\zeta }$ the transform domain modeling error. These are detailed in the following:

• $\mathit{g}={\mathit{g}}_0+\mathit{ϵ},$: $\mathit{ϵ}$ is assumed to be Gaussian:

$$p\left(\mathit{g}|{\mathit{g}}_0,v_{\mathit{ϵ}}\right)=\mathcal{N}\left(\mathit{g}|{\mathit{g}}_0,v_{\mathit{ϵ}}I\right),p\left(v_{\mathit{ϵ}}\right)=\mathcal{IG}\left(v_{\mathit{ϵ}}|{\alpha }_{{\mathit{ϵ}}_0},{\beta }_{{\mathit{ϵ}}_0}\right),$$

• ${\mathit{g}}_0=\mathit{H}\mathit{f}+\xi ,$$\xi$ is assumed to be Student-t:

$$\{\begin{array}{l}p\left({\mathit{g}}_0|\mathit{f},v_{\xi }\right)=\mathcal{N}\left({\mathit{g}}_0|\mathit{H}\mathit{f},{\mathit{V}}_{\xi }\right),{\mathit{V}}_{\xi }=diag\left[v_{\xi }\right],\\ p\left(v_{\xi }\right)={\Pi }_{i=1}^{M}p\left(v_{{\xi }_i}\right)={\Pi }_{i=1}^M\mathcal{IG}\left(v_{{\xi }_i}|{\alpha }_{{\xi }_{\mathit{z}}},{\beta }_{{\xi }_{\mathit{z}}}\right),\end{array}$$

• $\mathit{f}=\mathit{D}\mathit{z}+\mathit{\zeta },$$\mathit{\zeta }$ is assumed to be Gaussian:

$$p\left(\mathit{f}|\mathit{z},v_{\mathit{\zeta }}\right)=\mathcal{N}\left(\mathit{f}|\mathit{D}\mathit{z},v_{\mathit{\zeta }}I\right),p\left(v_{\mathit{\zeta }}\right)=\mathcal{IG}\left(v_{\mathit{\zeta }}|{\alpha }_{{\mathit{\zeta }}_0},{\beta }_{{\mathit{\zeta }}_0}\right),$$

• z is assumed to be sparse and thus modeled via Normal-IG:

$$\{\begin{array}{l}p\left(\mathit{z}|{\mathit{v}}_{\mathit{z}}\right)=\mathcal{N}\left(\mathit{z}|0,{\mathit{V}}_{\mathit{z}}\right),{\mathit{V}}_{\mathit{z}}=diag\left[{\mathit{v}}_{\mathit{z}}\right]\\ p\left({\mathit{v}}_{\mathit{z}}\right)={\Pi }_{j=1}^{N}p\left({\mathit{v}}_{{\mathit{z}}_j}\right)={\Pi }_{j=1}^N\mathcal{IG}\left({\mathit{v}}_{{\mathit{z}}_j}|{\alpha }_{{\mathit{z}}_0},{\beta }_{{\mathit{z}}_0}\right)\end{array}$$

which results in:

$$p\left(\mathit{f},{\mathit{g}}_0,\mathit{z},{\mathit{v}}_{\mathit{ϵ}},{\mathit{v}}_{\xi },{\mathit{v}}_{\mathit{z}}|\mathit{g}\right)\propto exp\left[-J\left(\mathit{f},{\mathit{g}}_0,\mathit{z},{\mathit{v}}_{\mathit{ϵ}},{\mathit{v}}_{\xi },{\mathit{v}}_{\mathit{z}}\right)\right]$$E70

with

$$\begin{array}{c}J\left(\mathit{f},{\mathit{g}}_0,\mathit{z},{\mathit{v}}_{ϵ},{\mathit{v}}_{\xi },{\mathit{v}}_{\mathit{z}}\right)=\frac{1}{2{\mathit{v}}_{ϵ}}{‖\mathit{g}-{\mathit{g}}_0‖}_2^2+\left({\alpha }_{{ϵ}_0}+1\right)ln{\mathit{v}}_{ϵ}+\frac{{\beta }_{{ϵ}_0}}{{\mathit{v}}_{ϵ}}\\ +\frac{1}{2}{‖{\mathit{V}}_{\xi }^{-1/2}\left({\mathit{g}}_0-\mathit{H}\mathit{f}\right)‖}_2^2+{\sum }_{i=1}^M\left[\left({\alpha }_{{\xi }_{\mathit{z}}}+1\right)ln{\mathit{v}}_{{\xi }_i}+\frac{{\beta }_{{\xi }_{\mathit{z}}}}{{\mathit{v}}_{{\xi }_i}}\right]\\ +\frac{1}{2{\mathit{v}}_{\mathit{\zeta }}}{‖\mathit{f}-\mathit{D}\mathit{z}‖}_2^2+\left({\alpha }_{{\mathit{\zeta }}_0}+1\right)ln{\mathit{v}}_{\xi }+\frac{{\beta }_{{\mathit{\zeta }}_0}}{{\mathit{v}}_{\xi }}\\ +\frac{1}{2}{‖{{\mathit{V}}_{\mathit{z}}}^{-1/2}\mathit{z}‖}_2^2+{\sum }_{j=1}^N\left[\left({\alpha }_{{\mathit{z}}_0}+1\right)ln{\mathit{v}}_{{\mathit{z}}_j}+\frac{{\beta }_{{\mathit{z}}_0}}{{\mathit{v}}_{{\mathit{z}}_j}}\right]\end{array}$$E71

Using then the JMAP approach with an alternate optimization strategy needs the following optimization steps:

• with respect to f: $J\left(\mathit{f}\right)=\frac{1}{2}{‖{\mathit{V}}_{\xi }^{-1/2}\left({\mathit{g}}_0-\mathit{H}\mathit{f}\right)‖}_2^2+\frac{1}{2{\mathit{v}}_{\mathit{\zeta }}}{‖\mathit{f}-\mathit{D}\mathit{z}‖}_2^2$

• with respect to g₀: $J\left({\mathit{g}}_0\right)=\frac{1}{2{\mathit{v}}_{ϵ}}{‖\mathit{g}-{\mathit{g}}_0‖}_2^2+\frac{1}{2}{‖{\mathit{V}}_{\xi }^{-1/2}\left({\mathit{g}}_0-\mathit{H}\mathit{f}\right)‖}_2^2$

• with respect to z: $J\left(\mathit{z}\right)=\frac{1}{2{\mathit{v}}_{\mathit{\zeta }}}{‖\mathit{f}-\mathit{D}\mathit{z}‖}_2^2+\frac{1}{2}{‖{\mathit{V}}_{\mathit{z}}^{-1/2}\mathit{z}‖}_2^2$

• with respect to ${\mathit{v}}_{ϵ}:$$J\left({\mathit{v}}_{ϵ}\right)=\frac{1}{2{\mathit{v}}_{ϵ}}{‖\mathit{g}-{\mathit{g}}_0‖}_2^2+\left({\alpha }_{{ϵ}_0}+1\right)ln{\mathit{v}}_{ϵ}+\frac{{\beta }_{{ϵ}_0}}{{\mathit{v}}_{ϵ}}$

• with respect to ${\mathit{v}}_{{\xi }_i}:$$J\left({\mathit{v}}_{{\xi }_i}\right)=\frac{1}{2}{‖{\mathit{V}}_{\xi }^{-1/2}\left({\mathit{g}}_0-\mathit{H}\mathit{f}\right)‖}_2^2+{\sum }_{i=1}^M\left[\left({\alpha }_{{\xi }_{\mathit{z}}}+1\right)ln{\mathit{v}}_{{\xi }_i}+\frac{{\beta }_{{\xi }_{\mathit{z}}}}{{\mathit{v}}_{{\xi }_i}}\right]$

• with respect to $v_{\mathit{\zeta }}:$$J\left(v_{\mathit{\zeta }}\right)=\frac{1}{2v_{\mathit{\zeta }}}{‖\mathit{f}-\mathit{D}\mathit{z}‖}_2^2+\left({\alpha }_{{\mathit{\zeta }}_0}+1\right)ln{v}_{\xi }+\frac{{\beta }_{{\mathit{\zeta }}_0}}{v_{\xi }}$

• to $v_{{\mathit{z}}_j}:$$J\left(v_{{\mathit{z}}_j}\right)=\frac{1}{2}{‖{{\mathit{V}}_{\mathit{z}}}^{-1/2}\mathit{z}‖}_2^2+{\sum }_{j=1}^N\left[\left({\alpha }_{{\mathit{z}}_0}+1\right)ln{v}_{{\mathit{z}}_j}+\frac{{\beta }_{{\mathit{z}}_0}}{v_{{\mathit{z}}_j}}\right]$

This approach has the following main advantages and limitations.

Advantages:

• All the optimization are either quadratic or explicit

• Quadratic optimizations can be done efficiently

• For great dimensional problems, the needed operators H, H′, D and D′ can be implemented on GPU

• For Computed Tomography, efficient GPU implementation of these operators have been done in our group for 2D and 3D CT.

Limitations:

• Huge amount of memory is needed for f, g₀, $v_{\xi }$ and v_z

• No easy way to study the global convergence of the algorithm.

• The number of hyper-hyperparameters $\left({\alpha }_{{ϵ}_0},{\beta }_{{ϵ}_0}\right),$$\left({\alpha }_{{\xi }_z},{\beta }_{{\xi }_z}\right),$$\left({\alpha }_{{\zeta }_0},{\beta }_{{\zeta }_0}\right),$$\left({\alpha }_{z_0},{\beta }_{z_0}\right)$ to be fixed is important. However, the results are not so sensitive to these parameters.

8.4 Gauss-Markov-Potts models

To introduce the Gauss-Markov-Potts model, let us have a look at the images in Figure 11.

Figure 11.

Different images with different characteristics in different imaging systems.

The question we want to answer is: which prior model can be more appropriate for these images? One way, to answer to this question is either look at the histogram of the pixels or the pixels of the gradient images. Then, if we take one typical line of the gradient or one typical line of the image itself and draw them as a 1D-signal, we obtain the cases in Figure 12.

Figure 12.

Different possible prior modeling in relation to the different images of Figure 11.

From these two figures, we see that for some cases, a Gaussian or generalized Gaussian may be very good models. But, for other cases, if we want to explicitly account for the presence of the contours, we can introduce a binary hidden variable to represent it. Finally, for the last example of the image in Figure 10 and its corresponding typical line in Figure 11, we need to introduce a hidden variable z which encodes the following fact that:

In NDT applications of CT, the objects are, in general, composed of a finite number of materials, and the voxels corresponding to each material are grouped in compact regions.

How to model this prior information?

To answer to this question, first consider such an image f (r) with its segmentation z(r) and contours q(r) as shown in Figure 13.

Figure 13.

An image of an object composed of homogeneous compact regions, its segmentation and the contours of thoses regions.

As it can be seen, we introduced two hidden variables z(r) and q(r), the first representing the segmentation and the second the contours of the image. z(r) takes the integer values $\left\{k=1,\cdots ,K\right\},$ each presented by a different color and q(r) a binary value {0, 1}. The second can easily be obtained from the first. So, from now, we consider only z(r).

As each value of z represents a homogeneous material, we can translate this by:

$$p\left(f\left(\mathit{r}\right)|z\left(\mathit{r}\right)=k,m_k,v_k\right)=\mathcal{N}\left(m_k,v_k\right)$$E72

encoding the fact that inside each homogeneous material, i.e.; all the pixels having $z\left(r\right)=k,$ represent a homogenous material characterized by the two parameters f (m_k, v_k). This results to:

$$p\left(f\left(\mathit{r}\right)\right)=\sum _{k}P\left(z\left(\mathit{r}\right)=k\right)\mathcal{N}\left(m_k,v_k\right)\text{Mixture of Gaussians}$$E73

which shows the mixture of Gaussian model of the pixel values. See also Figure 14.

Figure 14.

A metalic object with a default area inside it: Black pixels represent air, white pixels metal and gray pixels the defaults area. On left image these are codes by colors ($z=1$ represents air, $z=2$ represents metal and $z=3$ represents default area.

The next step is to propose a probability distribution for z. As we want a compactness of the regions, a Markov modeling is appropriate:

$$p\left(z\left(\mathit{r}\right)|z\left({\mathit{r}}^{\prime }\right),{\mathit{r}}^{\prime }\in \mathcal{V}\left(\mathit{r}\right)\right)\propto exp\left[-\gamma \sum _{{\mathit{r}}^{\prime }\in \mathcal{V}\left(\mathit{r}\right)}\delta \left(z\left(\mathit{r}\right)-z\left({\mathit{r}}^{\prime }\right)\right)\right]$$E74

A Potts Markov model is still more appropriate:

$$p\left(z\left(\mathit{r}\right),\mathit{r}\in \mathrm{\Omega }\right)\propto exp\left[-\gamma \sum _{\mathit{r}\in \mathrm{\Omega }{\mathit{r}}^{\prime }}\sum _{\in \mathcal{V}\left(\mathit{r}\right)}\delta \left(z\left(\mathit{r}\right)-z\left({\mathit{r}}^{\prime }\right)\right)\right]$$E75

where $\mathrm{\Omega }$ represents all pixels of the image.

Thus, to each pixel of the image is associated 2 variables f (r) and z(r) with the following possible properties:

• $\mathit{f}\mid \mathit{z}$ Gaussian iid, z iid: Mixture of Gaussians

• $\mathit{f}\mid \mathit{z}$ Gauss-Markov, z iid: Mixture of Gauss-Markov

• $\mathit{f}\mid \mathit{z}$ Gaussian iid, z Potts-Markov: Mixture of Independent Gaussians, (MIG with Hidden Potts)

• $\mathit{f}\mid \mathit{z}$ Markov, z Potts-Markov: Mixture of Gauss-Markov, (MGM with hidden Potts)

From these four different cases, we consider two which are illustrated in Figure 15.

Figure 15.

Two proposed gauss-Markov-Potts models used in many NDT applications.

Using the notations on this figure, and noting by f all the pixels of the image, by z all the pixels of the segmented image, and by $\mathit{\theta }$ all the parameters $\left\{v_{ϵ},\left({\alpha }_k,m_k,v_k\right),k=1,\cdots ,K\right\},$ we can write:

$$p\left(\mathit{f},\mathit{z},\mathit{\theta }|\mathit{g}\right)\propto p\left(\mathit{g}|\mathit{f},{\mathit{v}}_{ϵ}\right)p\left(\mathit{f}|\mathit{z},\mathit{m},\mathit{v}\right)p\left(\mathit{z}|\gamma ,\alpha \right)p\left(\mathit{\theta }\right)$$E76

where

$$\mathit{m}=\left\{m_k,,,k=1,,,\cdot ,,,K\right\},v=\left\{v_k,,,k=1,,,\cdot ,,,K\right\},\alpha =\left\{{\alpha }_k,,,k=1,,,\cdot ,,,K\right\},\mathit{\theta }=\left\{v_{ϵ},,,m,,,v,,,\mathit{\text{alphab}}\right\}$$

The expressions of $p\left(\mathit{g}|\mathit{f},v_{ϵ}\right),$$p\left(\mathit{f}|\mathit{z},\mathit{m},\mathit{\upsilon }\right)$ and $p\left(\mathit{z}|\gamma ,\mathit{\alpha }\right)$ have been given before. We need to define $p\left(\mathit{\theta }\right)$ which can be chosen as the conjugate priors: Dirichlet for $\alpha ,$ Gaussian for m and Inverse-Gamma for all the variances.

Direct computation and use of $p\left(\mathit{f},\mathit{z},\mathit{\theta }|\mathit{g};\mathcal{M}\right)$ is too complex, because we do not have analytical expression for the proportionality term of the joint probability law:

$$p\left(\mathit{f},\mathit{z},\mathit{\theta }|\mathit{g}\right)\propto p\left(\mathit{g}|\mathit{f},\mathit{z},\mathit{\theta }\right)p\left(\mathit{f}|\mathit{z},\mathit{\theta }\right)p\left(\mathit{z}\right)p\left(\mathit{\theta }\right)$$E77

As we have three sets of variables f, z and $\mathit{\theta },$ we can use different schemes, for example a Gibbs sampling scheme:

$$\widehat{\mathit{f}}\sim p\left(\mathit{f}|\widehat{\mathit{z}},\widehat{\mathit{\theta }},\mathit{g}\right)\to \widehat{\mathit{z}}\sim p\left(\mathit{z}|\widehat{\mathit{f}},\widehat{\mathit{\theta }},\mathit{g}\right)\to \widehat{\mathit{\theta }}\sim \left(\mathit{\theta }|\widehat{\mathit{f}},\widehat{\mathit{z}},\mathit{g}\right)$$E78

with:

• Sample f from $p\left(\mathit{f}|\widehat{\mathit{z}},\widehat{\mathit{\theta }},\mathit{g}\right)\propto p\left(\mathit{g}|\mathit{f},\mathit{\theta }\right)p\left(\mathit{f}|\widehat{\mathit{z}},\widehat{\mathit{\theta }}\right)$

Needs optimisation of a quadratic criterion.

• Sample z from $p\left(\mathit{z}|\widehat{\mathit{f}},\widehat{\mathit{\theta }},\mathit{g}\right)\propto p\left(\mathit{g}|\widehat{\mathit{f}},\widehat{\mathit{z}},\widehat{\mathit{\theta }}\right)p\left(\mathit{z}\right)$

Needs sampling of a Potts Markov field.

• Sample $\mathit{\theta }$ from $p\left(\mathit{\theta }|\widehat{\mathit{f}},\widehat{\mathit{z}},\mathit{g}\right)\propto p\left(\mathit{g}|\widehat{\mathit{f}},{\sigma }_{ϵ}^{2}I\right)p\left(\widehat{\mathit{f}}|\widehat{\mathit{z}},\left(m_k,v_k\right)\right)p\left(\mathit{\theta }\right)$

More details and other schemes such as JMAP and VBA can be found in Refs. [26].

To illustrate an example of application, we considered a NDT application, where a metalic object is tested to detect a default inside it. As, the problem was, not only to detect the default, but also to characterize its shape and size, an X-ray computed tomography (CT) with only two projections is proposed and used. This problem is illustrated in Figure 16.

Figure 16.

A non destructive testing (NDT) application where f (x, y) has to be reconstructed from its marginals g₁(x) and g₂(y).

The mathematical part of this very ill-posed inverse problem is the following:

Given the functions g₁(x) and g₂(y) find the image f (x, y).

This problem also arise in probability theory and statistics, where f (x, y) is a joint distribution and g₁(x) and g₂(y) its two marginals. We know that this problem has infinite number of solutions: $f\left(x,y\right)=g_1\left(x\right)g_2\left(y\right)\mathrm{\Omega }\left(x,y\right)$ where $\mathrm{\Omega }\left(x,y\right)$ is called a Copula:

$$\int \mathrm{\Omega }\left(x,y\right)\mathrm{d}x=1\text{and}\int \mathrm{\Omega }\left(x,y\right)\mathrm{d}y=1$$

So, any arbitrary copula function defines a solution. The problem is ill-posed and we need to use any possible prior information to try to obtain a unique or acceptable solution. The probabilistic solution we proposed is illustrated in Figure 17.

Figure 17.

Probabilistic Bayesian method for the NDT image resonstruction problem.

Unsupervised Bayesian estimation:

$$p\left(\mathit{f},\mathit{z},\mathit{\theta }|\mathit{g}\right)\propto p\left(\mathit{g}|\mathit{f},\mathit{z},\mathit{\theta }\right)p\left(\mathit{f}|\mathit{z},\mathit{\theta }\right)p\left(\mathit{\theta }\right)$$

A summary of the results is given in Figure 18 where the proposed method result.

Figure 18.

Probabilistic Bayesian method for the NDT image resonstruction problem. a) Shows the original image f, b) is the result of Back-projection, c) is the result of filtered Back-projection, d) and e) are the result of a Markov model with hidden line process, and f), g) and h) show the results of the gauss-Markov-Potts method.

9. Bayesian computation

As we could see, very often, we can find the expression of the posterior law $p\left(\mathit{f}|g\right),$ sometimes exactly as is the case of the linear models with Gaussian priors in the previous section, but often up to the normalization constant (the evidence term) p(g) in:

$$p\left(\mathit{f}|\mathit{g}\right)=\frac{1}{p\left(\mathit{g}\right)}p\left(\mathit{g}|\mathit{f}\right)p\left(\mathit{f}\right)=\frac{1}{p\left(\mathit{g}\right)}p\left(\mathit{g},\mathit{f}\right).$$E79

This term is not necessary for Maximum A Posteriori (MAP) but it is needed for Expected A Posteriori (EAP) and for doing any other expectation computation.

This is the case, almost in all Non-Gaussian prior models or Non-Gaussian noise models or the Non-Linear forward models. In this chapter, a few cases are considered more in detail. Even in the Gaussian and linear case which is the simplest case, and we have analytical expressions for almost everything, the computational cost for large scale problems brings us to search for approximate but fast solutions.

9.1 Large scale linear and Gaussian models

As we could see in previous chapter, the linear forward model $\mathit{g}=\mathit{H}\mathit{f}+\mathit{ϵ}$ with Gaussian noise and Gaussian prior is the simplest case where we can do all the computations analytically.

$$\{\begin{array}{l}p\left(\mathit{g}|\mathit{f}\right)=\mathcal{N}\left(\mathit{g}|\mathit{H}\mathit{f},v_{ϵ}\mathit{I}\right)\\ p\left(\mathit{f}\right)=\mathcal{N}\left(\mathit{f}|{\mathit{f}}_0,v_{\mathit{f}}\mathit{I}\right)\end{array}\to \{\begin{array}{l}p\left(\mathit{g}\right)=\mathcal{N}\left(\mathit{g}|\mathit{H}{\mathit{f}}_0,v_{\mathit{f}}\mathit{H}{\mathit{H}}^{\prime }+v_{ϵ}\mathit{I}\right),\\ p\left(\mathit{f}|\mathit{g}\right)=\mathcal{N}\left(\mathit{f}|\widehat{\mathit{f}},\widehat{\mathrm{\Sigma }}\right)\text{with}:\\ \widehat{\mathit{f}}={\mathit{f}}_0+{\left[{\mathit{H}}^{\prime }\mathit{H}+\mathit{\lambda I}\right]}^{-1}{\mathit{H}}^{\prime }\left(\mathit{g}-\mathit{H}{\mathit{f}}_0\right)\\ \widehat{\mathrm{\Sigma }}=v_{ϵ}{\left[{\mathit{H}}^{\prime }\mathit{H}+\lambda \mathit{I}\right]}^{-1}\lambda =\frac{v_{ϵ}}{v_{\mathit{f}}}\end{array}$$E80

The trick here is, for example, for computing $\widehat{\mathit{f}}$ to use the fact that

$$\{\begin{array}{l}p\left(\mathit{g}|\mathit{f}\right)\propto exp\left[-\frac{1}{2v_{ϵ}}{‖\mathit{g}-\mathit{H}\mathit{f}‖}_2^2\right]\\ p\left(\mathit{f}\right)\propto exp\left[-\frac{1}{2v_{\mathit{f}}}{‖\mathit{f}‖}_2^2\right]\\ p\left(\mathit{f}|\mathit{g}\right)\propto exp\left[-\frac{1}{2v_{ϵ}}J\left(\mathit{f}\right)\right]\text{with}J\left(\mathit{f}\right)=\frac{1}{2}{‖\mathit{g}-\mathit{H}\mathit{f}‖}^2+\lambda {‖\mathit{f}‖}_2^2,\lambda =\frac{v_{ϵ}}{v_{\mathit{f}}}\end{array}$$E81

and, as the mean and the mode of a Gaussian probability law are the same, we can use:

$$\widehat{\mathit{f}}=\underset{\mathit{f}}{\text{argmax}}\left\{J\left(\mathit{f}\right)\right\}\text{with}J\left(\mathit{f}\right)={‖\mathit{g}-\mathit{H}\mathit{f}‖}^2+\lambda {‖\mathit{f}‖}^2$$E82

and so the problem can be cast as an optimization problem of a quadratic criterion for which there are a great number of algorithms. Let here to show the simplest one which is the gradient based and so needs the expression of the gradient:

$$\nabla J\left(\mathit{f}\right)=-2{\mathit{H}}^{\prime }\left(\mathit{g}-\mathit{H}\mathit{f}\right)+2\mathit{\lambda }\mathit{f}$$E83

which can be summarized as follows:

$$\{\begin{array}{l}{\mathit{f}}^{\left(0\right)}=0\\ {\mathit{f}}^{\left(k+1\right)}={\mathit{f}}^{\left(k\right)}+\alpha \left[{\mathit{H}}^{\prime }\left(\mathit{g}-\mathit{H}{\mathit{f}}^{\left(k\right)}\right)+2\lambda {\mathit{f}}^{\left(k\right)}\right]\end{array}$$E84

As we can see, at each iteration, we need to be able to compute the forward operationH f and the backward operation${\mathit{H}}^{\prime }\mathit{\delta }\mathit{g}$ where $\mathit{\delta g}=\mathit{g}–\mathit{H}\mathit{f}.$ This optimization algorithm needs to write two programs:

• Forward operation H f

• Adjoint operation ${\mathit{H}}^{\prime }\mathit{\delta }\mathit{g}$

These two operations can be implemented using High Performance parallel processors such as Graphical Processor Units (GPU).

The computation of the posterior covariance is much more difficult. There are a few methods: The first category is the methods which use the particular structure of the matrix H and H′H or HH′ as we can use the matrix inversion lemma and see that

$$\widehat{\mathrm{\Sigma }}=v_{ϵ}{\left[{\mathit{H}}^{\prime }\mathit{H}+\mathit{\lambda }\mathit{I}\right]}^{-1}=v_{ϵ}\mathit{I}-{\mathit{H}}^{\prime }{\left[\mathit{H}{\mathit{H}}^{\prime }+{\lambda }^{-1}\mathit{I}\right]}^{-1}$$E85

For example, in a signal deconvolution problem, the matrix H has a Toeplitz structure and so have the matrices H′H and HH′ which can be approximated by Circulant matrices and be diagonalized using the Fourier Transform.

The second, more general, is to approximate $\widehat{\mathrm{\Sigma }}$ by a diagonal matrix, which can also be interpreted as to approximate the posterior law $p\left(\mathit{f}|\mathit{g}\right)$ by a separable $q\left(\mathit{f}\right)={\Pi }_{j}q\left({\mathit{f}}_j\right).$ This brings us naturally to the Approximate Bayesian Computation (ABC). But, before going to the details of ABC methods, let consider the case where the hyperparameters of the problem (parameters of the prior laws) are also unknown.

To be able to do the computation, we need mainly to compute the determinant of the matrix $v_f\mathit{H}{\mathit{H}}^{\prime }+v_{ϵ}\mathit{I}$ for p(g) and the inverse of the matrices $\left[{\mathit{H}}^{\prime }\mathit{H}+\mathit{\lambda }\mathit{I}\right]$ or $\left[\mathit{H}{\mathit{H}}^{\prime }+{\lambda }^{-1}\mathit{I}\right]$.

9.2 Large scale computation of the posterior covariance

Computing the determinant of the matrix $v_{\mathit{f}}\mathit{H}{\mathit{H}}^{\prime }+v_{ϵ}\mathit{I}$ for p(g) and the inverse of the matrices $\left[{\mathit{H}}^{\prime }\mathit{H}+\mathit{\lambda }\mathit{I}\right]$ or, $\left[\mathit{H}{\mathit{H}}^{\prime }+{\lambda }^{-1}\mathit{I}\right]$ which are needed for uncertainty quantification, are between the greatest subjects of open research for Big Data problems. Here, we consider a few cases.

9.2.1 Structured matrices

One solution is to use the particular structure of these matrices when possible. This is the case for deconvolution or image restoration, where these matrices have Toeplitz or Bloc-Toeplitz structures which can be well approximated by Circulant or Bloc-Circulant matrices and diagonalized using Fourier Transform (FT) and Fast FT (FFT). The main idea here is using the properties of the circulant matrices: If H is a circulant matrix, then

$$\mathit{H}=\mathit{F}\mathrm{\Lambda }{\mathit{F}}^{\prime }$$E86

where F is the DFT or FFT matrix and F′ the IDFT or IFFT and $\mathrm{\Lambda }$ is a diagonal matrix whose elements are the FT of the first line of the circulent matrix. As the first line of that circulent matrix contains the samples of the impulse response, the vector of the diagonal elements represents the spectrum of the impulse response (transfer function). Using this property, we have:

$${\left[{\mathit{H}}^{\prime }\mathit{H}+\mathit{\lambda }\mathit{I}\right]}^{-1}={\left[{\mathit{F}}^{\prime }\mathrm{\Lambda }\mathit{F}{\mathit{F}}^{\prime }\mathrm{\Lambda }+\lambda \right]}^{-1}={\left[{\mathit{F}}^{\prime }{\mathrm{\Lambda }}^2\mathit{F}+\mathit{\lambda }\mathit{I}\right]}^{-1}=\mathit{F}{\left[{\mathrm{\Lambda }}^2+\mathit{\lambda }\mathit{I}\right]}^{-1}{\mathit{F}}^{\prime }$$E87

9.2.2 Sampling based methods

Second solution is generating samples from the posterior law and use them to compute the variances and covariances. So, the problem is how to generate a sample from the posterior law

$$\{\begin{array}{l}p\left(\mathit{f}|\mathit{g}\right)=\mathcal{N}\left(\mathit{f}|\widehat{\mathit{f}},\widehat{\mathrm{\Sigma }}\right)\text{with}:\\ \widehat{\mathit{f}}={\mathit{f}}_0+{\left[{\mathit{H}}^{\prime }\mathit{H}+\mathit{\lambda }\mathit{I}\right]}^{-1}{\mathit{H}}^{\prime }\left(\mathit{g}-\mathit{H}{\mathit{f}}_0\right)\\ \widehat{\mathrm{\Sigma }}=v_{ϵ}{\left[{\mathit{H}}^{\prime }\mathit{H}+\mathit{\lambda }\mathit{I}\right]}^{-1},\lambda =\frac{v_{ϵ}}{v_{\mathit{f}}}\end{array}$$E88

One solution is to compute the Cholesky decomposition of the covariance matrix $\widehat{\mathrm{\Sigma }}=\mathit{A}{\mathit{A}}^{\prime },$ generate a vector, $\mathit{u}\sim \mathcal{N}\left(\mathit{u}|0,\mathit{I}\right)$ and then generate a sample $\mathit{f}=\mathit{Au}+\widehat{\mathit{f}}$ [27]. We can compute $\widehat{\mathit{f}}$ by optimizing

$$J\left(\mathit{f}\right)=\frac{1}{2}{‖\mathit{g}-\mathit{H}\mathit{f}‖}^2+\lambda {‖\mathit{f}-{\mathit{f}}_0‖}_2^2,\lambda =\frac{v_{ϵ}}{v_{\mathit{f}}},$$E89

but the main computational cost is the Cholesky factorization.

Another approach, called Perturbation-Optimization [28, 29] is based on the following property:

If we note $\mathit{x}=\mathit{f}+{\left[{\mathit{H}}^{\prime }\mathit{H}+\mathit{\lambda }\mathit{I}\right]}^{-1}{\mathit{H}}^{\prime }\left(\mathit{g}-\mathit{H}\mathit{f}\right)$ and look for its expected and covariance matrix, it can be shown that:

$$\{\begin{array}{l}\mathrm{E}\left\{\mathit{x}\right\}=\widehat{\mathit{f}}\\ \mathrm{Cov}\left[\mathit{x}\right]=\widehat{\mathrm{\Sigma }}\end{array}$$E90

So, to generate a sample from the posterior law, we can do the following:

• Generate two random vectors ${\mathit{ϵ}}_{\mathit{f}}\sim \mathcal{N}\left({\mathit{ϵ}}_{\mathit{f}}|0,v_{\mathit{f}}I\right)$ and ${\mathit{ϵ}}_{\mathit{g}}\sim \mathcal{N}\left({\mathit{ϵ}}_{\mathit{g}}|0,v_{\mathit{ϵ}}I\right);$

• Define $\stackrel{\sim }{\mathit{g}}=\mathit{g}+{\mathit{ϵ}}_{\mathit{g}}$ and $\stackrel{\sim }{\mathit{f}}=\mathit{f}+{\mathit{ϵ}}_{\mathit{f}}$ and optimize

$$J\left(\stackrel{\sim }{\mathit{f}}\right)=\frac{1}{2}{‖\stackrel{\sim }{\mathit{g}}-\mathit{H}\mathit{f}‖}^2+\lambda {‖\stackrel{\sim }{\mathit{f}}-{\mathit{f}}_0‖}_2^2$$E91

• The obtained solution ${\mathit{f}}^{\left(n\right)}=arg{min}_{\stackrel{\sim }{\mathit{f}}}\left\{J\left(\stackrel{\sim }{\mathit{f}}\right)\right\}$ is a sample from the desired posterior law.

By repeating this process for a great number of times, we can use them to obtain good approximations for the posterior mean $\widehat{\mathit{f}}$ and the posterior covariance $\widehat{\mathrm{\Sigma }}$ by computing their empirical mean values. We need however fast and accurate optimization algorithms.

10. References to examples of applications

The above mentioned methods have been used with success in different applications:

• Medical imaging and Computed tomography (CT) [30, 31, 32, 33, 34, 35, 36].

• Diffraction tomography and Microwave imaging [37, 38, 39, 40, 41, 42].

• 3D Computed Tomography [18, 22, 43]

• Acoustical imaging [44, 45, 46, 47, 48, 49]

• Hyperspectral imaging [50]

• Spectrometry [51]

• Eddy current tomography [52]

• Non destructive testing applications [53]

• Emission Tomography [54]

• SAR imaging [55]

• Chronobiological time series [56]

11. Conclusions

Mainly, in this chapter, first we described inverse problems and gave a few classical examples such as deconvolution, image restoration, computed tomography X-ray image reconstruction, Fourier synthesis inversion problem which arise in many imaging systems. Then, we mentioned that there are two classes of methods for inverse problems: deterministic regularization and Bayesian inference methods. Then, we started by describing the Bayesian parameter estimation. The main parts of the chapter is focused on Bayesian inference for inverse problems. We saw that the main difficulty is the great dimension of unknown quantities and the appropriate choice of the prior law. For this, first we described many simple and hierarchical prior models which are used in real applications. For the second main difficulty, which is the computational aspects, we described different approximate Bayesian computations (ABC) and in particular the variational Bayesian approximation (VBA) methods and showed how to use these methods, for example for hyperparameter estimation or for large scale inverse problems.