#
Entropy, Information Theory, Information Geometry and Bayesian Inference in Data, Signal and Image Processing and Inverse Problems^{ †}

^{†}

Next Article in Journal

Next Article in Special Issue

Next Article in Special Issue

Previous Article in Journal / Special Issue

Laboratoire des Signaux et Systèmes, UMR 8506 CNRS-SUPELEC-UNIV PARIS SUD, SUPELEC, Plateau de Moulon, 3 rue Juliot-Curie, 91192 Gif-sur-Yvette, France

This paper is an extended version of the paper published in Proceedings of the 34th International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering (MaxEnt 2014), Amboise, France, 21–26 September 2014.

Received: 20 November 2014
/
Revised: 4 May 2015
/
Accepted: 5 May 2015
/
Published: 12 June 2015

(This article belongs to the Special Issue Information, Entropy and Their Geometric Structures)

The main content of this review article is first to review the main inference tools using Bayes rule, the maximum entropy principle (MEP), information theory, relative entropy and the Kullback–Leibler (KL) divergence, Fisher information and its corresponding geometries. For each of these tools, the precise context of their use is described. The second part of the paper is focused on the ways these tools have been used in data, signal and image processing and in the inverse problems, which arise in different physical sciences and engineering applications. A few examples of the applications are described: entropy in independent components analysis (ICA) and in blind source separation, Fisher information in data model selection, different maximum entropy-based methods in time series spectral estimation and in linear inverse problems and, finally, the Bayesian inference for general inverse problems. Some original materials concerning the approximate Bayesian computation (ABC) and, in particular, the variational Bayesian approximation (VBA) methods are also presented. VBA is used for proposing an alternative Bayesian computational tool to the classical Markov chain Monte Carlo (MCMC) methods. We will also see that VBA englobes joint maximum a posteriori (MAP), as well as the different expectation-maximization (EM) algorithms as particular cases.

As this paper is an overview and an extension of my tutorial paper in MaxEnt 2014 workshop [1], this Introduction gives a summary of the content of this paper.

The qualification Bayesian refers to the influence of Thomas Bayes [2], who introduced what is now known as Bayes’ rule, even if the idea had been developed independently by Pierre-Simon de Laplace [3]. For this reason, I am asking a question of the community if we shall change Bayes to Laplace and Bayesian to Laplacian or at least mention them both. Whatever the answer, we assume that the reader knows what probability means in a Bayesian or Laplacian framework. The main idea is that a probability law P(X) assigned to a quantity X represents our state of knowledge that we have about it. If X is a discrete valued variable, {P(X = x_{n}) = p_{n}; n = 1, ⋯, N} with mutually exclusive values x_{n} is its probability distribution. When X is a continuous valued variable, p(x) is its probability density function from which we can compute
$P(a\le X<b)={\displaystyle {\int}_{a}^{b}p(x)\phantom{\rule{0.2em}{0ex}}\mathrm{d}x}$ or any other probabilistic quantity, such as its mode, mean, median, region of high probabilities, etc.

In science, it happens very often that a quantity cannot be directly observed or, even when it can be observed, the observations are uncertain (commonly said to be noisy), by uncertain or noisy, here, I mean that, if we repeat the experiences with the same practical conditions, we obtain different data. However, in the Bayesian approach, for a given experiment, we have to use the data as they are, and we want to infer it from those observations. Before starting the observation and gathering new data, we have very incomplete knowledge about it. However, this incomplete knowledge can be translated in probability theory via an a priori probability law. We will discuss this point later on regarding how to do this. For now, we assume that this can be done. When a new observation (data D) on X becomes available (direct or indirect), we gain some knowledge via the likelihood P(D|X). Then, our state of knowledge is updated combining P(D|X) and P(X) to obtain an a posteriori law P(X|D), which represents the new state of knowledge on X. This is the main esprit of the Bayes rule, which can be summarized as:
As P(X|D) has to be a probability law, we have:
This relation can be extended to the continuous case. Some more details will be given in Section 2.

$$P(X|D)=P(D|X)P(X)/P(D).$$

$$P(D)={\displaystyle \sum _{X}P(D|X)P(X).}$$

Associated with a probability law is the quantity of information it contains. Shannon [4] introduced the notion of quantity of information I_{n} associated with one of the possible values of x_{n} of X with probabilities P(X = x_{n}) = p_{n} to be
${I}_{n}=\mathrm{ln}\frac{1}{{p}_{n}}=-\mathrm{ln}{p}_{n}$ and the entropy H as the expected value

$$H=-{\displaystyle \sum _{n=1}^{N}{p}_{n}\mathrm{ln}{p}_{n}.}$$

The word entropy has also its roots in thermodynamics and physics. However, this notion of entropy has no direct link with entropy in physics, even if for a particular physical system, we may attribute a probability law to a quantity of interest of that system and then define its entropy. This information definition of Shannon entropy became the main basis of information theory in many data analyses and the science of communication. More details and extensions about this subject will be given in Section 3.

As we can see up to now, we did not yet discuss how to assign a probability law to a quantity. For the discrete value variable, when X can take one of the N values {x_{1}, ⋯, x_{N}} and when we do not know anything else about it, Laplace proposed the “Principe d’indifférence”, where
$P(X={x}_{n})={p}_{n}=\frac{1}{N},\forall n=1,\cdots ,N$, a uniform distribution. However, what if we know more, but not enough to be able to assign the probability law {p_{1}, ⋯, p_{N}} completely?

For example, if we know that the expected value is
${\sum}_{n}{x}_{n}{p}_{n}=d$, this problem can be handled by considering this equation as a constraint on the probability distribution {p_{1}, ⋯, p_{N}}. If we have a sufficient number of constraints (at least N), then we may obtain a unique solution. However, very often, this is not the case. The question now is how to assign a probability distribution {p_{1}, ⋯, p_{N}} that satisfies the available constraints. This question is an ill-posed problem in the mathematical sense of Hadamard [5] in the sense that the solution is not unique. We can propose many probability distributions that satisfy the constraint imposed by this problem. To answer this question, Jaynes [6–8] introduced the maximum entropy principle (MEP) as a tool for assigning a probability law to a quantity on which we have some incomplete or macroscopic (expected values) information. Some more details about this MEP, the mathematical optimization problem, the expression of the solution and the algorithm to compute it will be given in Sections 3 and 4.

Kullback [9] was interested in comparing two probability laws and introduced a tool to measure the increase of information content of a new probability law with respect to a reference one. This tool is called either the Kullback–Leibler (KL) divergence, cross entropy or relative entropy. It has also been used to update a prior law when new pieces of information in the form of expected values are given. As we will see later, this tool can also be used as an extension to the MEP of Jaynes. Furthermore, as we will see later, this criterion of comparison of two probability laws is not symmetric: one of the probability laws has to be chosen to be the reference, and then, the second is compared to this reference. Some more details and extensions will be given in Section 5.

Fisher [10] wanted to measure the amount of information that a random variable X carries about an unknown parameter θ upon which its probability law p(x|θ) depends. The partial derivative with respect to θ of the logarithm of this probability law, called the log-likelihood function for θ, is called the score. He showed that the first order moment of the score is zero, but its second order moment is positive and is also equivalent to the expected values of the second derivative of log-likelihood function with respect to θ. This quantity is called the Fisher information. It is also been shown that for the small variations of θ, the Fisher information induces locally a distance in the space of parameters Θ, if we had to compare two very close values of θ. In this way, the notion of the geometry of information is introduced [11,12]. We must be careful here that this geometrical property is related to the space of the parameters Θ for small changes of the parameter for a given family of parametric probability law p(x|θ) and not in the space of probabilities. However, for two probability laws p_{1}(x) = p(x|θ_{1}) and p_{2}(x) = p(x|θ_{2}) in the same exponential family, the Kullback–Leibler divergence KL [p_{1} : p_{2}] induces a Bregman divergence B[θ_{1} : θ_{2}] between the two parameters [13,14]. More details will be given in Section 8.

At this stage, we have almost introduced all of the necessary tools that we can use for different levels of data, signal and image processing. In the following, we give some more details for each of these tools and their inter-relations. Then, we review a few examples of their use in different applications. As examples, we demonstrate how these tools can be used in independent components analysis (ICA) and source separation, data model selection, in spectral analysis of the signals and in inverse problems, which arise in many sciences and engineering applications. At the end, we focus more on the Bayesian approach for inverse problems. We present some details concerning unsupervised methods, where the hyper parameters of the problem have to be estimated jointly with the unknown quantities (hidden variables). Here, we will see how the Kullback–Leibler divergence can help approximate Bayesian computation (ABC). In particular, some original materials concerning variational Bayesian approximation (VBA) methods are presented.

Let us introduce things very simply. If we have two discrete valued related variables X and Y, for which we have assigned probability laws P(X) and P(Y), respectively, and their joint probability law P(X, Y), then from the sum and product rule, we have:
where P(X, Y) is the joint probability law,
$P(X)={\displaystyle {\sum}_{Y}P(X,Y)}$ and
$P(Y)={\displaystyle {\sum}_{X}P(X,Y)}$ are the marginals and
$P(X|Y)=P(X,Y)/P(Y)$ and
$P(Y|X)=P(X,Y)/P(X)$ are the conditionals. Now, consider the situation where Y can be observed, but not X. Because these two quantities are related, we may want to infer X from the observations on Y. Then, we can use:
which is called the Bayes rule.

$$P(X,Y)=P(X|Y)P(Y)=P(Y|X)P(X)$$

$$P(X|Y)=\frac{P(Y|X)P(X)}{P(Y)}$$

This relation is extended to the continuous valued variables using the measure theory [15,16]:
with:
More simply, the Bayes rule is often written as:
This writing can be used when we want to use p(x|y) to compute quantities that are only dependent on the shape of p(x|y), such as the mode, the median or quantiles. However, we must be careful that the denominator is of importance in many other cases, for example when we want to compute expected values. There is no need for more sophisticated mathematics here if we want to use this approach.

$$p(x|y)=\frac{p(y|x)p(x)}{p(y)}$$

$$p(y)={\displaystyle \int p(y|x)p(x)\phantom{\rule{0.2em}{0ex}}\mathrm{d}x}.$$

$$p(x|y)\phantom{\rule{0.2em}{0ex}}\propto \phantom{\rule{0.2em}{0ex}}p(y|x)p(x).$$

As we mentioned, the main use of this rule is in particular when X can not be observed (unknown quantity), but Y is observed and we want to infer X. In this case, the term p(y|x) is called the likelihood (of unknown quantity X in the observed data y), p(x) is called a priori and p(x|y) a posteriori. The likelihood is assigned using the link between the observed Y and the unknown X, and p(x) is assigned using the prior knowledge about it. The Bayes rule then is a way to do state of knowledge fusion. Before taking into account any observation, our state of knowledge is represented by p(x), and after the observation of Y, it becomes p(x|y).

However, in this approach, two steps are very important. The first step is the assigning of p(x) and p(y|x) before being able to use the Bayes rule. As noted in the Introduction and as we will see later, we need other tools for this step. The second important step is after: how to use p(x|y) to summarize it. When X is just a scalar value variable, we can do this computation easily. For example, we can compute the probability that X is in the interval [a, b] via:
However, when the unknown becomes a high dimensional vectorial variable **X**, as is the case in many signal and image processing applications, this computation becomes very costly [17]. We may then want to summarize p(x|y) by a few interesting or significant point estimates. For example, compute the maximum a posteriori (MAP) solution:
the expected a posteriori (EAP) solution:
the domains of X which include an integrated probability mass of more than some desired value (0.95 for example):
or any other questions, such as median or any α-quantiles:
As we see, computation of MAP needs an optimization algorithm, while these last three cases need integration, which may become very complicated for high dimensional cases [17].

$$P(a\le X<b|y)={\displaystyle {\int}_{a}^{b}p(x|y)}\phantom{\rule{0.2em}{0ex}}\mathrm{d}x.$$

$${\widehat{x}}_{\text{MAP}}=\underset{x}{\mathrm{arg}\mathrm{max}}\phantom{\rule{0.2em}{0ex}}\{p(x|y)\},$$

$${\widehat{x}}_{\text{EAP}}={\displaystyle \int x\phantom{\rule{0.2em}{0ex}}p(x|y)\phantom{\rule{0.2em}{0ex}}\mathrm{d}x},$$

$$\left[{x}_{1},{x}_{2}\right]:{\displaystyle {\int}_{{x}_{1}}^{{x}_{2}}p(x|y)\phantom{\rule{0.2em}{0ex}}\mathrm{d}x=.95},$$

$${x}_{q}:{\displaystyle {\int}_{-\infty}^{{x}_{q}}p(x|y)\phantom{\rule{0.2em}{0ex}}\mathrm{d}x=(1-\alpha )}.$$

We can also just explore numerically the whole space of the distribution using the Markov chain Monte Carlo (MCMC) [18–26] or any other sampling techniques [17]. In the scalar case (one dimension), all of these computations can be done numerically very easily. For the vectorial case, when the dimensions become large, we need to develop specialized approximation methods, such as VBA and algorithms to do these computations. We give some more details about these when using this approach for inverse problems in real applications.

$$\mathrm{E}\{X\}={E}_{p}\{X\}=<X>=<X{>}_{p}={\displaystyle \int x\phantom{\rule{0.2em}{0ex}}p(x)\phantom{\rule{0.2em}{0ex}}\mathrm{d}x}$$

$$\mathrm{E}\{h(X)\}={E}_{p}\{h(X)\}=<h(X)>=<h(X){>}_{p}={\displaystyle \int h(x)\phantom{\rule{0.2em}{0ex}}p(x)\phantom{\rule{0.2em}{0ex}}\mathrm{d}x}.$$

$$H[p]=\mathrm{E}\{-\mathrm{ln}(p(X))\}={\mathrm{E}}_{p}\{-\mathrm{ln}p(X)\}=<-\mathrm{ln}p(X)>=<-\mathrm{ln}p(X){>}_{p}=-{\displaystyle \int p(x)\phantom{\rule{0.2em}{0ex}}\mathrm{ln}p(x)\phantom{\rule{0.2em}{0ex}}\mathrm{d}x}.$$

$$\mathrm{E}\{X|y\}=\mathrm{E}{}_{p(x|y)}\{X\}=<X|y>=<X{>}_{p(x|y)}={\displaystyle \int x\phantom{\rule{0.2em}{0ex}}p(x|y)\phantom{\rule{0.2em}{0ex}}\mathrm{d}x}$$

$$\mathrm{E}\{h(X)|y\}={\mathrm{E}}_{p(x|y)}\{h(X)\}=<h(X)|y>=<h(X){>}_{p(x|y)}={\displaystyle \int h(x)\phantom{\rule{0.2em}{0ex}}p(x|)\phantom{\rule{0.2em}{0ex}}\mathrm{d}x}.$$

$$D[p|q]={\mathrm{E}}_{p}\left\{-\mathrm{ln}\frac{p(X)}{q(X)}\right\}=<-\mathrm{ln}\frac{p(X)}{q(X)}{>}_{p(x)}-{\displaystyle \int p(x)\mathrm{ln}\frac{p(x)}{q(x)}\mathrm{d}x}$$

$$D[p|q]={\mathrm{E}}_{p}\left\{-\mathrm{ln}\frac{p}{q}\right\}=<-\mathrm{ln}\frac{p}{q}{>}_{p}=-{\displaystyle \int p\mathrm{ln}\frac{p}{q}}.$$

$$\mathrm{E}\{h(X)|y\}={\mathrm{E}}_{p(x|y)}\{h(X)\}={\mathrm{E}}_{X|Y}\{h(X)\}=<h(X)|y>=<h(X){>}_{p(x|y)}={\displaystyle \int h(x)\phantom{\rule{0.2em}{0ex}}p(x|y)\phantom{\rule{0.2em}{0ex}}\mathrm{d}x}.$$

To introduce the quantity of information and the entropy, Shannon first considered a discrete valued variable X taking values {x_{1}, ⋯, x_{N}} with probabilities {p_{1}, ⋯, p_{N}} and defined the quantities of information associated with each of them as
${I}_{n}=\mathrm{ln}\frac{1}{{p}_{n}}=-\mathrm{ln}{p}_{n}$ and its expected value as the entropy:

$$\mathrm{H}[X]=-{\displaystyle \sum _{i=1}^{N}{p}_{i}\mathrm{ln}{p}_{i}}.$$

Later, this definition is extended to the continuous case by:
By extension, if we consider two related variables (X, Y) with the probability laws, joint p(x, y), marginals, p(x), p(y), and conditionals, p(y|x), p(x|y), we can define, respectively, the joint entropy:
as well as H[X], H[Y], H[Y|X] and H[X|Y].

$$\mathrm{H}[X]=-{\displaystyle \int p(x)\mathrm{ln}p(x)\phantom{\rule{0.2em}{0ex}}\mathrm{d}x}.$$

$$\mathrm{H}[X,Y]=-{\displaystyle \iint p(x,y)\mathrm{ln}p(x,y)\phantom{\rule{0.2em}{0ex}}\mathrm{d}x\phantom{\rule{0.2em}{0ex}}\mathrm{d}y},$$

Therefore, for any well-defined probability law, we can have an expression for its entropy. H[X], H[Y], H[Y|X], H[X|Y] and H[X, Y], which should better be noted as H[p(x)], H[p(y)], H[p(y|x)], H[p(x|y)] and H[p(x, y)].

Entropy is also a property of thermodynamical systems introduced by Clausius [27]. For a closed homogeneous system with reversible transformation, the differential entropy δS is related to δQ the incremental reversible transfer of heat energy into that system by δS = δQ/T with T being the uniform temperature of the closed system.

It is very hard to establish a direct link between these two entropies. However, in statistical mechanics, thanks to Boltzmann, Gibbs and many others, we can establish some link if we consider the microstates (for example, the number, positions and speeds of the particles) and the macrostates (for example, the temperature T, pressure P, volume V and energy E) of the system and if we assign a probability law to microstates and consider the macrostates as the average (expected values) of some functions of those microstates. Let us give a very brief summary of some of those interpretations.

The interpretation of entropy in statistical mechanics is the measure of uncertainty that remains about the state of a system after its observable macroscopic properties, such as temperature (T), pressure (P) and volume (V), have been taken into account. For a given set of macroscopic variables T, P and V, the entropy measures the degree to which the probability of the system is spread out over different possible microstates. In contrast to the macrostate, which characterizes plainly observable average quantities, a microstate specifies all atomic details about the system, including the position and velocity of every atom. Entropy in statistical mechanics is a measure of the number of ways in which the microstates of the system may be arranged, often taken to be a measure of “disorder” (the higher the entropy, the higher the disorder). This definition describes the entropy as being proportional to the natural logarithm of the number of possible microscopic configurations of the system (microstates), which could give rise to the observed macroscopic state (macrostate) of the system. The proportionality constant is the Boltzmann constant.

Boltzmann described the entropy as a measure of the number of possible microscopic configurations Ω of the individual atoms and molecules of the system (microstates) that comply with the macroscopic state (macrostate) of the system. Boltzmann then went on to show that k ln Ω was equal to the thermodynamic entropy. The factor k has since been known as Boltzmann’s constant.

In particular, Boltzmann showed that the entropy S of an ideal gas is related to the number of states of the molecules (microstates Ω) with a given temperature (macrostate):

$$S=k\mathrm{ln}\mathrm{\Omega}$$

The macroscopic state of the system is defined by a distribution on the microstates that are accessible to a system in the course of its thermal fluctuations. Therefore, the entropy is defined over two different levels of description of the given system. The entropy is given by the Gibbs entropy formula, named after J. Willard Gibbs. For a classical system (i.e., a collection of classical particles) with a discrete set of microstates, if E_{i} is the energy of microstate i and p_{i} is its probability that it occurs during the system’s fluctuations, then the entropy of the system is:
where k is again the physical constant of Boltzmann, which, like the entropy, has units of heat capacity. The logarithm is dimensionless. It is interesting to note that Relation (17) can be obtained from Relation (18) when the probability distribution is uniform over the volume Ω [28–30].

$$S=-k{\displaystyle \sum _{i=1}^{N}{p}_{i}\mathrm{ln}{p}_{i}}.$$

Kullback wanted to compare the relative quantity of information between two probability laws p_{1} and p_{2} on the same variable X. Two related notions have been defined:

- Relative Entropy of p
_{1}with respect to p_{2}:$$D[{p}_{1}:{p}_{2}]=-{\displaystyle \int {p}_{1}(x)\mathrm{ln}\frac{{p}_{1}(x)}{{p}_{2}(x)}}\phantom{\rule{0.2em}{0ex}}\mathrm{d}x$$ - Kullback–Leibler divergence of p
_{1}with respect to p_{2}:$$\mathrm{KL}[{p}_{1}:{p}_{2}]=-D[{p}_{1}:{p}_{2}]={\displaystyle \int {p}_{1}(x)\mathrm{ln}\frac{{p}_{1}(x)}{{p}_{2}(x)}}\phantom{\rule{0.2em}{0ex}}\mathrm{d}x$$

- KL [q : p] ≥ 0,
- KL [q : p] = 0, if q = p and
- KL [q : p
_{0}] ≥ KL [q : p_{1}] + KL [p_{1}: p_{0}]. - KL [q : p] is invariant with respect to a scale change, but is not symmetric.
- A symmetric quantity can be defined as:$$\mathrm{J}[q,p]=\frac{1}{2}(\mathrm{KL}[q:p]+\mathrm{KL}[p:q]).$$

The purpose of mutual information is to compare two related variables Y and X. It can be defined as the expected amount of information that one gains about Xif we observe the value of Y, and vice versa. Mathematically, the mutual information between X and Y is defined as:
or equivalently as:
With this definition, we have the following properties:
and:
We may also remark on the following property:

$$\mathrm{I}[Y,X]=\mathrm{H}[X]-\mathrm{H}[X|Y]=\mathrm{H}[Y]-\mathrm{H}[Y|X]$$

$$\mathrm{I}[Y,X]=D[p(X,Y):p(X)p(Y)].$$

$$\mathrm{H}[X,Y]=\mathrm{H}[X]+\mathrm{H}[Y|X]=\mathrm{H}[Y]+\mathrm{H}[X|Y]=\mathrm{H}[X]+\mathrm{H}[Y]-\mathrm{I}[Y,X]$$

$$\begin{array}{l}\mathrm{I}[Y,X]={\mathrm{E}}_{X}\{D[p(Y|X)]:p(Y)\}\stackrel{\mathrm{\Delta}}{=}{\displaystyle \int D[p(y|x):p(y)]p(x)\phantom{\rule{0.2em}{0ex}}\mathrm{d}x}\\ \phantom{\rule{3em}{0ex}}={E}_{Y}\{D[p(X|Y)]:p(X)\}\stackrel{\mathrm{\Delta}}{=}{\displaystyle \int D[p(x|y):p(x)]p(y)\phantom{\rule{0.2em}{0ex}}\mathrm{d}y}.\end{array}$$

- I[Y, X] is a concave function of p(y) when p(x|y) is fixed and a convex function of p(x|y) when p(y) is fixed.
- I[Y, X] ≥ 0 with equality only if X and Y are independent.

The first step before applying any probability rules for inference is to assign a probability law to a quantity. Very often, the available knowledge on that quantity can be described mathematically as the constraints on the desired probability law. However, in general, those constraints are not enough to determine in a unique way that probability law. There may exist many solutions that satisfy those constraints. We need then a tool to select one.

Jaynes introduced the MEP [8], which can be summarized as follows: When we do not have enough constraints to determine a probability law that satisfies those constraints, we may select between them the one with maximum entropy.

Let us be now more precise. Let us assume that the available information on that quantity X is the form of:
where ϕ_{k} are any known functions. First, we assume that such probability laws exist by defining:
with ϕ_{0} = 1 and d_{0} = 1 for the normalization purpose. Then, the MEP is written as an optimization problem:
whose solution is given by:
where Z(**λ**), called the partition function, is given by:
$Z(\mathbf{\lambda})={\displaystyle \int \mathrm{exp}[-{\displaystyle {\sum}_{k=1}^{K}{\lambda}_{k}{\varphi}_{k}(x)}]\phantom{\rule{0.2em}{0ex}}\mathrm{d}x}$ and
$\mathbf{\lambda}={\left[{\lambda}_{1},\dots ,{\lambda}_{K}\right]}^{\prime}$ have to satisfy:
which can also be written as −∇_{λ} ln Z(**λ**) = **d**. Different algorithms have been proposed to compute numerically the ME distributions. See, for example, [31–37]

$$E\left\{{\varphi}_{k}(X)\right\}={d}_{k},\phantom{\rule{1em}{0ex}}k=1,\dots ,K.$$

$$\mathcal{P}=\left\{p(x):{\displaystyle \int {\varphi}_{k}(x)p(x)\phantom{\rule{0.2em}{0ex}}\mathrm{d}x=\mathrm{d}k,\phantom{\rule{1em}{0ex}}k=0,\dots ,K}\right\}$$

$${p}_{ME}(x)=\underset{p\in \mathcal{P}}{\mathrm{arg}\mathrm{max}}\left\{H[p]=-{\displaystyle \int p(x)\mathrm{ln}p(x)\phantom{\rule{0.2em}{0ex}}\mathrm{d}x}\right\}$$

$${p}_{ME}(x)=\frac{1}{Z(\mathbf{\lambda})}\mathrm{exp}\left[-{\displaystyle \sum _{k=1}^{K}{\lambda}_{k}{\varphi}_{k}(x)}\right]$$

$$-\frac{\partial \mathrm{ln}Z(\mathbf{\lambda})}{\partial {\lambda}_{k}}={d}_{k},\phantom{\rule{1em}{0ex}}k=1,\dots ,K$$

The maximum value of entropy reached is given by:
This optimization can easily be extended to the use of relative entropy by replacing H(p) by D[p : q] where q(x) is a given reference of a priori law. See [9,38,39] and [34,40–42] for more details.

$${H}_{\mathrm{max}}=\mathrm{ln}Z(\mathbf{\lambda})+{\mathbf{\lambda}}^{\prime}\mathit{d}.$$

Consider the problem of the parameter estimation **θ** of a probability law p(x|**θ**) from an n-element sample of data **x** = {x_{1}, ⋯, x_{n}}.

The log-likelihood of **θ** is defined as:
Maximizing L(**θ**) with respect to **θ** gives what is called the maximum likelihood (ML) estimate of θ.

$$L(\mathit{\theta})=\mathrm{ln}{\displaystyle \prod _{i=1}^{n}p({x}_{i}|\theta )}={\displaystyle \sum _{i=1}^{n}\mathrm{ln}p({x}_{i}|\theta )}.$$

Noting that L(**θ**) depends on n, we may consider
$\frac{1}{n}L(\mathbf{\theta})$ and define:
where **θ**^{*} is the right answer and p(x|**θ**^{*}) its corresponding probability law. We may then remark that:
The first term in the right-hand side being a constant, we derive that:
In this way, there is a link between the maximum likelihood and maximum relative entropy solutions [24].

$$\overline{L}(\mathit{\theta})=\underset{x\mapsto \infty}{\mathrm{lim}}\frac{1}{n}L(\mathit{\theta})=\mathrm{E}\{\mathrm{ln}p(x|\mathit{\theta})\}={\displaystyle \int p(x|\mathit{\theta}*)\mathrm{ln}p(x|\mathit{\theta})\phantom{\rule{0.2em}{0ex}}\mathrm{d}x},$$

$$D[p(x|\mathit{\theta}*):p(x|\mathit{\theta})]=-{\displaystyle \int p(x|\mathit{\theta}*)\mathrm{ln}\frac{p(x|\mathit{\theta})}{p(x|\mathit{\theta}*)}\phantom{\rule{0.2em}{0ex}}\mathrm{d}x}=-{\displaystyle \int p(x|\mathit{\theta}*)\mathrm{ln}p(x|\mathit{\theta}*)\phantom{\rule{0.2em}{0ex}}\mathrm{d}x+\overline{L}(\mathit{\theta})}.$$

$$\underset{\mathit{\theta}}{\mathrm{arg}\mathrm{max}}\{D[p(x|\mathit{\theta}*);p(x|\mathit{\theta})]\}=\underset{\mathit{\theta}}{\mathrm{arg}\mathrm{max}}\{\overline{L}(\mathit{\theta})\}.$$

Fisher [10] was interested in measuring the amount of information that samples of a variable X carries about an unknown parameter θ upon which its probability law p(x|θ) depends. For a given sample of observation x and its probability law p(x|θ), the function L(θ) = p(x|θ) is called the likelihood of θ in the sample x. He called the score of x over θ the partial derivative with respect to θ of the logarithm of this function:
He also showed that the first order moment of the score is zero:
but its second order moment is positive and is also equivalent to the expected values of the second derivative of the log-likelihood function with respect to θ.
This quantity is called the Fisher information [14].

$$S(x|\theta )=\frac{\partial \mathrm{ln}p(x|\theta )}{\partial \theta}$$

$$\mathrm{E}\{S(X|\theta )\}=\mathrm{E}\left\{\frac{\partial \mathrm{ln}p(x|\theta )}{\partial \theta}\right\}=0$$

$$\mathrm{E}\{{S}^{2}(X|\theta )\}=\mathrm{E}\left\{{\left|\frac{\partial \mathrm{ln}p(x|\theta )}{\partial \theta}\right|}^{2}\right\}=\mathrm{E}\left\{\frac{{\partial}^{2}\mathrm{ln}p(x|\theta )}{\partial {\theta}^{2}}\right\}=F$$

It is also shown that for the small variations of θ, the Fisher information induces locally a distance in the space of parameters Θ, if we had to compare two very close values of θ. In this way, the notion of the geometry of information is introduced. The main steps for introducing this notion are the following: Consider D [p(x|**θ**^{*}) : p(x|**θ**^{*} + ∆**θ**)] and assume that ln p(x|**θ**) can be developed in a Taylor series. Then, keeping the terms up to the second order, we obtain:
where **F** is the Fisher information:

$$D[p(x|\mathit{\theta}*):p(x|\mathit{\theta}*+\mathrm{\Delta}\mathit{\theta})]\simeq \frac{1}{2}\mathrm{\Delta}{\mathit{\theta}}^{\prime}\mathit{F}(\theta *)\mathrm{\Delta}\mathit{\theta}.$$

$$\mathit{F}(\theta *)=\mathrm{E}\left\{\frac{{\partial}^{2}\mathrm{ln}p(x|\mathit{\theta})}{\partial {\mathit{\theta}}^{\prime}\partial \mathit{\theta}}{|}_{\mathit{\theta}=\mathit{\theta}*}\right\}.$$

We must be careful here that this geometry property is related to the space of the parameters Θ for a given family of parametric probability law p(x|θ) and not in the space of probabilities. However, for two probability laws p_{1}(x) = p(x|θ_{1}) and p_{2}(x) = p(x|θ_{2}) in the same exponential family, the Kullback–Leibler divergence KL [p_{1} : p_{2}] induces a Bregman divergence B[θ_{1}|θ_{2}] between the two parameters [14,45–48].

To go further into detail, let us extend the discussion about the link between Fisher information and KL divergence, as well as other divergences, such as f-divergences, Rényi’s divergences and Bregman divergences.

- f-divergences:The f-divergences, which are a general class of divergences, indexed by convex functions f, that include the KL divergence as a special case. Let f: (0, ∞) ↦
**R**be a convex function for which f(1) = 0. The f-divergence between two probability measures P and Q is defined by:$${D}_{f}[P:Q]={\displaystyle \int q\phantom{\rule{0.2em}{0ex}}f\left(\frac{p}{q}\right)}$$Every f-divergence can be viewed as a measure of distance between probability measures with different properties. Some important special cases are:- – f(x) = x ln x gives KL divergence: $\mathrm{KL}[P:Q]={\displaystyle \int p\mathrm{ln}\left(\frac{p}{q}\right)}$.
- – $f(x)=|x-1|/2$ gives total variation distance: $\mathrm{TV}[P,Q]={\displaystyle \int |p-q|/2}$.
- – $f(x)={\left(\sqrt{x}-1\right)}^{2}$ gives the square of the Hellinger distance: ${H}^{2}[P,Q]={\displaystyle \int {\left(\sqrt{p}-\sqrt{q}\right)}^{2}}$.
- – f(x) = (x − 1)
^{2}gives the chi-squared divergence: ${\chi}^{2}[P:Q]={\displaystyle \int \frac{{(p-q)}^{2}}{q}}$.

- Rényi divergences:These are another generalization of the KL divergence. The Rényi divergence between two probability distributions P and Q is:$${D}_{\alpha}[P:Q]=\frac{1}{\alpha -1}\mathrm{ln}{\displaystyle \int {p}^{\alpha}{q}^{1-\alpha}}.$$
_{α}[P : Q] converges to KL [P : Q].${D}_{1/2}[P,Q]=-2\mathrm{ln}{\displaystyle \int \sqrt{pq}}$ is called Bhattacharyya divergence (closely related to Hellinger distance). Interestingly, this quantity is always smaller than KL:$${D}_{1/2}[P:Q]\le \mathrm{KL}[P:Q].$$As a result, it is sometimes easier to derive risk bounds with ${D}_{1/2}$ as the loss function as opposed to KL. - Bregman divergences:The Bregman divergences provide another class of divergences that are indexed by convex functions and include both the Euclidean distance and the KL divergence as special cases. Let ϕ be a differentiable strictly convex function. The Bregman divergence B
_{ϕ}is defined by:$${B}_{\varphi}[\mathbf{x}:\mathbf{y}]=\varphi (\mathbf{x})-\varphi (\mathbf{y})-\langle \mathbf{x}-\mathbf{y},\nabla \varphi (\mathbf{y})\rangle $$**x**and**y**and where the domain of ϕ is a space where convexity and differentiability make sense (e.g., whole or a subset of**R**^{d}or an L_{p}space). For example, $\varphi (\mathbf{x})={\Vert \mathbf{x}\Vert}^{2}$ and**R**^{d}gives the Euclidean distance:$${B}_{\varphi}[\mathbf{x}:\mathbf{y}]=\varphi (\mathbf{x})-\varphi (\mathbf{y})-\langle \mathbf{x}-\mathbf{y},\nabla \varphi (\mathbf{y})\rangle ={\Vert \mathbf{x}\Vert}^{2}-{\Vert \mathbf{y}\Vert}^{2}-\langle \mathbf{x}-\mathbf{y},2\mathbf{y}\rangle ={\Vert \mathbf{x}-\mathbf{y}\Vert}^{2}$$**R**^{d}gives the KL divergence:$${B}_{\varphi}[\mathbf{x}:\mathbf{y}]={\displaystyle \sum _{j}{x}_{j}\mathrm{ln}{x}_{j}}-{\displaystyle \sum _{j}{y}_{j}\mathrm{ln}{y}_{j}}-{\displaystyle \sum _{j}({x}_{j}-{y}_{j})(1+\mathrm{ln}{y}_{j})}={\displaystyle \sum _{j}{x}_{j}\mathrm{ln}\frac{{x}_{j}}{{y}_{j}}}=\mathrm{KL}[\mathbf{x}:\mathbf{y}]$$

Let X be a quantity taking values in the domain of ϕ with a probability distribution function p(x).

Then, E_{p}_{(}_{x}_{)} {B_{ϕ}(X, m)} is minimized over m in the domain of ϕ at m = E {X}:
Moreover, this property characterizes Bregman divergence. When applied to the Bayesian approach, this means that, using the Bregman divergence as the loss function, the Bayes estimator is the posterior mean. This point is detailed in the following.

$$\widehat{m}=\underset{m}{\mathrm{arg}\mathrm{max}}\{{B}_{\varphi}(X,m)\}=\mathrm{E}\{X\}.$$

Links between all of these through an example:

Let us consider the Bayesian parameter estimation where we have some data **y**, a set of parameters **x**, a likelihood p(y|x) and a prior π(x), which gives the posterior
$p(x|y)\propto p(y|x)\pi (x)$. Let us also consider a cost function
$C[x,\tilde{x}]$ in the parameter space x ∈ **X**. The classical Bayesian point estimation of **x** is expressed as the minimizer of an expected risk:
where:
It is very well known that the mean squared error estimator, which corresponds to
$C[\mathbf{x},\tilde{\mathbf{x}}]={\Vert \mathbf{x}-\tilde{\mathbf{x}}\Vert}^{2}$, is the posterior mean. It is now interesting to know that choosing
$C[\mathbf{x},\tilde{\mathbf{x}}]$ to be any Bregman divergence
${B}_{\varphi}[x,\tilde{x}]$, we obtain also the posterior mean:
where:

$$\widehat{\mathbf{x}}=\underset{\tilde{\mathbf{x}}}{\mathrm{arg}\mathrm{min}}\left\{\overline{C}(\tilde{\mathbf{x}})\right\}$$

$$\overline{C}(\tilde{\mathbf{x}})={\mathrm{E}}_{p(x|y)}\{C[\mathbf{x},\tilde{\mathbf{x}}]\}={\displaystyle \int C[\mathbf{x},\tilde{\mathbf{x}}]p(\mathbf{x}|\mathbf{y})\phantom{\rule{0.2em}{0ex}}\mathrm{d}x}$$

$$\widehat{\mathbf{x}}=\underset{\tilde{\mathbf{x}}}{\mathrm{arg}\mathrm{min}}\{{\overline{B}}_{\varphi}(\tilde{\mathbf{x}})\}={\mathrm{E}}_{p(\mathbf{x}|\mathbf{y})}\left\{{\displaystyle \int \mathbf{x}\phantom{\rule{0.2em}{0ex}}p(\mathbf{x}|\mathbf{y})\phantom{\rule{0.2em}{0ex}}\mathrm{d}\mathbf{x}}\right\}$$

$${\overline{B}}_{\varphi}(\tilde{\mathbf{x}})={\mathrm{E}}_{p(\mathbf{x}|\mathbf{y})}\{{D}_{\varphi}[\mathbf{x},\tilde{\mathbf{x}}]\}={\displaystyle \int {B}_{\varphi}[\mathbf{x},\tilde{\mathbf{x}}]p(\mathbf{x}|\mathbf{y})\phantom{\rule{0.2em}{0ex}}\mathrm{d}\mathbf{x}}$$

Consider now that we have two prior probability laws π_{1}(**x**) and π_{2}(**x**), which give rise to two posterior probability laws p_{1}(**x**|**y**) and p_{2}(**x**|**y**). If the prior laws and the likelihood are in the exponential families, then the posterior laws are also in the exponential family. Let us note them as p_{1}(**x**|**y**; **θ**_{1}) and p_{2}(**x**|**y**; **θ**_{2}), where **θ**_{1} and **θ**_{2} are the parameters of those posterior laws. We then have the following properties:

- KL [p
_{1}: p_{2}] is expressed as a Bregman divergence B[**θ**_{1}:**θ**_{2}]. - A Bregman divergence B[
**x**_{1}:**x**_{2}] is induced when KL [p_{1}: p_{2}] is used to compare the two posteriors.

The extension of the scalar variable to the finite dimensional vectorial case is almost immediate. In particular, for the Gaussian case
$p(\mathit{x})=\mathcal{N}(\mathit{x}|\mu ,\mathit{R})$, the mean becomes a vector μ = E{**X**}, and the variances are replaced by a covariance matrix:
$\mathit{R}=\mathrm{E}\left\{(\mathit{X}-\mu ){\left(\mathit{X}-\mu \right)}^{\prime}\right\}$; and almost all of the quantities can be defined immediately. For example, for a Gaussian vector
$p(\mathit{x})=\mathcal{N}(\mathit{x}|0,\mathit{R})$, the entropy is given by [49]:
and the relative entropy of
$\mathcal{N}(x|\mathbf{0},\mathit{R})$$\mathcal{N}(\mathit{x}|0,\mathit{R})$with respect to
$\mathcal{N}(\mathit{x}|0,\mathbf{S})$ is given by:

$$H=\frac{n}{2}\mathrm{ln}(2\mathrm{\pi})+\frac{1}{2}\mathrm{ln}\left(|\mathrm{det}\phantom{\rule{0.2em}{0ex}}(\mathit{R})|\right)$$

$$D=\frac{1}{2}\left(\mathrm{tr}(\mathit{R}{\mathit{S}}^{-1})-\mathrm{log}\frac{|\mathrm{det}(\mathit{R})|}{\mathrm{det}(\mathit{S})}-n\right).$$

The notion of time series or processes needs extra definitions. For example, for a random time series X(t), we can define p(X(t)), ∀t, the expected value time series
$\overline{x}(t)=\mathrm{E}\left\{X(t)\right\}$ and what is called the autocorrelation function Γ(t, τ) = E {X(t) X(t + τ)}. A time series is called stationary when these quantities does not depend on t, i.e.,
$\overline{x}(t)=m$ and Γ(t, τ) = Γ(τ) [50]. Another quantity of interest for a stationary time series is its power spectral density (PSD) function:

$$S(\omega )=\mathrm{FT}\{\Gamma (\tau )\}={\displaystyle \int \Gamma (\tau )\mathrm{exp}[-j\omega \tau ]\phantom{\rule{0.2em}{0ex}}\mathrm{d}\tau}.$$

When X(t) is observed on times t = nΔT with ΔT = 1, we have X(n), and for a sample {X(1),⋯, X(N)}, we may define the mean **μ** = E {**X**} and the covariance matrix
$\sum =\mathrm{E}\{(\mathit{X}-\mu ){(\mathit{X}-\mu )}^{\prime}\}$.

With these definitions, it can easily been shown that the covariance matrix of a stationary Gaussian process is Toeplitz [49]. It is also possible to show that the entropy of such a process can be expressed as a function of its PSD function:

$$\underset{n\to \infty}{\mathrm{lim}}\frac{1}{n}H(p)=\frac{1}{2\pi}{\displaystyle {\int}_{-\pi}^{\pi}\mathrm{ln}S(\omega )\phantom{\rule{0.2em}{0ex}}\mathrm{d}\omega}.$$

For two stationary Gaussian processes with two spectral density functions S_{1}(ω) and S_{2}(ω), we have:
where we find the Itakura–Saito distance in the spectral analysis literature [50–53].

$$\underset{n\to \infty}{\mathrm{lim}}\frac{1}{n}D[{p}_{1}:{p}_{2}]=\frac{1}{4\pi}{\displaystyle {\int}_{-\pi}^{\pi}\left(\frac{{S}_{1}(\omega )}{{S}_{2}(\omega )}-\mathrm{ln}\frac{{S}_{1}(\omega )}{{S}_{2}(\omega )}-1\right)}\phantom{\rule{0.2em}{0ex}}\mathrm{d}\omega $$

These definitions and expressions have often been used in time series analysis. In what follows, we give a few examples of the different ways these notions and quantities have been used in different applications of data, signal and image processing.

Given a vector of time series **x**(t), the independent component analysis (ICA) consists of finding a separating matrix **B**, such that the components y(t) = **Bx**(t) are as independent as possible. The notion of entropy is used here as a measure of independence. For example, to find **B**, we may choose
$D\left[p(y):{\displaystyle {\prod}_{j}{p}_{j}({y}_{j})}\right]$ as a criterion of independence of the components y_{j}. The next step is to choose a probability law p(**x**) from which we can find an expression for p(**y**) from which we can find an expression for
$D\left[p(y):{\displaystyle {\prod}_{j}{p}_{j}({y}_{j})}\right]$ as a function of the matrix **B**, which can be optimized to obtain it.

The ICA problem has a tight link with the source separation problem, where it is assumed that the measured time series **x**(t) is a linear combination of the sources s(t), i.e., **x**(t) = **As**(t), with **A** being the mixing matrix. The objective of source separation is then to find the separating matrix **B** = **A**^{−}^{1}.

To see how the entropy is used here, let us note **y** = **Bx**. Then,
H(**y**) is used as a criterion for ICA or source separation. As the objective in ICA is to obtain **y** in such a way that its components become as independent as possible, the separating matrix **B** has to maximize H(**y**). Many ICA algorithms are based on this optimization [54–65]

$${p}_{Y}(\mathbf{y})=\frac{1}{|\partial \mathbf{y}/\partial \mathit{x}|}{p}_{X}(\mathit{x})\to H(\mathbf{y})=-\mathrm{E}\{\mathrm{ln}{p}_{Y}(\mathbf{y})\}=\mathrm{E}\left\{\mathrm{ln}|\partial \mathbf{y}/\partial \mathrm{x}|\right\}-H(\mathit{x}).$$

Determining the order of a model, i.e., the dimension of the vector parameter **θ** in a probabilistic model p(**x**|**θ**), is an important subject in many data and signal processing problems. As an example, in autoregressive (AR) modeling:
where **θ** = {θ_{1}⋯, θ_{K}}, we may want to compare two models with two different values of K.

$$x(n)={\displaystyle \sum _{k=1}^{K}{\theta}_{k}x(n-k)+\epsilon (n)}$$

When the order K is fixed, the estimation of the parameters **θ** is a very well-known problem, and there are likelihood based [66] or Bayesian approaches for that [67]. The determination of the order is however more difficult [68]. Between the tools, we may mention here the Bayesian methods [69–74], but also the use of relative entropy D [p(**x**|**θ**^{*}): p(**x**|**θ**)], where **θ**^{*} represents the vector of the parameters of dimension K^{*} and **θ** and the vector **θ** with dimension K ≤ K^{*}. In such cases, even if the two probability laws to be compared have parameters with different dimensions, we can always use the KL [p(**x**|**θ**^{*}): p(**x**|**θ**)] to compare them. The famous criterion of Akaike [75–78] uses this quantity to determine the optimal order. For a linear parameter model with Gaussian probability laws and likelihood-based methods, there are analytic solutions for it [68].

Entropy and MEP have been used in different ways in the spectral analysis problem. It has been an important subject of signal processing for the decades. Here, we are presenting, in a brief way, these different approaches.

A classical one is Burg’s entropy method [79], which can be summarized as follows: Let X(n) be a stationary, centered process, and assume we have as data a finite number of samples (lags) of its autocorrelation function:
The task is then to estimate its power spectral density function:
As we can see, due to the fact that we have only the elements of the right-hand for k = −K, ⋯, +K, the problem is ill posed. To obtain a probabilistic solution, we may start by assigning a probability law p(**x**) to the vector
$\underset{\xaf}{X}=[X(0),\dots ,X(N-1){]}^{\prime}$. For this, we can use the principle of maximum entropy (PME) with the data as constraints (54). As these constraints are the second order moments, the PME solution is a Gaussian probability law:
$\mathcal{N}(\mathit{x}|0,\mathit{R})$. For a stationary Gaussian process, when the number of samples N → ∞, the expression of the entropy becomes:
This expression is called Burg’s entropy [79]. Thus, Burg’s method consists of maximizing H subject to the constraints (54). The solution is:
where
$\mathbf{\lambda}=[{\lambda}_{0},\cdots ,{\lambda}_{K}{]}^{\prime}$, the Lagrange multipliers associated with the constraints (54), are here equivalent to the AR modeling of the Gaussian process X(n).

$$r(k)=\mathrm{E}\{X(n)X(n+k)\}=\frac{1}{2\pi}{\displaystyle {\int}_{-\pi}^{\pi}S(\omega )\mathrm{exp}[jk\omega ]\phantom{\rule{0.2em}{0ex}}\mathrm{d}\omega},\phantom{\rule{1em}{0ex}}k=0,\dots ,K.$$

$$S(\omega )={\displaystyle \sum _{k=-\infty}^{\infty}r(k)\mathrm{exp}[-jk\omega ]}$$

$$H={\displaystyle {\int}_{-\pi}^{\pi}\mathrm{ln}S(\omega )\phantom{\rule{0.2em}{0ex}}\mathrm{d}\omega}.$$

$$S(\omega )=\frac{1}{{\left|{\displaystyle \sum _{k=-K}^{K}{\mathrm{\lambda}}_{k}\mathrm{exp}[\mathrm{j}k\omega ]}\right|}^{2}},$$

We may note that, in this particular case, we have an analytical expression for **λ**, which provides the possibility to give an analytical expression for S(ω) as a function of the data {r(k), k = 0,⋯, K}:
where **Γ** = Toeplitz(r(0),⋯, r(K)) is the correlation matrix and **δ** and **e** are two vectors defined by **δ** = [1, 0,⋯, 0]′ and **e** = [1, e^{−}^{j}^{ω}, e^{−}^{j2}^{ω},⋯, e^{−}^{j}^{Kω}]′.

$$S(\omega )=\frac{\mathit{\delta}{\Gamma}^{-1}\mathit{\delta}}{e{\Gamma}^{-1}e},$$

We may note that we first used MEP to choose a probability law for X(n). With the prior knowledge that we have second order moments, the MEP results in a Gaussian probability density function. Then, as for a stationary Gaussian process, the expression of the entropy is related to the power spectral density S(ω), and as this is related to the correlation data by a Fourier transform, an ME solution could be computed easily.

The second approach consists of maximizing the relative entropy D [p(**x**): p_{0}(**x**)] or minimizing KL [p(**x**) : p_{0}(**x**)] where p_{0}(**x**) is an a priori law. The choice of the prior is important. Choosing a uniform p_{0}(**x**), we retrieve the previous case [77].

However, choosing a Gaussian law for p_{0}(**x**), the expression to maximize becomes:
when N ↦ ∞ and where S_{0}(ω) corresponds to the power spectral density of the reference process p_{0}(**x**). Now, the problem becomes: minimize D [p(**x**): p_{0}(**x**)] subject to the constraints (54).

$$D[p(\mathit{x}):{p}_{0}(\mathit{x})]=\frac{1}{4\pi}{\displaystyle {\int}_{-\pi}^{\pi}\left(\frac{S(\omega )}{{S}_{0}(\omega )}-\mathrm{ln}\frac{S(\omega )}{{S}_{0}(\omega )}-1\right)}\phantom{\rule{0.2em}{0ex}}\mathrm{d}\omega $$

Another approach is to decompose first the process X(n) on the Fourier basis {cos kωt, sin kωt}, to consider ω to be the variable of interest and S(ω), normalized properly, to be considered as its probability distribution function. Then, the problem can be reformulated as the determination of the S(ω), which maximizes the entropy:
subject to the linear constraints (54). The solution is in the form of:
which can be considered as the most uniform power spectral density that satisfies those constraints.

$$-{\displaystyle {\int}_{-\pi}^{\pi}S(\omega )\mathrm{ln}S(\omega )\phantom{\rule{0.2em}{0ex}}\mathrm{d}\omega}$$

$$S(\omega )=\mathrm{exp}\left[{\displaystyle \sum _{k=-K}^{K}{\mathrm{\lambda}}_{k}\mathrm{exp}[\mathrm{j}k\omega ]}\right].$$

In this approach, we consider S(ω) as the expected value Z(ω) for which we have a prior law μ(z), and we are looking to assign p(z), which maximizes the relative entropy D [p(z) : μ(z)] subject to the constraints (54).

When p(z) is determined, the solution is given by:
The expression of S(ω) depends on μ(z). When μ(z) is Gaussian, we obtain the Rényi entropy:
If we choose a Poisson measure for μ(z), we obtain the Shannon entropy:
and if we choose a Lebesgue measure over [0, ∞], we obtain Burg’s entropy:
$$H={\displaystyle {\int}_{-\pi}^{\pi}\mathrm{ln}S(\omega )\phantom{\rule{0.2em}{0ex}}\mathrm{d}\omega}.$$
When this step is done, the next step becomes maximizing these entropies subject to the constraints of the correlations. The obtained solutions are very different. For more details, see [39,79–85].

$$S(\omega )=\mathrm{E}\{Z(\omega )\}={\displaystyle \int Z(\omega )p(z)\phantom{\rule{0.2em}{0ex}}\mathrm{d}z}.$$

$$H={\displaystyle {\int}_{-\pi}^{\pi}{S}^{2}(\omega )\phantom{\rule{0.2em}{0ex}}\mathrm{d}\omega}.$$

$$H=-{\displaystyle {\int}_{-\pi}^{\pi}S(\omega )\mathrm{ln}S(\omega )\phantom{\rule{0.2em}{0ex}}\mathrm{d}\omega},$$

A general way to introduce inverse problems is the following: Infer an unknown signal f(t), image f(x, y) or any multi-variable function f(**r**) through an observed signal g(t), image g(x, y) or any multi-variable observable function g(**s**), which are related through an operator
$\mathscr{H}:f\mapsto g$. This operator can be linear or nonlinear. Here, we consider only linear operators
$g=Hf$:
where h(**r**, **s**) is the response of the measurement system. Such linear operators are very common in many applications of signal and image processing. We may mention a few examples of them:

$$g(\mathbf{s})={\displaystyle \int h(\mathbf{r},\mathbf{s})f(\mathbf{r})\phantom{\rule{0.2em}{0ex}}\mathrm{d}\mathbf{r}}$$

- Convolution operations g = h * f in 1D (signal):$$g(t)={\displaystyle \int h(t-{t}^{\prime})f({t}^{\prime})\phantom{\rule{0.2em}{0ex}}\mathrm{d}{t}^{\prime}}$$$$g(x,y)={\displaystyle \iint h(x-{x}^{\prime},y-{y}^{\prime})f({x}^{\prime},{y}^{\prime})\phantom{\rule{0.2em}{0ex}}\mathrm{d}{x}^{\prime}\phantom{\rule{0.2em}{0ex}}\mathrm{d}{y}^{\prime}}$$
- Radon transform (RT) in computed tomography (CT) in the 2D case [86]:$$g(r,\varphi )={\displaystyle \iint \delta (r-x\mathrm{cos}\varphi -y\mathrm{sin}\varphi )f(x,y)\phantom{\rule{0.2em}{0ex}}\mathrm{d}x\phantom{\rule{0.2em}{0ex}}\mathrm{d}y}$$
- Fourier transform (FT) in the 2D case:$$g(u,v)={\displaystyle \iint \mathrm{exp}[-j(ux+uy)]f(x,y)\phantom{\rule{0.2em}{0ex}}\mathrm{d}x\phantom{\rule{0.2em}{0ex}}\mathrm{d}y}$$

No matter the category of the linear transforms, when the problem is discretized, we arrive at the relation:
where **f** = [f_{1},⋯, f_{n}]′ represents the unknowns, **g** = [g_{1},⋯, g_{m}]′ the observed data, **ϵ** = [ϵ_{1},⋯, ϵ_{m}]′ the errors of modeling and measurement and **H** the matrix of the system response.

$$\mathit{g}=\mathit{Hf}+\mathit{\u03f5}.$$

Let us consider first the simple no noise case:
where **H** is a matrix of dimensions (M × N), which is in general singular or very ill conditioned. Even if the cases M > N or M = N may appear easier, they have the same difficulties as those of the underdetermined case M < N that we consider here. In this case, evidently the problem has an infinite number of solutions, and we need to choose one.

$$\mathit{g}=\mathit{Hf},$$

Between the numerous methods, we may mention the minimum norm solution, which consists of choosing between all of the possible solutions:
the one that has the minimum norm:
This optimization problem can be solved easily in this case, and we obtain:
In fact, we may choose any other convex criterion Ω(**f**) and satisfy the uniqueness of the solution. For example:
which can be interpreted as the entropy when f_{j} > 0 and ∑ f_{j} = 1, thus considering f_{j} as a probability distribution f_{j} = P (U = u_{j}). The variable U can correspond (or not) to a physical quantity. Ω(**f**) is the entropy associated with this variable.

$$\mathrm{F}=\{\mathit{f}:\mathit{Hf}=\mathit{g}\}$$

$$\mathrm{\Omega}(\mathit{f})={\Vert \mathit{f}\Vert}_{2}^{2}={\displaystyle \sum _{j}{f}_{j}^{2}}.$$

$${\widehat{f}}_{NM}=\underset{\mathit{f}\in \mathcal{F}}{\mathrm{arg}\mathrm{min}}\{\mathrm{\Omega}(\mathit{f})={\Vert \mathit{f}\Vert}_{2}^{2}\}={\mathit{H}}^{\prime}{(\mathit{H}{\mathit{H}}^{\prime})}^{-1}\mathit{g}.$$

$$\mathrm{\Omega}(\mathit{f})=-{\displaystyle \sum _{j}{f}_{j}\mathrm{ln}{f}_{j}}$$

If we consider f_{j} > 0 to represent the power spectral density of a physical quantity, then the entropy becomes:
and we can use it as criterion to select a solution to the problem (71).

$$\mathrm{\Omega}(\mathit{f})={\displaystyle \sum _{j}{f}_{j}\mathrm{ln}{f}_{j}}$$

As we can see, any convex criterion Ω(**f**) can be used. Here, we mentioned four of them with different interpretations.

- L
_{2}or quadratic:$$\mathrm{\Omega}(\mathit{f})={\displaystyle \sum _{j}{f}_{j}^{2}}$$ - L
_{β}:$$\mathrm{\Omega}(\mathit{f})={\displaystyle \sum _{j}{|{f}_{j}|}^{\beta}}$$ - Shannon entropy:$$\mathrm{\Omega}(\mathit{f})=-{\displaystyle \sum _{j}{f}_{j}\mathrm{ln}{f}_{j}}$$
_{j}< 1, - The Burg entropy:$$\mathrm{\Omega}(\mathit{f})={\displaystyle \sum _{j}\mathrm{ln}{f}_{j}}$$
_{j}> 0.

A second approach consists of considering f_{j} = E {U_{j}} or **f** = E {**U**} [41,41,42]. Again, here, U_{j} or **U** can, but need not, correspond to some physical quantities. In any case, we now want to assign a probability law
$\widehat{p}(\mathit{u})$ to it. Noting that the data **g** = **H f** = **HE** {**U**} = E {**HU**} can be considered as the constraints on it, we may need again a criterion to determine
$\widehat{p}(\mathit{u})$. Assuming then having some prior μ(**u**), we may maximize the relative entropy as that criterion. The mathematical problem then becomes:
The solution is:
where:
When
$\widehat{p}(\mathit{u})$ is obtained, we may be interested in computing:
which is the required solution.

$$\text{minimize}\phantom{\rule{0.2em}{0ex}}D[p(\mathit{u}):\mu (\mathit{u})]\phantom{\rule{0.2em}{0ex}}\text{subject}\phantom{\rule{0.2em}{0ex}}\mathrm{to}\phantom{\rule{0.2em}{0ex}}{\displaystyle \int \mathit{H}\phantom{\rule{0.2em}{0ex}}\mathit{u}\phantom{\rule{0.2em}{0ex}}p(\mathit{u})\phantom{\rule{0.2em}{0ex}}\mathrm{d}\mathit{u}=\mathit{g}}$$

$$\widehat{p}(\mathit{u})=\frac{1}{Z(\mathbf{\lambda})}\mu (\mathit{u})\mathrm{exp}[-{\mathbf{\lambda}}^{\prime}\mathit{H}\mathit{u}]$$

$$Z(\mathbf{\lambda})={\displaystyle \int \mu (\mathit{u})\mathrm{exp}[-{\mathbf{\lambda}}^{\prime}\mathit{H}\mathit{u}]\phantom{\rule{0.2em}{0ex}}\mathrm{d}\mathit{u}}.$$

$$\widehat{f}=\mathrm{E}\left\{U\right\}={\displaystyle \int \mathit{u}\widehat{p}(\mathit{u})\phantom{\rule{0.2em}{0ex}}\mathrm{d}\mathit{u}}$$

Interestingly, if we focus on
$\widehat{f}=\mathrm{E}\left\{U\right\}$, we will see that its expression depends on the choice of the prior μ(**u**). When μ(**u**) is separable:
$\mu (\mathit{u})={\displaystyle {\prod}_{j}{\mu}_{j}({u}_{j})}$, the expression of
$\widehat{p}(\mathit{u})$ will also be separable.

To go a little more into the details, let us introduce
$\mathit{s}=\mathit{H}{\mathbf{\lambda}}^{\prime}$ and define:
and its conjugate convex:
It can be shown easily that
$\widehat{f}=\mathrm{E}\left\{U\right\}$ can be obtained either via the dual
$\widehat{\mathbf{\lambda}}$ variables:
where
$\widehat{\mathbf{\lambda}}$ is obtained by:
or directly:
D(**λ**) is called the dual criterion and F (**f**) primal. However, it is not always easy to obtain an analytical expression for G(**s**) and its gradient G′(**s**). The functions F (**f**) and G(**s**) are conjugate convex.

$$G(\mathit{s})=\mathrm{ln}{\displaystyle \int \mu (\mathit{u})\mathrm{exp}[-{\mathit{s}}^{\prime}\mathit{u}]}\phantom{\rule{0.2em}{0ex}}\mathrm{d}\mathit{u}$$

$$F(\mathit{f})=\underset{\mathit{s}}{\mathrm{sup}}\{{\mathit{f}}^{\prime}\mathit{s}-G(\mathit{s})\}.$$

$$\widehat{\mathit{f}}={G}^{\prime}(\mathit{H}\widehat{\mathbf{\lambda}})$$

$$\widehat{\mathbf{\lambda}}=\mathrm{arg}\mathrm{min}\phantom{\rule{0.2em}{0ex}}\{D(\mathbf{\lambda})=\mathrm{ln}Z(\mathbf{\lambda})+{\mathbf{\lambda}}^{\prime}g\},$$

$$\widehat{\mathit{f}}=\underset{\{\mathit{f}:\mathit{Hf}=\mathit{g}\}}{\mathrm{arg}\mathrm{min}}\phantom{\rule{0.2em}{0ex}}\{F(\mathit{f})\}.$$

For the computational aspect, unfortunately, the cases where we may have analytical expressions for Z(**λ**) or G(**s**) = ln Z or F (**f**) are very limited. However, when there is analytical expressions for them, the computations can be done very easily. In Table 1, we summarizes some of those solutions:

In this section, we present in a brief way the Bayesian approach for the inverse problems in signal and image processing.

The different steps to find a solution to an inverse problem using the Bayesian approach can be summarized as follows:

- Assign a prior probability law p(
**ϵ**) to the modeling and observation errors, here**ϵ**. From this, find the expression of the likelihood p(**g**|**f**,**θ**_{1}). As an example, consider the Gaussian case:$$p(\in )=\text{N}(\in |0,{v}_{\in}\mathit{I})\to p(\mathit{g}|\mathit{f}=\mathcal{N}(\mathit{g}|\mathit{H}\phantom{\rule{0.2em}{0ex}}\mathit{f},{v}_{\in}\mathit{I}).$$_{1}in this case is the noise variance v_{ϵ}. - Assign a prior probability law p(
**f**|**θ**_{2}) to the unknown**f**to translate your prior knowledge on it. Again, as an example, consider the Gaussian case:$$p(\mathit{f})=\mathcal{N}(\mathit{f}|0,{v}_{f}\mathit{I})$$_{2}in this case is the variance v_{f}. - Apply the Bayes rule to obtain the expression of the posterior law:$$p(\mathit{f}|\mathit{g},{\mathit{\theta}}_{1},{\mathit{\theta}}_{2})=\frac{p(\mathit{g}|\mathit{f},{\mathit{\theta}}_{1})p(\mathit{f}|{\mathit{\theta}}_{2})}{p(\mathit{g}|{\mathit{\theta}}_{1},{\mathit{\theta}}_{2})}\propto p(\mathit{g}|\mathit{f},{\mathit{\theta}}_{1})p(\mathit{f}|{\mathit{\theta}}_{2}),$$
**g**|**f**,**θ**_{1}) is the likelihood, p(**f**|**θ**_{2}) the prior model,**θ**= [**θ**_{1},**θ**_{2}]′ their corresponding parameters (often called the hyper-parameters of the problem) and p(**g**|**θ**_{1},**θ**_{2}) is called the evidence of the model. - Use p(
**f**|**g**,**θ**_{1},**θ**_{2}) to infer any quantity dependent of**f**.

For the expressions of likelihood in (90) and the prior in (91), we obtain very easily the expression of the posterior:

$$p(\mathit{f}|\mathit{g},{v}_{\in},{v}_{f})=\mathcal{N}(\mathit{f}|\widehat{\mathit{f}},\widehat{\mathit{V}})\phantom{\rule{0.2em}{0ex}}\text{with}\phantom{\rule{0.2em}{0ex}}\widehat{\mathit{V}}={(\mathit{H}{\mathit{H}}^{\prime}+\frac{{v}_{\in}}{{v}_{f}}\mathit{I})}^{-1}and\phantom{\rule{0.2em}{0ex}}\widehat{\mathit{f}}=\widehat{\mathit{V}}{\mathit{H}}^{\prime}\mathit{g}$$

When the hyper-parameters **θ** can be fixed a priori, the problem is easy. In practice, we may use some summaries, such as:

- MAP:$${\widehat{\mathit{f}}}_{\mathit{MAP}}=\underset{\mathit{f}}{\mathrm{arg}\mathrm{max}}\{p(\mathit{f}|\mathit{g},\mathit{\theta})\}$$
- EAP or posterior mean (PM):$${\widehat{\mathit{f}}}_{EAP}={\displaystyle \int \mathit{f}p(\mathit{f}|\mathit{g},\mathit{\theta})}\mathrm{d}\mathit{f}$$

For the Gaussian case of (91), the MAP and EAP are the same and can be obtained by noting that:

$${\widehat{\mathit{f}}}_{\text{MAP}}=\underset{\mathit{f}}{{\displaystyle \mathrm{arg}\mathrm{min}}}\{J(\mathit{f})\}\phantom{\rule{0.2em}{0ex}}\text{with}\phantom{\rule{0.2em}{0ex}}J(\mathit{f})={\Vert \mathit{g}-\mathit{Hf}\Vert}_{2}^{2}+\mathrm{\lambda}{\Vert \mathit{f}\Vert}_{2}^{2},\phantom{\rule{0.2em}{0ex}}\text{where}\phantom{\rule{0.2em}{0ex}}\mathrm{\lambda}={v}_{\in}/{v}_{f}.$$

However, in real applications, the computation of even these simple point estimators may need efficient algorithm:

- For MAP, we need optimization algorithms, which can handle the huge dimensional criterion J(
**f**) = − ln p(**f**|**g**,**θ**). Very often, we may be limited to using gradient-based algorithms. - For EAP, we need integration algorithms, which can handle huge dimensional integrals. The most common tool here is the MCMC methods [24]. However, for real applications, very often, the computational costs are huge. Recently, different methods, called approximate Bayesian computation (ABC) [96–100] or VBA, have been proposed [74,96,98,101–107].

When the hyperparameters **θ** have also to be estimated, a prior p(**θ**) is assigned to them, and the expression of the joint posterior:
is obtained, which can then be used to infer them jointly. Very often, the expression of this joint posterior law is complex, and any computation may become very costly. The VBA methods try to approximate p(**f**, **θ**|**g**) by a simpler distribution, which can be handled more easily. Two particular and extreme cases are:
Obtaining the expressions of these approximated separable probability laws has to be done via a criterion. The natural criterion with some geometrical interpretation for the probability law manifolds is the Kullback–Leibler (KL) criterion:

$$p(\mathit{f},\mathit{\theta}|\mathit{g})=\frac{p(\mathit{g}|\mathit{f},{\mathit{\theta}}_{1})p(\mathit{f}|{\mathit{\theta}}_{2})p(\mathit{\theta})}{p(\mathit{g})}$$

- Bloc separable, such as q(
**f**,**θ**) = q_{1}(**f**) q_{2}(**θ**) or - Completely separable, such as $q(\mathit{f},\mathit{\theta})={\displaystyle {\prod}_{j}{q}_{1j}}({f}_{j}){\displaystyle {\prod}_{k}{q}_{2k}(}{\theta}_{k})$.

$$q(\mathit{f},\mathit{\theta})={q}_{1}(\mathit{f}){\displaystyle \prod _{k}{q}_{2k}({\theta}_{k})}$$

$$\mathrm{KL}[q:p]={\displaystyle \int q\mathrm{ln}\frac{q}{p}}={\langle \mathrm{ln}\frac{q}{p}\rangle}_{q}.$$

For hierarchical prior models with hidden variables z, the problem becomes more complex, because we have to give the expression of the joint posterior law:
and then approximate it by separable ones:
and then use them for estimation. See more discussions in [9,31,38,108–110]

$$p(\mathit{f},\mathit{z},\mathit{\theta}|\mathit{g})\propto p(\mathit{g}|\mathit{f},{\mathit{\theta}}_{1})p(\mathit{f}|\mathit{z},{\mathit{\theta}}_{2})p(\mathit{z}|{\mathit{\theta}}_{3})p(\mathit{\theta})$$

$$q(\mathit{f},z,\mathit{\theta}|\mathit{g})={q}_{1}(\mathit{f}){q}_{2}(z){q}_{3}(\mathit{\theta})\phantom{\rule{0.2em}{0ex}}\mathrm{or}\phantom{\rule{0.2em}{0ex}}q(\mathit{f},\mathit{\theta})={\displaystyle \prod _{j}{q}_{1j}({f}_{j}|{z}_{{f}_{j}})}{\displaystyle \prod _{j}{q}_{{2}_{j}}({z}_{{f}_{j}})}{\displaystyle \prod _{k}{q}_{3k}({\theta}_{k})}$$

In the following, first the general VBA method is detailed for the inference problems with hierarchical prior models. Then, a particular class of prior model (Student t) is considered, and the details of VBA algorithms for that are given.

To illustrate the basic ideas and tools, let us consider a vector **X** and its probability density function p(**x**), which we want to approximate by q(**x**) = ∏_{j} q_{j}(x_{j}). Using the KL criterion:
where we used the notation: 〈ln p(**x**)〉_{q} = ∫ q(**x**) ln p(**x**) d**x** and q_{−j}(**x**) = ∏_{i≠j} q_{i} (x_{i})

$$\begin{array}{l}\mathrm{KL}[q:p]={\displaystyle \int q(x)\mathrm{ln}\frac{q(x)}{p(x)}\phantom{\rule{0.2em}{0ex}}\mathrm{d}x}={\displaystyle \int q(x)\mathrm{ln}q(x)\phantom{\rule{0.2em}{0ex}}\mathrm{d}x-}{\displaystyle \int q(x)\mathrm{ln}p(x)\phantom{\rule{0.2em}{0ex}}\mathrm{d}x}\\ \phantom{\rule{3.9em}{0ex}}={\displaystyle \sum _{j}{\displaystyle \int {q}_{j}({x}_{j})\mathrm{ln}{q}_{j}({x}_{j})}\mathrm{d}{x}_{j}-{\langle \mathrm{ln}p(x)\rangle}_{q}}\\ \phantom{\rule{3.9em}{0ex}}={\displaystyle \sum _{j}{\displaystyle \int {q}_{j}({x}_{j})\mathrm{ln}{q}_{j}({x}_{j})\mathrm{d}{x}_{j}}-}{\displaystyle \int {q}_{j}({x}_{j})<\mathrm{ln}p(x){>}_{{q}_{-j}}}d{x}_{j}\end{array}$$

From here, trying to find the solution q_{i}, the basic method is an alternate optimization algorithm:
As we can see, the expression of q_{j}(x_{j}) depends on q (x_{i}), i ≠ j. It is not always possible to obtain analytical expressions for q_{j}(x_{j}). It is however possible to show that, if p(**x**) is a member of exponential families, then q_{j}(x_{j}) are also members of exponential families. These iterations then become much simpler, because at each iteration, we need to update the parameters of the exponential families. To go a little more into the details, let us consider some particular simple cases.

$${q}_{j}({x}_{j})\propto \mathrm{exp}[<\mathrm{ln}p(\mathit{x}){>}_{{q}_{-j}}].$$

In the case of two variables **x** = [x_{1}, x_{2}]′, we have:
As an illustrative example, consider the case where we want to approximate p(x_{1}, x_{2}) by q(x_{1}, x_{2}) = q_{1}(x_{1}) q_{2}(x_{2}) to be able to compute the expected values:
which need double integrations when p(x_{1}, x_{2}) is not separable in its two variables. If we can do that separable approximation, then, we can compute:
which needs only 1D integrals. Let us see if
$({\tilde{\mu}}_{1},{\tilde{\mu}}_{2})$ will converge to (m_{1}, m_{2}). To illustrate this, let us consider the very simple case of the Gaussian:
It is then easy to see that
${q}_{1}({x}_{1})=\mathcal{N}({x}_{1}|{\tilde{\mu}}_{1},{\tilde{v}}_{1})$ and
${q}_{2}({x}_{2})=\mathcal{N}({x}_{2}|{\tilde{\mu}}_{2},{\tilde{v}}_{2})$ and that:
with:
See [111] for details and where we showed that, initializing the algorithm with
${\tilde{\mu}}_{1}^{(0)}=0$ and
${\tilde{\mu}}_{2}^{(0)}=0$, the means converges to the right values m_{1} and m_{2}, However, we may be careful about the convergence of the variances.

$$\{\begin{array}{c}{q}_{1}\left({x}_{1}\right)\propto \mathrm{exp}\left[<\mathrm{ln}p(\mathit{x})>{q}_{2}({x}_{2})\right]\\ {q}_{2}\left({x}_{2}\right)\propto \mathrm{exp}\left[<\mathrm{ln}p(\mathit{x})>{q}_{1}({x}_{1})\right]\end{array}$$

$$\{\begin{array}{c}{m}_{1}=E\left\{{x}_{1}\right\}={\displaystyle \iint {x}_{1}}p\left({x}_{1},{x}_{2}\right)\mathrm{d}{x}_{1}\phantom{\rule{0.2em}{0ex}}\mathrm{d}{x}_{2}\\ {m}_{2}=\left\{{x}_{2}\right\}={\displaystyle \iint {x}_{2}p\left({x}_{1},{x}_{2}\right)}\mathrm{d}{x}_{1}\phantom{\rule{0.2em}{0ex}}\mathrm{d}{x}_{2}\end{array}$$

$$\{\begin{array}{c}{\tilde{\mu}}_{1}=\mathrm{E}\left\{{x}_{1}\right\}={\displaystyle \int {x}_{1}{q}_{1}\left({x}_{1}\right)}\mathrm{d}{x}_{1}\\ {\tilde{\mu}}_{2}=\mathrm{E}\left\{{x}_{2}\right\}={\displaystyle \int {x}_{2}{q}_{2}\left({x}_{2}\right)}d{x}_{2}\end{array}$$

$$p({x}_{1},{x}_{2})=\mathcal{N}\left(\left[\begin{array}{c}{x}_{1}\\ {x}_{2}\end{array}\right]|\left[\begin{array}{c}{m}_{1}\\ {m}_{2}\end{array}\right],\left[\begin{array}{c}{v}_{1}\\ \rho \sqrt{{v}_{1}{v}_{2}}\end{array}\begin{array}{c}\rho \sqrt{{v}_{1}{v}_{2}}\\ {v}_{2}\end{array}\right]\right).$$

$$\{\begin{array}{c}{q}_{1}^{(k+1)}({x}_{1})=p({x}_{1})|{x}_{2}={\tilde{\mu}}_{2}^{(k)}=\mathcal{N}\left({x}_{1}|{\tilde{\mu}}_{1}^{(k)},{\tilde{v}}_{1}^{(k)}\right)\\ {q}_{1}^{(k+1)}({x}_{2})=p({x}_{2}|{x}_{1}={\tilde{\mu}}_{1}^{(k)}=\mathcal{N}\left({x}_{1}|{\tilde{\mu}}_{2}^{(k)},{\tilde{v}}_{2}^{(k)}\right)\end{array}$$

$$\{\begin{array}{lll}{\tilde{\mu}}_{1}^{(k+1)}\hfill & =\hfill & {m}_{1}+\rho \sqrt{{v}_{1}/{v}_{2}}\left({\tilde{\mu}}_{2}^{(k)}-{m}_{2}\right)\hfill \\ {\tilde{v}}_{1}^{(k+1)}\hfill & =\hfill & \left(1-{\rho}^{2}\right){v}_{1}\hfill \\ {\tilde{\mu}}_{2}^{(k+1)}\hfill & =\hfill & {m}_{2}+\rho \sqrt{{v}_{2}/{v}_{1}}\left({\tilde{\mu}}_{1}^{(k)}-{m}_{1}\right)\hfill \\ {\tilde{v}}_{2}^{(k+1)}\hfill & =\hfill & (1-{\rho}^{2}){v}_{2}\hfill \end{array}.$$

As we could see, to be able to use such an algorithm in practical cases, we need to be able to compute < ln p(**x**) >_{q}_{2(}_{x}_{2)} and < ln p(**x**) >_{q}_{1(}_{x}_{1)}. Only for a few cases can we can do this analytically. Different algorithms can be obtained depending on the choice of a particular family for q_{j}(x_{j}) [103,112–120].

To show this, let us consider the exponential family:
where **θ** is a vector of parameter and g(**θ**) and **u**(**x**) are known functions.

$$p(\mathit{x}|\mathit{\theta})=g(\mathit{\theta})\mathrm{exp}[{\mathit{\theta}}^{\prime}\mathit{u}(\mathit{x})]$$

This parametric exponential family has the following conjugacy property: For a given prior p(**θ**) in the family:
the corresponding posterior:
is in the same family.

$$p(\mathit{\theta}|\eta ,\mathit{\nu})=h(\eta ,\mathit{\nu})g{(\mathit{\theta})}^{\eta}\mathrm{exp}[{\mathit{\nu}}^{\prime}\mathit{\theta}]$$

$$\begin{array}{c}p(\mathit{\theta}|\mathit{x})\propto p(\mathit{x}|\mathit{\theta})p(\mathit{\theta}|\eta ,\mathit{\nu})\\ \propto g{(\mathit{\theta})}^{\eta +1}\mathrm{exp}[\mathit{\nu}+\mathit{u}(\mathit{x}){]}^{\prime}\mathit{\theta}]\\ \propto p(\mathit{\theta}|\eta +1,\mathit{\nu}+\mathit{u}(\mathit{x}))\end{array}$$

For this family, we have:
It is then easy to show that:
which are in the same exponential family. This simplifies greatly the computations, thanks to the fact that, in each iteration, we only need to compute
$\tilde{\mathit{u}}(\mathit{x})=\langle \mathit{u}{\left(\mathit{x}\right)}_{{q}_{-j}}\rangle $ and update the parameters.

$${\langle \mathrm{ln}p\left(\mathit{x}|\mathit{\theta}\right)\rangle}_{q}=\mathrm{ln}g\left(\mathit{\theta}\right)+{\mathit{\theta}}^{\prime}{\langle \mathit{u}\left(\mathit{x}\right)\rangle}_{q}.$$

$${q}_{j}({x}_{j})\propto g(\mathit{\theta})\mathrm{exp}\left[{\mathit{\theta}}^{\prime}{\langle \mathit{u}\left(\mathit{x}\right)\rangle}_{{q}_{-j}}\right]$$

Now, if we consider:
$$p(\mathit{x}|\mathit{\theta})=g(\mathit{\theta})\mathrm{exp}[{\mathit{\theta}}^{\prime}\mathit{u}(\mathit{x})]$$
with a prior on **θ**:
$$p(\mathit{\theta}|\eta ,\mathit{\nu})=h(\eta ,\mathit{\nu})g{(\mathit{\theta})}^{\eta}\mathrm{exp}[{\mathit{\nu}}^{\prime}\mathit{\theta}]$$
and the joint p(x, **θ**|η, **ν**) = p(**x**|**θ**) p(**θ**|η, **ν**), which is not separable in **x** and **θ**, and we want to approximate it with the separable q(**x**, **θ**) = q_{1}(**x**) q_{2}(**θ**), then we will have:
where
$\tilde{\mathit{\mu}}={\langle \mathit{u}\left(\mathit{x}\right)\rangle}_{{q}_{1}\left(\mathit{x}\right)}$.

$$\{\begin{array}{c}q(\mathit{\theta})=h(\tilde{\eta},\tilde{\mathit{\nu}})g{(\mathit{\theta})}^{\tilde{\eta}}\mathrm{exp}[{\tilde{\mathit{\nu}}}^{\prime}\mathit{\theta}]\\ q(\mathit{x})=g(\tilde{\mathit{\theta}})\mathrm{exp}[{\tilde{\mathit{\theta}}}^{\prime}\mathit{u}(x)]\end{array}\phantom{\rule{0.2em}{0ex}}\text{with}\phantom{\rule{0.2em}{0ex}}\{\begin{array}{lll}\tilde{\eta}\hfill & =\hfill & \eta +1\hfill \\ \tilde{\mathit{\nu}}\hfill & =\hfill & \mathit{\nu}+\tilde{\mathit{u}}(\mathit{x})\hfill \\ \tilde{\mathit{\theta}}\hfill & =\hfill & \tilde{\mathit{\nu}}\hfill \end{array}$$

Before going into the details and for similarity with the notations in the next sections, we replace x by f, such that now we are trying to approximate p(**f**, **θ**) = p(**f**|**θ**) p(**θ**) by a separable q(**f**, **θ**) = q_{1}(**f**) q_{2}(**θ**). Interestingly, depending on the choice of the family laws for q_{1} and q_{2}, we obtain different algorithms:

- ${q}_{1}(\mathit{f})=\delta (\mathit{f}-\tilde{\mathit{f}})$ and ${q}_{2}(\mathit{\theta})=\delta (\mathit{\theta}-\tilde{\mathit{\theta}})$. In this case, we have:$$\{\begin{array}{c}{q}_{1}(\mathit{f})\propto \mathrm{exp}[<\mathrm{ln}p(\mathit{f},\mathit{\theta}){>}_{{q}_{2}}]\propto \mathrm{exp}\left[\mathrm{ln}p(\mathit{f},\tilde{\mathit{\theta}})\right]\propto p(\mathit{f},\mathit{\theta}=\tilde{\mathit{\theta}})\propto p(\mathit{f}|\mathit{\theta}=\tilde{\mathit{\theta}})\\ {q}_{2}(\mathit{\theta})\propto \mathrm{exp}[<\mathrm{ln}p(\mathit{f},\mathit{\theta}){>}_{{q}_{1}}]\propto \mathrm{exp}\left[\mathrm{ln}p(\tilde{\mathit{f}},\mathit{\theta})\right]\propto p(\mathit{f}=\tilde{\mathit{f}},\mathit{\theta})\propto p(\mathit{\theta}|\mathit{f}=\tilde{\mathit{f}})\end{array}$$$$\{\begin{array}{c}\tilde{\mathit{f}}=\mathrm{arg}{\mathrm{max}}_{\mathit{f}}\left\{p(\mathit{f},\mathit{\theta})=\tilde{\mathit{\theta}}\right\}\\ \tilde{\mathit{\theta}}=\mathrm{arg}{\mathrm{max}}_{\mathit{\theta}}\left\{p(\mathit{f}=\tilde{\mathit{f}},\mathit{\theta})\right\}\end{array}$$$$(\tilde{\mathit{f}},\tilde{\mathit{\theta}})=\underset{(\mathit{f},\mathit{\theta})}{\mathrm{arg}\mathrm{max}}\{p(\mathit{f},\mathit{\theta})\}.$$
**f**are not used for the estimation of**θ**and the uncertainties of**θ**are not used for the estimation of**f**. - q
_{1}(**f**) is free form and ${q}_{2}(\mathit{\theta})=\delta (\mathit{\theta}-\tilde{\mathit{\theta}})$ In the same way, this time we obtain:$$\{\begin{array}{l}<\mathrm{ln}p(\mathit{f},\mathit{\theta}){>}_{{q}_{2}(\mathit{\theta})}=\mathrm{ln}p(\mathit{f},\tilde{\mathit{\theta}})\hfill \\ <\mathrm{ln}p(\mathit{f},\mathit{\theta}){>}_{{q}_{1}(\mathit{f})}=<\mathrm{ln}p(\mathit{f},\mathit{\theta}){>}_{{q}_{1}(\mathit{f}|\tilde{\mathit{\theta}})}=Q(\mathit{\theta},\tilde{\mathit{\theta}})\hfill \end{array}$$$$\{\begin{array}{l}\begin{array}{l}{q}_{1}(\mathit{f})\propto \mathrm{exp}\left[\mathrm{ln}p(\mathit{f},\mathit{\theta}=\tilde{\mathit{\theta}})\right]\propto p(\mathit{f},\tilde{\mathit{\theta}})\hfill \\ {q}_{2}(\mathit{\theta})\propto \mathrm{exp}\left[Q(\mathit{\theta},\tilde{\mathit{\theta}})\right]\to \tilde{\mathit{\theta}}=\mathrm{arg}{\mathrm{max}}_{\mathit{\theta}}\left\{Q(\mathit{\theta},\tilde{\mathit{\theta}})\right\}\hfill \end{array}\hfill \end{array}$$**f**are used for the estimation of**θ**, but the uncertainties of**θ**are not used for the estimation of**f**. - ${q}_{1}(\mathit{f})=\delta (\mathit{f}-\tilde{\mathit{f}})$ and q
_{2}(**θ**) is free form. In the same way, this time we obtain:$$\{\begin{array}{l}<\mathrm{ln}p(\mathit{f},\mathit{\theta}){>}_{{q}_{1}(\mathit{f})}=\mathrm{ln}p(\mathit{f}=\tilde{\mathit{f}},\mathit{\theta})\hfill \\ <\mathrm{ln}p(\mathit{f},\mathit{\theta}){>}_{{q}_{2}(\mathit{\theta})}=<\mathrm{ln}p(\mathit{f},\mathit{\theta}){>}_{p(\mathit{\theta}|\mathit{f}=\tilde{\mathit{f}})}=Q(\tilde{\mathit{f}},\mathit{f})\hfill \end{array}$$$$\{\begin{array}{l}{q}_{2}(\theta )\propto \mathrm{ln}p(\mathit{f}=\tilde{\mathit{f}},\mathit{\theta})=p(\mathit{\theta}|\mathit{f}=\tilde{\mathit{f}})\hfill \\ {q}_{1}(\mathit{f})\propto \mathrm{exp}\left[Q(\tilde{\mathit{f}},\mathit{\theta})\right]\to \tilde{\mathit{\theta}}=\mathrm{arg}{\mathrm{max}}_{\mathit{\theta}}\left\{Q(\mathit{f}=\tilde{\mathit{f}},\mathit{\theta})\right\}\hfill \end{array}$$**f**are used for the estimation of**θ**, but the uncertainties of**θ**are not used for the estimation of**f**. - Both q
_{1}(f) and q_{2}(θ) have free form. The main difficulty here is that, at each iteration, the expression of q_{1}and q_{2}may change. However, if p(**f**,**θ**) is in the generalized exponential family, the expressions of q_{1}(**f**) and q_{2}(**θ**) will also be in the same family, and we have only to update the parameters at each iteration.

As a simple example, consider the Gaussian case where
$p\left(\mathit{g}|\mathit{f},{\theta}_{1}\right)=\mathcal{N}\left(\mathit{g}|\mathit{Hf},\left(1/{\theta}_{1}\right)\mathit{I}\right),p\left(\mathit{f}|{\theta}_{2}\right)=\mathcal{N}\left(\mathit{f}|0,\left(1/{\theta}_{2}\right)\mathit{I}\right)$ and
$p\left({\theta}_{1}\right)=\mathcal{G}({\theta}_{1}|{\alpha}_{10},{\beta}_{10})p({\theta}_{2})=\mathcal{G}({\theta}_{2}|{\alpha}_{20},{\beta}_{20})$, and so, we have:
From this expression J(**f**, θ_{1}, θ_{2}) = ln p(**f**, θ_{1}, θ_{2}|**g**), it is easy to obtain the equations of an alternate JMAP algorithm by computing the derivatives of it with respective to its arguments and equating them to zero:

$$\begin{array}{c}\mathrm{ln}p(\mathit{f},{\theta}_{1,}{\theta}_{2}|\mathit{g})=\frac{\mathrm{M}}{2}\mathrm{ln}{\theta}_{1}-\frac{{\theta}_{1}}{2}\Vert \mathit{g}-\mathit{Hf}\Vert {}_{2}^{2}+\frac{\mathcal{N}}{2}\mathrm{ln}{\theta}_{2}-\frac{{\theta}_{2}}{2}\Vert \mathit{f}\Vert {}_{2}^{2}\\ +({\alpha}_{10}-1)\mathrm{ln}{\theta}_{1}-{\beta}_{10}{\theta}_{1}+({\alpha}_{20}-1)\mathrm{ln}{\theta}_{2}-{\beta}_{20}{\theta}_{2}.\end{array}$$

$$\begin{array}{l}\frac{\partial J}{\partial \mathit{f}}=0\to \mathit{f}={({\mathit{H}}^{\prime}\mathit{H}+\mathrm{\lambda}\mathit{I})}^{-1}{\mathit{H}}^{\prime}\mathit{g}\phantom{\rule{0.2em}{0ex}}\mathrm{with}\phantom{\rule{0.2em}{0ex}}\mathrm{\lambda}=\frac{{\theta}_{2}}{{\theta}_{1}}\hfill \\ \frac{\partial J}{\partial {\theta}_{1}}=0\to {\theta}_{1}=\frac{{\tilde{\alpha}}_{1}}{{\tilde{\beta}}_{1}}\phantom{\rule{0.2em}{0ex}}\text{with}\phantom{\rule{0.2em}{0ex}}{\tilde{\alpha}}_{1}=({\alpha}_{10}-1)+\frac{\mathrm{M}}{2}\phantom{\rule{0.2em}{0ex}}\text{and}\phantom{\rule{0.2em}{0ex}}{\tilde{\beta}}_{1}={\beta}_{10}+\frac{1}{2}\Vert \mathit{g}-\mathit{Hf}\Vert {}_{2}^{2}\hfill \\ \frac{\partial J}{\partial {\theta}_{2}}=0\to {\theta}_{1}=\frac{{\tilde{\alpha}}_{2}}{{\tilde{\beta}}_{2}}\phantom{\rule{0.2em}{0ex}}\text{with}\phantom{\rule{0.2em}{0ex}}{\tilde{\alpha}}_{2}=({\tilde{\alpha}}_{20}-1)+\frac{\mathrm{M}}{2}\phantom{\rule{0.2em}{0ex}}\text{and}\phantom{\rule{0.2em}{0ex}}{\tilde{\beta}}_{2}={\beta}_{20}+\frac{1}{2}\Vert \mathit{f}\Vert {}_{2}^{2}\hfill \end{array}$$

From the expression of the joint probability law p(**f**, θ_{1}, θ_{2}|**g**), we can also obtain the expressions of the conditionals:
However, obtaining analytical expressions of the marginals p(**f**|**g**), p(θ_{1}|**g**) and p(θ_{2}|**g**) is not easy. We can then obtain approximate expressions q_{1}(**f**|**g**), q_{2}(θ_{1}|**g**) and q_{3}(θ_{2}|**g**) using the VBA method. For this case, thanks to the conjugacy property, we have:
We can then compare the three algorithms in Table 2:

$$\{\begin{array}{l}\begin{array}{l}p(\mathit{f}|\mathit{g},{\theta}_{1},{\theta}_{2})=\mathcal{N}(\mathit{f}|\tilde{\mathit{f}},\tilde{\mathit{V}})\\ \phantom{\rule{0.2em}{0ex}}\text{with}\phantom{\rule{0.2em}{0ex}}\tilde{\mathit{V}}={({\mathit{H}}^{\prime}\mathit{H}+\mathrm{\lambda}\mathit{I})}^{-1},\phantom{\rule{1.2em}{0ex}}\tilde{\mathit{f}}=\tilde{\mathit{V}}{\mathit{H}}^{\prime}\mathit{g},\phantom{\rule{0.2em}{0ex}}\mathrm{\lambda}=\frac{{\theta}_{2}}{{\theta}_{1}}\end{array}\hfill \\ \begin{array}{l}p({\theta}_{1}|\mathit{g},\mathit{f},{\theta}_{2})=\mathcal{G}({\theta}_{1}|{\tilde{\alpha}}_{1},{\tilde{\beta}}_{1})\\ \phantom{\rule{0.2em}{0ex}}\text{with}\phantom{\rule{0.2em}{0ex}}{\tilde{\alpha}}_{1}=({\alpha}_{10}-1)+\frac{\mathrm{M}}{2},\phantom{\rule{1.2em}{0ex}}{\tilde{\beta}}_{1}={\beta}_{10}+\frac{1}{2}\Vert \mathit{g}-\mathit{Hf}\Vert {}_{2}^{2}\end{array}\hfill \\ \begin{array}{l}p({\theta}_{2}|\mathit{g},\mathit{f},{\theta}_{1})=\mathcal{G}({\theta}_{2}|{\tilde{\alpha}}_{2},{\tilde{\beta}}_{2})\\ \phantom{\rule{0.2em}{0ex}}\text{with}\phantom{\rule{0.2em}{0ex}}{\tilde{\alpha}}_{2}=({\alpha}_{20}-1)+\frac{M}{2},\phantom{\rule{1.2em}{0ex}}{\tilde{\beta}}_{2}={\beta}_{20}+\frac{1}{2}\Vert \mathit{f}\Vert {}_{2}^{2}\end{array}\hfill \end{array}$$

$$\{\begin{array}{l}\begin{array}{l}q(f)=\mathcal{N}(\mathit{f}|\tilde{\mathit{f}},\tilde{\mathit{V}})\phantom{\rule{0.2em}{0ex}}\\ \text{with}\phantom{\rule{0.2em}{0ex}}\tilde{\mathit{V}}={({\mathit{H}}^{\prime}\mathit{H}+{\lambda}^{\prime}\mathit{I})}^{-1},\phantom{\rule{1.2em}{0ex}}\tilde{\mathit{f}}=\tilde{\mathit{V}}{\mathit{H}}^{\prime}\mathit{g},\phantom{\rule{0.2em}{0ex}}\tilde{\lambda}=\frac{<{\theta}_{2}>}{<{\theta}_{1}>};\end{array}\hfill \\ \begin{array}{l}q({\theta}_{1}=\mathcal{G}({\theta}_{1}|{\tilde{\alpha}}_{1},{\tilde{\beta}}_{1})\phantom{\rule{0.2em}{0ex}}\\ \text{with}\phantom{\rule{0.2em}{0ex}}{\tilde{\alpha}}_{1}=({\alpha}_{10}-1)+\frac{\mathrm{M}}{2},\phantom{\rule{1.2em}{0ex}}{\tilde{\beta}}_{1}={\beta}_{10}+\frac{1}{2}<\Vert \mathit{g}-\mathit{Hg}\Vert {}_{2}^{2}>\end{array}\hfill \\ \begin{array}{l}p({\theta}_{2}|\mathit{g},\mathit{f})=\mathcal{G}({\theta}_{2}|{\tilde{\alpha}}_{2},{\tilde{\beta}}_{2})\phantom{\rule{0.2em}{0ex}}\\ \text{with}\phantom{\rule{0.2em}{0ex}}{\tilde{\alpha}}_{2}=({\alpha}_{20}-1)+\frac{\mathcal{N}}{2},\phantom{\rule{1.2em}{0ex}}{\tilde{\beta}}_{2}={\beta}_{20}+\frac{1}{2}<\Vert \mathit{f}\Vert {}_{2}^{2}>\end{array}\hfill \end{array}$$

It is important to remark that, in JMAP, the computation of **f** can be done via the optimization of the criterion J(**f**, θ_{1}, θ_{2}) = ln p(**f**, θ_{1}, θ_{2}|**g**), which does not need explicitly the matrix inversion of
$\tilde{\mathit{V}}={({\mathit{H}}^{\prime}\mathit{H}+{\lambda}^{\prime}\mathit{I})}^{-1}$ However, in BEM and VBA, we need to compute it due to the following requirements:

$$\begin{array}{l}\phantom{\rule{3.1em}{0ex}}<f{>}_{q}=\tilde{f},\\ \phantom{\rule{2.1em}{0ex}}<\Vert \mathit{f}\Vert {}^{2}{>}_{q}=\mathrm{tr}\left(<\tilde{\mathit{f}}{\tilde{\mathit{f}}}^{\prime}{>}_{q}\right)=\mathrm{tr}\left(\tilde{\mathit{f}}{\tilde{\mathit{f}}}^{\prime}+\tilde{\mathit{V}}\right)+\Vert \tilde{\mathit{f}}\Vert {}^{2}+\mathrm{tr}\left(\tilde{\mathit{V}}\right),\\ \phantom{\rule{2.7em}{0ex}}<{f}_{j}^{2}{>}_{q}={\left[\tilde{\mathit{V}}\right]}_{jj}+{\tilde{f}}_{j}^{2},\\ <\Vert \mathit{g}-\mathit{Hf}\Vert {}^{2}{>}_{q}=\left[{\mathit{g}}^{\prime}\mathit{g}-2{\langle \mathit{f}\rangle}^{\prime}{}_{q}{\mathit{H}}^{\prime}\mathit{g}+{\mathit{H}}^{\prime}{\langle \tilde{\mathit{f}}\mathit{f}\rangle}_{q}\mathit{H}\right]\\ \phantom{\rule{5.8em}{0ex}}=\left[{\mathit{g}}^{\prime}\mathit{g}-2{\tilde{\mathit{f}}}^{\prime}\tilde{\mathit{H}}\mathit{g}+{\mathit{H}}^{\prime}\left(\tilde{\mathit{V}}+\tilde{\mathit{f}}{\tilde{\mathit{f}}}^{\prime}\right)\mathit{H}\right]\\ \phantom{\rule{5.8em}{0ex}}=\Vert \mathit{g}-\mathit{H}\tilde{\mathit{f}}\Vert {}^{2}+\mathrm{tr}\left({\mathit{H}}^{\prime}\tilde{\mathit{V}}\mathit{H}\right)\end{array}$$

For some extensions and more details, see [111].

For a linear inverse problem:
with an assigned likelihood p(**g**|**f**, **θ**_{1}), when a hierarchical prior model p(**f**|**z**, **θ**_{2}) p(**z**|**θ**_{3}) is used and when the estimation of the hyper-parameters **θ** = [**θ**_{1}, **θ**_{2}, **θ**_{3}]′ has to be considered, the joint posterior law of all the unknowns becomes:
The main idea behind the VBA is to approximate this joint posterior by a separable one, for example: q(**f**, **z**, **θ**|**g**) = q_{1}(**f**) q_{2}(**z**) q_{3}(**θ**) and where the expressions of q(**f**, **z**, **θ**|**g**) are obtained by minimizing the Kullback–Leibler divergence (99), as explained in previous section. This approach can also be used for model selection based on the evidence of the model ln p(**g**) [121] where:
Interestingly, it is easy to show that:
where
$\mathcal{F}\left(q\right)$ is the free energy associated with q defined as:
Therefore, for a given model
$\mathcal{M}$, minimizing KL [q : p] is equivalent to maximizing
$\mathcal{F}\left(q\right)$ and when optimized,
$\mathcal{F}\left(q*\right)$ gives a lower bound for ln p(**g**). Indeed, the name variational approximation is due to the fact that
$\mathrm{ln}p(g)\ge \mathcal{F}(q)$, and so,
$\mathcal{F}\left(q\right)$ is a lower bound to the evidence ln p(**g**).

$$\mathrm{M}:g=\mathit{Hf}+\in $$

$$p(\mathit{f},\mathit{z},\mathit{\theta}|\mathit{g})=\frac{p(\mathit{f},\mathit{z},\mathit{\theta}|\mathit{g})}{p(\mathit{g})}=\frac{p(\mathit{g},{\mathit{\theta}}_{1})p(\mathit{f}|\mathit{z},{\mathit{\theta}}_{2})p(\mathit{z}|{\mathit{\theta}}_{3})p(\mathit{\theta})}{p(\mathit{g})}$$

$$p(\mathit{g})={\displaystyle \iiint p(\mathit{f},\mathit{z},\mathit{\theta},\mathit{g})}\phantom{\rule{0.2em}{0ex}}\mathrm{d}f\phantom{\rule{0.2em}{0ex}}\mathrm{d}z\phantom{\rule{0.2em}{0ex}}\mathrm{d}\theta .$$

$$\mathrm{ln}p(\mathit{g})=\mathrm{KL}\left[q:p\right]+\mathcal{F}(q)$$

$$\mathcal{F}(q)={\langle \mathrm{ln}\frac{p(\mathit{f},\mathit{z},\mathit{\theta},\mathit{g})}{q(\mathit{f},\mathit{z},\mathit{\theta})}\rangle}_{q}$$

Without any other constraint than the normalization of q, an alternate optimization of
$\mathcal{F}\left(q\right)$ with respect to q_{1}, q_{2} and q_{3} results in:
Note that these relations represent an implicit solution for q_{1}(**f**), q_{2}(z) and q_{3}(**θ**), which need, at each iteration, the expression of the expectations in the right hand of exponentials. If p(**g**|**f**, **z**, **θ**_{1}) is a member of an exponential family and if all of the priors p(**f**|**z**, **θ**_{2}), p(**z**|**θ**_{3}), p(**θ**_{1}), p(**θ**_{2}) and p(**θ**_{3}) are conjugate priors, then it is easy to see that these expressions lead to standard distributions for which the required expectations are easily evaluated. In that case, we may note:
where the tilded quantities
$\tilde{\mathit{z}}$,
$\tilde{\mathit{f}}$ and
$\tilde{\mathit{\theta}}$are, respectively, functions of (
$\tilde{\mathit{f}}$,
$\tilde{\mathit{\theta}}$), (
$\tilde{\mathit{f}}$,
$\tilde{\mathit{z}}$) and (
$\tilde{\mathit{f}}$,
$\tilde{\mathit{z}}$). This means that the expression of
${q}_{1}(\mathit{f}|\tilde{\mathit{z}},\tilde{\mathit{\theta}})$ depends on (
$\tilde{\mathit{f}}$,
$\tilde{\mathit{z}}$), the expression of
${q}_{2}(\mathit{z}|\tilde{\mathit{f}},\tilde{\mathit{\theta}})$ depends on (
$\tilde{\mathit{z}}$,
$\tilde{\mathit{\theta}}$) and the expression of
${q}_{3}(\mathit{\theta}|\tilde{\mathit{f}},\tilde{\mathit{z}})$ depends on (
$\tilde{\mathit{f}}$,
$\tilde{\mathit{z}}$). With this notation, the alternate optimization results in alternate updating of the parameters (
$\tilde{\mathit{z}}$,
$\tilde{\mathit{\theta}}$) of q_{1}, the parameters (
$\tilde{\mathit{f}}$,
$\tilde{\mathit{\theta}}$) of q_{2} and the parameters (
$\tilde{\mathit{f}}$,
$\tilde{\mathit{z}}$) of q_{3}. Finally, we may note that, to monitor the convergence of the algorithm, we may evaluate the free energy:
Other decompositions for q(**f**, **z, θ**) are also possible. For example: q(**f**, **z, θ**) = q_{1}(**f**|**z**) q_{2}(**z**) q_{3}(**θ**) or even:
$q(\mathit{f},\mathit{z},\mathit{\theta})={\displaystyle {\prod}_{j}{q}_{1j}({f}_{j})}\phantom{\rule{0.2em}{0ex}}{\displaystyle {\prod}_{j}{q}_{2j}({z}_{fj})\phantom{\rule{0.2em}{0ex}}{\displaystyle {\prod}_{l}{q}_{3l}({\theta}_{l})}}$. Here, we consider the first case and give some more details on it.

$$\{\begin{array}{l}{q}_{1}(\mathit{f})\propto \mathrm{exp}\left[-{\langle \mathrm{ln}p\left(\mathit{f},\mathit{z},\mathit{\theta},\mathit{g}\right)\rangle}_{q{(\mathit{z})}_{q}(\mathit{\theta})}\right],\hfill \\ {q}_{2}(\mathit{z})\propto \mathrm{exp}\left[-{\langle \mathrm{ln}p\left(\mathit{f},\mathit{z},\mathit{\theta},\mathit{g}\right)\rangle}_{q{(\mathit{f})}_{q}(\mathit{\theta})}\right],\hfill \\ {q}_{2}(\mathit{\theta})\propto \mathrm{exp}\left[-{\langle \mathrm{ln}p\left(\mathit{f},\mathit{z},\mathit{\theta},\mathit{g}\right)\rangle}_{q{(\mathit{f})}_{q}(\mathit{z})}\right],\hfill \end{array}$$

$$q(\mathit{f},\mathit{z},\mathit{\theta})={q}_{1}(\mathit{f}|\tilde{\mathit{z}},\tilde{\mathit{\theta}}){q}_{2}(\mathit{z}|\tilde{\mathit{f}},\tilde{\mathit{\theta}}){q}_{3}(\mathit{\theta}|\tilde{\mathit{f}},\tilde{\mathit{z}})$$

$$\begin{array}{l}\mathcal{F}(q)={\langle \mathrm{ln}p(\mathit{f},\mathit{z},\mathit{\theta},\mathit{g}\phantom{\rule{0.2em}{0ex}})\rangle}_{q}-{\langle \mathrm{ln}q(\mathit{f},\mathit{z},\mathit{\theta})\rangle}_{q}\\ \phantom{\rule{2.6em}{0ex}}={\langle \mathrm{ln}p(\mathit{g}|\mathit{f},\mathit{z},\mathit{\theta})\rangle}_{q}+{\langle \mathrm{ln}p(\mathit{f}|\mathit{z},\mathit{\theta})\rangle}_{q}+{\langle \mathrm{ln}p(\mathit{z}|\mathit{\theta})\rangle}_{q}+{\langle \mathrm{ln}p(\mathit{\theta})\rangle}_{q}-{\langle \mathrm{ln}q(\mathit{f})\rangle}_{q}-{\langle \mathrm{ln}q(\mathit{z})\rangle}_{q}-{\langle \mathrm{ln}q(\mathit{\theta})\rangle}_{q}.\end{array}$$

The Student t model is:
The Cauchy model is obtained when ν = 1. Knowing that:
we can write this model via the positive hidden variables z_{fj}:

$$p(\mathit{f}|\nu )={\displaystyle \prod _{j}\mathcal{S}t}({f}_{j}|\nu )\phantom{\rule{0.2em}{0ex}}\text{with}\phantom{\rule{0.2em}{0ex}}\mathcal{S}t({f}_{j}|\nu )=\frac{1}{\sqrt{\pi \nu}}\frac{\Gamma ((\nu +1)/2)}{\Gamma (\nu /2)}{\left(1+{f}_{j}^{2}/\nu \right)}^{-(\nu +1)/2}$$

$$\mathcal{S}t({f}_{j}|\nu )={\displaystyle {\int}_{0}^{\infty}\mathcal{N}}({f}_{j}|0,1/{z}_{{f}_{j}})\mathcal{G}({z}_{{f}_{j}}|\nu /2,\nu /2)\mathrm{d}{z}_{{f}_{j}}$$

$$\{\begin{array}{l}p({f}_{j}|{z}_{{j}_{f}})=\mathcal{N}({f}_{j}|0,1/{z}_{{f}_{j}})\propto \mathrm{exp}\left[-\frac{1}{2}{z}_{{f}_{j}}{f}_{j}^{2}\right]\hfill \\ p({z}_{{f}_{j}}|\alpha ,\beta )=\mathcal{G}({z}_{{j}_{j}}|\alpha ,\beta )\propto {z}_{{f}_{j}}^{(\alpha -1)}\mathrm{exp}\left[-\beta {z}_{{f}_{j}}\right]\hfill \end{array}$$

Now, let us consider the forward model **g** = **Hf** + **ϵ** and assign a Gaussian law with unknown variance
${\upsilon}_{{\in}_{i}}$ to the noise ϵ_{i}, which results in
$p(\in )=\mathcal{N}(g|0,{\mathrm{V}}_{\in})$ with
${\mathrm{V}}_{\in}=\text{diag}[{\upsilon}_{e}]$ with
${\upsilon}_{\in}=[{\upsilon}_{{\in}_{1}},\cdots ,{\upsilon}_{{\in}_{M}}]$, and so:
Let us also note by
${z}_{{\in}_{i}}=1/{\upsilon}_{{e}_{i}}$,
${z}_{\in}=[{z}_{{\in}_{1}},\cdots ,{z}_{{\in}_{M}}]$and
${\mathrm{Z}}_{\in}=\text{diag}\phantom{\rule{0.2em}{0ex}}[{z}_{\in}]={\mathrm{V}}_{\in}^{-1}$ and assign a prior on it
$p({\upsilon}_{{\in}_{i}}|{\alpha}_{{\in}_{0}},{\beta}_{\in}{}_{0})=\mathcal{I}\mathcal{G}({\upsilon}_{{\in}_{i}}|{\alpha}_{{\in}_{0}},{\beta}_{{\in}_{0}})$ or equivalently:

$$p(\mathit{g}|\mathit{f},{\mathbf{\upsilon}}_{\in})=\mathcal{N}(\mathit{g}|\mathit{Hf},{\mathit{V}}_{\in})\propto \mathrm{exp}\left[-\frac{1}{2}(\mathit{g}-\mathit{Hf}){\mathit{V}}_{\in}^{-1}(\mathit{g}-\mathit{Hf})\right].$$

$$p({z}_{{\in}_{i}}|{\alpha}_{{\in}_{0}},{\beta}_{{\in}_{0}})=\mathcal{G}({z}_{{\in}_{i}}|{\alpha}_{{\in}_{0}},{\beta}_{{\in}_{0}})\phantom{\rule{0.2em}{0ex}}and\phantom{\rule{0.2em}{0ex}}p({z}_{\in}|{\alpha}_{{\in}_{0}},{\beta}_{{\in}_{0}})={\displaystyle \prod _{i}\mathcal{G}({z}_{{\in}_{i}}|{\alpha}_{{\in}_{0}},{\beta}_{{\in}_{0}})}.$$

Let us also note
${\upsilon}_{f}=[{\upsilon}_{f{}_{1}},\cdots ,{\upsilon}_{{f}_{N}}]$,
${\mathrm{V}}_{f}=\text{diag}\phantom{\rule{0.2em}{0ex}}[{\upsilon}_{f}]$,
${z}_{{f}_{i}}=1/{\upsilon}_{{f}_{i}}$,
${\mathrm{Z}}_{f}=\text{diag}\phantom{\rule{0.2em}{0ex}}[{z}_{f}]={\mathrm{V}}_{f}^{-1}$ and note:
and finally,

$$p(\mathit{f}|{\upsilon}_{f})={\displaystyle \prod _{j}\mathcal{N}}({f}_{j}|0,{\upsilon}_{{f}_{j}})={\displaystyle \prod _{j}\mathcal{N}({f}_{j}|0,{\upsilon}_{{f}_{j}})}=\mathcal{N}(\mathit{f}|0,{\mathit{V}}_{f})$$

$$p({\upsilon}_{f}|{\alpha}_{{f}_{0}},{\beta}_{{f}_{0}})={\displaystyle \prod _{j}\mathcal{G}({\upsilon}_{{f}_{i}}|{\alpha}_{{f}_{0}},{\beta}_{{f}_{0}})}.$$

Then, we obtain the following expressions for the VBA:
where:
We have implemented these algorithms for many linear inverse problems [102], such as periodic components estimation in time series [122] or computed tomography [123], blind deconvolution [124], blind image separation [125,126] and blind image restoration [89].

$$\{\begin{array}{l}{q}_{1}(\mathit{f}|\tilde{\mathit{\mu}},{\tilde{\mathit{V}}}_{f})=\mathcal{N}(\mathit{f}|\tilde{\mathit{\mu}},\tilde{\mathit{V}})\phantom{\rule{0.2em}{0ex}}\text{with}\phantom{\rule{0.2em}{0ex}}\tilde{\mathit{\mu}}=\tilde{\mathit{V}}{\mathit{H}}^{\prime}\mathit{g},\phantom{\rule{0.2em}{0ex}}\tilde{\mathit{V}}={({\mathit{H}}^{\prime}{\tilde{\mathit{V}}}_{\in}^{-1}\mathit{H}+{\tilde{\mathit{Z}}}_{f})}^{-1};\hfill \\ {q}_{2j}({z}_{{f}_{j}})=\mathcal{G}({z}_{{f}_{j}}|{\tilde{\alpha}}_{j},{\tilde{\beta}}_{j})\phantom{\rule{0.2em}{0ex}}\text{with}\phantom{\rule{0.2em}{0ex}}{\tilde{\alpha}}_{j}={\alpha}_{00}+1/2,{\tilde{\beta}}_{j}={\beta}_{00}+<{f}_{j}^{2}>/2;\hfill \\ {q}_{3}({z}_{{\in}_{i}})=\mathcal{G}(z{}_{{\in}_{i}}|{\tilde{\alpha}}_{{\in}_{i}},{\tilde{\beta}}_{{\in}_{i}})\phantom{\rule{0.2em}{0ex}}\text{with}\phantom{\rule{0.2em}{0ex}}{\tilde{\alpha}}_{{\in}_{i}}={\alpha}_{{\in}_{0}}+(\mathcal{N}+1)/2,\phantom{\rule{0.2em}{0ex}}{\tilde{\beta}}_{{\in}_{i}}={\beta}_{{\in}_{0}}+\frac{1}{2}<{g}_{i}-{[\mathit{Hf}]}_{i}{|}^{2}>;\hfill \end{array}$$

$$\begin{array}{l}<{|{g}_{i}-{[\mathit{Hf}]}_{i}|}^{2}>={|{g}_{i}-\mathit{H}<\mathit{f}>{]}_{i}|}^{2}+{[{\mathit{H}}^{\prime}\tilde{\mathit{V}}\mathit{H}]}_{\mathit{ii}},\\ <\mathit{f}>=\tilde{\mathit{\mu}},\phantom{\rule{0.2em}{0ex}}<\mathit{f}{\mathit{f}}^{\prime}>=\tilde{\mathit{V}}+\tilde{\mathit{\mu}}{\tilde{\mathit{\mu}}}^{\prime},\\ <{f}_{j}^{2}>={[\tilde{\mathit{V}}]}_{\mathit{jj}}+{\tilde{\mu}}_{j}^{2}\end{array}$$

The main conclusions of this paper can be summarized as follows:

- A probability law is a tool for representing our state of knowledge about a quantity.
- The Bayes or Laplace rule is an inference tool for updating our state of knowledge about an inaccessible quantity when another accessible, related quantity is observed.
- Entropy is a measure of information content in a variable with a given probability law.
- The maximum entropy principle can be used to assign a probability law to a quantity when the available information about it is in the form of a limited number of constraints on that probability law.
- Relative entropy and Kullback–Leibler divergence are tools for updating probability laws in the same context.
- When a parametric probability law is assigned to a quantity and we want to measure the amount of information gain about the parameters when some direct observations of that quantity is available, we can use the Fisher information. The structure of the Fisher information geometry in the space of parameters is derived from the relative entropy by a second order Taylor series approximation.
- All of these rules and tools are used currently in different ways in data and signal processing. In this paper, a few examples of the ways these tools are used in data and signal processing problems are presented. One main conclusion is that each of these tools has to be used in appropriate contexts. The example in spectral estimation shows that it is very important to define the problems very clearly at the beginning and to use appropriate tools and interpret the results appropriately.
- The Laplacian or Bayesian inference is the appropriate tool for proposing satisfactory solutions to inverse problems. Indeed, the expression of the posterior probability law represents the combination of the state of the knowledge in the forward model and the data and the state of the knowledge before using the data.
- The Bayesian approach can also easily be used to propose unsupervised methods for the practical application of these methods.
- One of the main limitation of those sophisticated methods is the computational cost. For this, we proposed to use VBA as an alternative to MCMC methods to propose realistic algorithms in huge dimensional inverse problems where we want to estimate an unknown signal (1D), image (2D), volume (3D) or even more (3D + time or 3D + wavelength), etc.

The author would like to thank the reviewers who, by their true review work and their extensive comments and remarks, helped to improve this review paper greatly.

^{†}This paper is an extended version of the paper published in Proceedings of the 34th International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering (MaxEnt 2014), Amboise, France, 21–26 September 2014.

The author declares no conflict of interest.

- Mohammad-Djafari, A. Bayesian or Laplacian inference, entropy and information theory and information geometry in data and signal processing. AIP Conf. Proc
**2014**, 1641, 43–58. [Google Scholar] - Bayes, T. An Essay toward Solving a Problem in the Doctrine of Chances. Philos. Trans
**1763**, 53, 370–418, By the late Rev. Mr. Bayes communicated by Mr. Price, in a Letter to John Canton. [Google Scholar] - De Laplace, P. S. Mémoire sur la probabilité des causes par les évènements. Mémoires de l’Academie Royale des Sciences Presentés par Divers Savan
**1774**, 6, 621–656. [Google Scholar] - Shannon, C. A Mathematical Theory of Communication. Bell Syst. Tech. J
**1948**, 27, 379–423. [Google Scholar] - Hadamard, J. Mémoire sur le problème d’analyse relatif à l’équilibre des plaques élastiques encastrées; Mémoires présentés par divers savants à l’Académie des sciences de l’Institut de France; Imprimerie nationale, 1908. [Google Scholar]
- Jaynes, E.T. Information Theory Statistical Mechanics. Phys. Rev
**1957**, 106, 620–630. [Google Scholar] - Jaynes, E.T. Information Theory and Statistical Mechanics II. Phys. Rev
**1957**, 108, 171–190. [Google Scholar] - Jaynes, E.T. Prior Probabilities. IEEE Trans. Syst. Sci. Cybern
**1968**, 4, 227–241. [Google Scholar] - Kullback, S. Information Theory and Statistics; Wiley: New York, NY, USA, 1959. [Google Scholar]
- Fisher, R. On the mathematical foundations of theoretical statistics. Philos. Trans. R. Stat. Soc. A
**1922**, 222, 309–368. [Google Scholar] - Rao, C. Information and accuracy attainable in the estimation of statistical parameters. Bull. Culcutta Math. Soc
**1945**, 37, 81–91. [Google Scholar] - Sindhwani, V.; Belkin, M.; Niyogi, P. The Geometric basis for Semi-supervised Learning. In Semi-supervised Learning; Chapelle, O., Schölkopf, B., Zien, A., Eds.; MIT press: Cambridge, MA, USA, 2006; pp. 209–226. [Google Scholar]
- Lin, J. Divergence Measures Based on the Shannon Entropy. IEEE Trans. Inf. Theory
**1991**, 37, 145–151. [Google Scholar] - Johnson, O.; Barron, A.R. Fisher Information Inequalities and the Central Limit Theorem. Probab. Theory Relat. Fields
**2004**, 129, 391–409. [Google Scholar] - Berger, J. Statistical Decision Theory and Bayesian Analysis, 2nd ed; Springer-Verlag: New York, NY, USA, 1985. [Google Scholar]
- Gelman, A.; Carlin, J.B.; Stern, H.S.; Rubin, D.B. Bayesian Data Analysis, 2nd ed; Chapman & Hall/CRC Texts in Statistical Science; Chapman and Hall/CRC: Boca Raton, FL, USA, 2003. [Google Scholar]
- Skilling, J. Nested Sampling. In Bayesian Inference and Maximum Entropy Methods in Science and Engineering; Proceedings of 24th International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering, Garching, Germany, 25–30 July 2004, Fischer, R., Preuss, R., Toussaint, U.V., Eds.; pp. 395–405.
- Metropolis, N.; Rosenbluth, A.W.; Rosenbluth, M.N.; Teller, A.H.; Teller, E. Equation of State Calculations by Fast Computing Machines. J. Chem. Phys
**1953**, 21, 1087–1092. [Google Scholar] - Hastings, W.K. Monte Carlo Sampling Methods using Markov Chains and their Applications. Biometrika
**1970**, 57, 97–109. [Google Scholar] - Gelfand, A.E.; Smith, A.F.M. Sampling-Based Approaches to Calculating Marginal Densities. J. Am. Stat. Assoc
**1990**, 85, 398–409. [Google Scholar] - Gilks, W.R.; Richardson, S.; Spiegelhalter, D.J. Introducing Markov Chain Monte Carlo. In Markov Chain Monte Carlo in Practice; Gilks, W.R., Richardson, S., Spiegelhalter, D.J., Eds.; Chapman and Hall: London, UK, 1996; pp. 1–19. [Google Scholar]
- Gilks, W.R. Strategies for Improving MCMC. In Markov Chain Monte Carlo in Practice; Gilks, W.R., Richardson, S., Spiegelhalter, D.J., Eds.; Chapman and Hall: London, UK, 1996; pp. 89–114. [Google Scholar]
- Roberts, G.O. Markov Chain Concepts Related to Sampling Algorithms. In Markov Chain Monte Carlo in Practice; Gilks, W.R., Richardson, S., Spiegelhalter, D.J., Eds.; Chapman and Hall: London, UK, 1996; pp. 45–57. [Google Scholar]
- Tanner, M.A. Tools for Statistical Inference: Methods for the Exploration of Posterior Distributions and Likelihood Functions; Springer series in Statistics; Springer: New York, NY, USA, 1996. [Google Scholar]
- Djurić, P.M.; Godsill, S.J. (Eds.) Special Issue on Monte Carlo Methods for Statistical Signal Processing; IEEE: New York, NY, USA, 2002.
- Andrieu, C.; de Freitas, N.; Doucet, A.; Jordan, M.I. An Introduction to MCMC for Machine Learning. Mach. Learn
**2003**, 50, 5–43. [Google Scholar] - Clausius, R. On the Motive Power of Heat, and on the Laws Which Can be Deduced From it for the Theory of Heat; Poggendorff’s Annalen der Physick, LXXIX, Dover Reprint: New York, NY, USA, 1850; ISBN ISBN 0-486-59065-8. [Google Scholar]
- Caticha, A. Maximum Entropy, fluctuations and priors.
- Giffin, A.; Caticha, A. Updating Probabilities with Data and Moments.
- Caticha, A.; Preuss, R. Maximum Entropy and Bayesian Data Analysis: Entropic Priors Distributions. Phys. Rev. E
**2004**, 70, 046127. [Google Scholar] - Akaike, H. On Entropy Maximization Principle. In Applications of Statistics; Krishnaiah, P.R., Ed.; North-Holland: Amsterdam, The Netherlands, 1977; pp. 27–41. [Google Scholar]
- Agmon, N.; Alhassid, Y.; Levine, D. An Algorithm for Finding the Distribution of Maximal Entropy. J. Comput. Phys
**1979**, 30, 250–258. [Google Scholar] - Jaynes, E.T. Where do we go from here? In Maximum-Entropy and Bayesian Methods in Inverse Problems; Smith, C.R., Grandy, W.T., Jr, Eds.; Springer: Dordrecht, The Netherlands, 1985; pp. 21–58. [Google Scholar]
- Borwein, J.M.; Lewis, A.S. Duality relationships for entropy-like minimization problems. SIAM J. Control Optim
**1991**, 29, 325–338. [Google Scholar] - Elfwing, T. On some Methods for Entropy Maximization and Matrix Scaling. Linear Algebra Appl
**1980**, 34, 321–339. [Google Scholar] - Eriksson, J. A note on Solution of Large Sparse Maximum Entropy Problems with Linear Equality Constraints. Math. Program
**1980**, 18, 146–154. [Google Scholar] - Erlander, S. Entropy in linear programs. Math. Program
**1981**, 21, 137–151. [Google Scholar] - Jaynes, E.T. On the Rationale of Maximum-Entropy Methods. Proc. IEEE
**1982**, 70, 939–952. [Google Scholar] - Shore, J.E.; Johnson, R.W. Properties of Cross-Entropy Minimization. IEEE Trans. Inf. Theory
**1981**, 27, 472–482. [Google Scholar] - Mohammad-Djafari, A. Maximum d’entropie et problèmes inverses en imagerie. Traitement Signal
**1994**, 11, 87–116. [Google Scholar] - Bercher, J. Développement de critères de nature entropique pour la résolution des problèmes inverses linéaires. Ph.D. Thesis, Université de Paris–Sud, Orsay, France, 1995. [Google Scholar]
- Le Besnerais, G. Méthode du maximum d’entropie sur la moyenne, critère de reconstruction d’image et synthèse d’ouverture en radio astronomie. In Ph.D. Thesis; Université de Paris-Sud: Orsay, France, 1993. [Google Scholar]
- Caticha, A.; Giffin, A. Updating Probabilities. [CrossRef]
- Caticha, A. Entropic Inference.
- Costa, S.I.R.; Santos, S.A.; Strapasson, J.E. Fisher information distance: A geometrical reading
**2012**. arXiv: 1210.2354. - Rissanen, J. Fisher Information and Stochastic Complexity. IEEE Trans. Inf. Theory
**1996**, 42, 40–47. [Google Scholar] - Shimizu, R. On Fisher’s amount of information for location family. In A Modern Course on Statistical Distributions in Scientific Work; D. Reidel Dordrecht: The Netherlands, 1975; Volume 3, pp. 305–312. [Google Scholar]
- Nielsen, F.; Nock, R. Sided and Symmetrized Bregman Centroids. IEEE Trans. Inf. Theory
**2009**, 55, 2048–2059. [Google Scholar] - Rasmussen, C.E.; Williams, C.K.I. Gaussian Processes for Machine Learning; MIT Press: Cambridge, MA, USA, 2006. [Google Scholar]
- Schroeder, M.R. Linear prediction, entropy and signal analysis. IEEE ASSP Mag
**1984**, 1, 3–11. [Google Scholar] - Itakura, F.; Saito, S. A Statistical Method for Estimation of Speech Spectral Density and Formant Frequencies. Electron. Commun. Jpn
**1970**, 53-A, 36–43. [Google Scholar] - Kitagawa, G.; Gersch, W. Smoothness Priors Analysis of Time Series; Lecture Notes in Statistics; Volume 116, Springer: New York, NY, USA, 1996. [Google Scholar]
- Rue, H.; Held, L. Gaussian Markov Random Fields: Theory and Applications; CRC Press: New York, NY, USA, 2005. [Google Scholar]
- Amari, S.; Cichocki, A.; Yang, H.H. A new learning algorithm for blind source separation. 757–763.
- Amari, S. Neural learning in structured parameter spaces—Natural Riemannian gradient. 127–133.
- Amari, S. Natural gradient works efficiently in learning. Neural Comput
**1998**, 10, 251–276. [Google Scholar] - Knuth, K.H. Bayesian source separation and localization. SPIE Proc
**1998**, 3459. [Google Scholar] [CrossRef] - Knuth, K.H. A Bayesian approach to source separation. 283–288.
- Attias, H. Independent Factor Analysis. Neural Comput
**1999**, 11, 803–851. [Google Scholar] - Mohammad-Djafari, A. A Bayesian approach to source separation. 221–244.
- Choudrey, R.A.; Roberts, S. Variational Bayesian Mixture of Independent Component Analysers for Finding Self-Similar Areas in Images. 107–112.
- Lopes, H.F.; West, M. Bayesian Model Assessment in Factor Analysis. Statsinica
**2004**, 14, 41–67. [Google Scholar] - Ichir, M.; Mohammad-Djafari, A. Bayesian Blind Source Separation of Positive Non Stationary Sources. 493–500.
- Mohammad-Djafari, A. Bayesian Source Separation: Beyond PCA and ICA.
- Comon, P.; Jutten, C. (Eds.) Handbook of Blind Source Separation: Independent Component Analysis and Applications; Academic Press: Burlington, MA, USA, 2010.
- Yuan, M.; Lin, Y. Model selection and estimation in the Gaussian graphical model. Biometrika
**2007**, 94, 19–35. [Google Scholar] - Fitzgerald, W. Markov Chain Monte Carlo methods with Applications to Signal Processing. Signal Process
**2001**, 81, 3–18. [Google Scholar] - Matsuoka, T.; Ulrych, T. Information theory measures with application to model identification. IEEE Trans. Acoust. Speech Signal Process
**1986**, 34, 511–517. [Google Scholar] - Bretthorst, G.L. Bayesian Model Selection: Examples Relevant to NMR. In Maximum Entropy and Bayesian Methods; Springer: Dordrecht, The Netherlands, 1989; pp. 377–388. [Google Scholar]
- Gelfand, A.E.; Dey, D.K. Bayesian model choice: Asymptotics and exact calculations. J. R. Stat. Soc. Ser. B
**1994**, 56, 501–514. [Google Scholar] - Mohammad-Djafari, A. Model selection for inverse problems: Best choice of basis function and model order selection.
- Clyde, M.A.; Berger, J.O.; Bullard, F.; Ford, E.B.; Jefferys, W.H.; Luo, R.; Paulo, R.; Loredo, T. Current Challenges in Bayesian Model Choice. 71, 224–240.
- Wyse, J.; Friel, N. Block clustering with collapsed latent block models. Stat. Comput
**2012**, 22, 415–428. [Google Scholar] - Giovannelli, J.F.; Giremus, A. Bayesian noise model selection and system identification based on approximation of the evidence. 125–128.
- Akaike, H. A new look at the statistical model identification. IEEE Trans. Automat. Control
**1974**, AC-19, 716–723. [Google Scholar] - Akaike, H. Power spectrum estimation through autoregressive model fitting. Ann. Inst. Stat. Math
**1969**, 21, 407–419. [Google Scholar] - Farrier, D. Jaynes’ principle and maximum entropy spectral estimation. IEEE Trans. Acoust. Speech Signal Process
**1984**, 32, 1176–1183. [Google Scholar] - Wax, M. Detection and Estimation of Superimposed Signals. Ph.D. Thesis, Standford University, CA, USA, March 1985. [Google Scholar]
- Burg, J.P. Maximum Entropy Spectral Analysis.
- McClellan, J.H. Multidimensional spectral estimation. Proc. IEEE
**1982**, 70, 1029–1039. [Google Scholar] - Lang, S.; McClellan, J.H. Multidimensional MEM spectral estimation. IEEE Trans. Acoust. Speech Signal Process
**1982**, 30, 880–887. [Google Scholar] - Johnson, R.; Shore, J. Which is Better Entropy Expression for Speech Processing:-SlogS or logS? IEEE Trans. Acoust. Speech Signal Process
**1984**, ASSP-32, 129–137. [Google Scholar] - Wester, R.; Tummala, M.; Therrien, C. Multidimensional Autoregressive Spectral Estimation Using Iterative Methods. 1. [CrossRef]
- Picinbono, B.; Barret, M. Nouvelle présentation de la méthode du maximum d’entropie. Traitement Signal
**1990**, 7, 153–158. [Google Scholar] - Borwein, J.M.; Lewis, A.S. Convergence of best entropy estimates. SIAM J. Optim
**1991**, 1, 191–205. [Google Scholar] - Mohammad-Djafari, A. (Ed.) Inverse Problems in Vision and 3D Tomography; digital signal and image processing series; ISTE: London, UK; Wiley: Hoboken, NJ, USA, 2010.
- Mohammad-Djafari, A.; Demoment, G. Tomographie de diffraction and synthèse de Fourier à maximum d’entropie. Rev. Phys. Appl. (Paris)
**1987**, 22, 153–167. [Google Scholar] - Féron, O.; Chama, Z.; Mohammad-Djafari, A. Reconstruction of piecewise homogeneous images from partial knowledge of their Fourier transform. 68–75.
- Ayasso, H.; Mohammad-Djafari, A. Joint NDT Image Restoration and Segmentation Using Gauss–Markov–Potts Prior Models and Variational Bayesian Computation. IEEE Trans. Image Process
**2010**, 19, 2265–2277. [Google Scholar] - Ayasso, H.; DuchÃłne, B.; Mohammad-Djafari, A. Bayesian inversion for optical diffraction tomography. J. Mod. Opt
**2010**, 57, 765–776. [Google Scholar] - Burch, S.; Gull, S.F.; Skilling, J. Image Restoration by a Powerful Maximum Entropy Method. Comput. Vis. Graph. Image Process
**1983**, 23, 113–128. [Google Scholar] - Gull, S.F.; Skilling, J. Maximum entropy method in image processing. IEE Proc. F
**1984**, 131, 646–659. [Google Scholar] - Gull, S.F. Developments in maximum entropy data analysis. In Maximum Entropy and Bayesian Methods; Skilling, J., Ed.; Springer: Dordrecht, The Netherlands, 1989; pp. 53–71. [Google Scholar]
- Jones, L.K.; Byrne, C.L. General entropy criteria for inverse problems with application to data compression, pattern classification and cluster analysis. IEEE Trans. Inf. Theory
**1990**, 36, 23–30. [Google Scholar] - Macaulay, V.A.; Buck, B. Linear inversion by the method of maximum entropy. Inverse Probl
**1989**, 5. [Google Scholar] [CrossRef] - Rue, H.; Martino, S. Approximate Bayesian inference for hierarchical Gaussian Markov random field models. J. Stat. Plan. Inference
**2007**, 137, 3177–3192. [Google Scholar] - Wilkinson, R. Approximate Bayesian computation (ABC) gives exact results under the assumption of model error
**2009**. arXiv:0811.3355. - Rue, H.; Martino, S.; Chopin, N. Approximate Bayesian Inference for Latent Gaussian Models Using Integrated Nested Laplace Approximations. J. R. Stat. Soc. Ser. B
**2009**, 71, 319–392. [Google Scholar] - Fearnhead, P.; Prangle, D. Constructing Summary Statistics for Approximate Bayesian Computation: Semi-automatic ABC
**2011**. arxiv:1004.1112v2. - Turner, B.M.; van Zandt, T. A tutorial on approximate Bayesian computation. J. Math. Psych
**2012**, 56, 69–85. [Google Scholar] - MacKay, D.J.C. A Practical Bayesian Framework for Backpropagation Networks. Neural Comput
**1992**, 4, 448–472. [Google Scholar] - Mohammad-Djafari, A. Variational Bayesian Approximation for Linear Inverse Problems with a hierarchical prior models. 8085, 669–676.
- Likas, C.L.; Galatsanos, N.P. A Variational Approach For Bayesian Blind Image Deconvolution. IEEE Trans. Signal Process
**2004**, 52, 2222–2233. [Google Scholar] - Beal, M.; Ghahramani, Z. Variational Bayesian learning of directed graphical models with hidden variables. Bayesian Stat
**2006**, 1, 793–832. [Google Scholar] - Kim, H.; Ghahramani, Z. Bayesian Gaussian Process Classification with the EM-EP Algorithm. IEEE Trans. Pattern Anal. Mach. Intell
**2006**, 28, 1948–1959. [Google Scholar] - Jordan, M.I.; Ghahramani, Z.; Jaakkola, T.S.; Saul, L.K. An introduction to variational methods for graphical models. Mach. Learn
**2006**, 37, 183–233. [Google Scholar] - Forbes, F.; Fort, G. Combining Monte Carlo and Mean-Field-Like Methods for Inference in Hidden Markov Random Fields. IEEE Trans. Image Process
**2007**, 16, 824–837. [Google Scholar] - Dempster, A.P.; Laird, N.M.; Rubin, D.B. Maximum Likelihood from incomplete data via the EM algorithm. J. R. Stat. Soc. (B)
**1977**, 39, 1–38. [Google Scholar] - Miller, M.I.; Snyder, D.L. The Role of Likelihood and Entropy in Incomplete-Data Problems: Applications to Estimating Point-Process Intensities and Toeplitz Constrained Covariances. Proc. IEEE
**1987**, 75, 892–907. [Google Scholar] - Snoussi, H.; Mohammad-Djafari, A. Information geometry of Prior Selection.
- Mohammad-Djafari, A. Approche variationnelle pour le calcul bayésien dans les problémes inverses en imagerie
**2009**. arXiv:0904.4148. - Beal, M. Variational Algorithms for Approximate Bayesian Inference. In Ph.D. Thesis; Gatsby Computational Neuroscience Unit, University College London: UK, 2003. [Google Scholar]
- Winn, J.; Bishop, C.M.; Jaakkola, T. Variational message passing. J. Mach. Learn. Res
**2005**, 6, 661–694. [Google Scholar] - Chatzis, S.; Varvarigou, T. Factor Analysis Latent Subspace Modeling and Robust Fuzzy Clustering Using t-DistributionsClassification of binary random Patterns. IEEE Trans. Fuzzy Syst
**2009**, 17, 505–517. [Google Scholar] - Park, T.; Casella, G. The Bayesian Lasso. J. Am. Stat. Assoc
**2008**, 103, 681–686. [Google Scholar] - Mohammad-Djafari, A. A variational Bayesian algorithm for inverse problem of computed tomography. In Mathematical Methods in Biomedical Imaging and Intensity-Modulated Radiation Therapy (IMRT); Censor, Y., Jiang, M., Louis, A.K., Eds.; Publications of the Scuola Normale Superiore/CRM Series; Edizioni della Normale: Rome, Italy, 2008; pp. 231–252. [Google Scholar]
- Mohammad-Djafari, A.; Ayasso, H. Variational Bayes and mean field approximations for Markov field unsupervised estimation. 1–6.
- Tipping, M.E. Sparse Bayesian learning and the relevance vector machine. J. Mach. Learn. Res
**2001**, 1, 211–244. [Google Scholar] - He, L.; Chen, H.; Carin, L. Tree-Structured Compressive Sensing With Variational Bayesian Analysis. IEEE Signal Process. Lett
**2010**, 17, 233–236. [Google Scholar] - Fraysse, A.; Rodet, T. A gradient-like variational Bayesian algorithm, Proceedings of 2011 IEEE Conference on Statistical Signal Processing Workshop (SSP), Nice France, 28–30 June 2011; pp. 605–608.
- Johnson, V.E. On Numerical Aspects of Bayesian Model Selection in High and Ultrahigh-dimensional Settings. Bayesian Anal
**2013**, 8, 741–758. [Google Scholar] - Dumitru, M.; Mohammad-Djafari, A. Estimating the periodic components of a biomedical signal through inverse problem modeling and Bayesian inference with sparsity enforcing prior. AIP Conf. Proc
**2015**, 1641, 548–555. [Google Scholar] - Wang, L.; Gac, N.; Mohammad-Djafari, A. Bayesian 3D X-ray computed tomography image reconstruction with a scaled Gaussian mixture prior model. AIP Conf. Proc
**2015**, 1641, 556–563. [Google Scholar] - Mohammad-Djafari, A. Bayesian Blind Deconvolution of Images Comparing JMAP, EM and VBA with a Student-t a priori Model. 98–103.
- Su, F.; Mohammad-Djafari, A. An Hierarchical Markov Random Field Model for Bayesian Blind Image Separation.
- Su, F.; Cai, S.; Mohammad-Djafari, A. Bayesian blind separation of mixed text patterns. 1373–1378.

$\mu (u)\propto \mathrm{exp}\left[{-}_{2}^{1}{\sum}_{j}{u}_{j}^{2}\right]$ | $\widehat{\mathit{f}}={\mathit{H}}^{\prime}\mathbf{\lambda}$ | $\widehat{\mathit{f}}={\mathit{H}}^{\prime}{(\mathit{H}{\mathit{H}}^{\prime})}^{-1}\mathit{g}$ |

$\mu (u)\propto \mathrm{exp}\left[-{\sum}_{j}|{u}_{j}|\right]$ | $\widehat{\mathit{f}}=1./({\mathit{H}}^{\prime}\mathbf{\lambda}\pm 1)$ | $\mathit{H}\widehat{\mathit{f}}=\mathit{g}$ |

$\mu (u)\propto \mathrm{exp}\left[-{\sum}_{j}{u}_{j}^{\mathrm{\alpha}-1}\mathrm{exp}\left[-\beta {u}_{j}\right]\right],\phantom{\rule{0.2em}{0ex}}{u}_{j}>0$ | $\widehat{\mathit{f}}=\alpha 1./({\mathit{H}}^{\prime}\mathbf{\lambda}+\beta 1)$ | $\mathit{H}\widehat{\mathit{f}}=\mathit{g}$ |

JMPA | BEM | VBA |
---|---|---|

$q(\mathit{f})=\delta (\mathit{f}-\tilde{\mathit{f}})$ | $q(\mathit{f})=N(\mathit{f}|\tilde{\mathit{f}},\tilde{\mathit{V}})$ | $q(\mathit{f})=N(\mathit{f}|\tilde{\mathit{f}},\tilde{\mathit{V}})$ |

$\tilde{\mathit{V}}={({\mathit{H}}^{\prime}\mathit{H}+\tilde{\lambda}\mathit{I})}^{-1}$ | $\tilde{\mathit{V}}={({\mathit{H}}^{\prime}\mathit{H}+\tilde{\lambda}\mathit{I})}^{-1}$ | $\tilde{\mathit{V}}={({\mathit{H}}^{\prime}\mathit{H}+\tilde{\lambda}\mathit{I})}^{-1}$ |

$\tilde{\mathit{f}}=\tilde{\mathit{V}}{\mathit{H}}^{\prime}\mathit{g}$ | $\tilde{\mathit{f}}=\tilde{\mathit{V}}{\mathit{H}}^{\prime}\mathit{g}$ | $\tilde{\mathit{f}}=\tilde{\mathit{V}}{\mathit{H}}^{\prime}\mathit{g}$ |

$q({\theta}_{1})=\delta ({\theta}_{1}-{\tilde{\theta}}_{1})$ | $q({\theta}_{1})=\delta ({\theta}_{1}-{\tilde{\theta}}_{1})$ | $q({\theta}_{1})=\mathcal{G}({\theta}_{1}|{\tilde{\alpha}}_{1},{\tilde{\beta}}_{1})$ |

${\tilde{\alpha}}_{1}=({\alpha}_{10}-1)+\frac{\mathrm{M}}{2}$ | ${\tilde{\alpha}}_{1}=({\alpha}_{10}-1)+\frac{\mathrm{M}}{2}$ | ${\tilde{\alpha}}_{1}=({\alpha}_{10}-1)+\frac{\mathrm{M}}{2}$ |

${\tilde{\beta}}_{1}={\beta}_{10}+\frac{1}{2}\Vert \mathit{g}-\mathit{Hf}\Vert {}_{2}^{2}$ | ${\tilde{\beta}}_{1}={\beta}_{10}+\frac{1}{2}<\Vert \mathit{g}-\mathit{Hf}\Vert {}_{2}^{2}>$ | ${\tilde{\beta}}_{1}={\beta}_{10}+\frac{1}{2}<\Vert \mathit{g}-\mathit{Hf}\Vert {}_{2}^{2}>$ |

${\tilde{\theta}}_{1}=\frac{{\tilde{\alpha}}_{1}}{{\beta}_{1}}$ | ${\tilde{\theta}}_{1}=\frac{{\tilde{\alpha}}_{1}}{{\tilde{\beta}}_{1}}$ | ${\tilde{\theta}}_{1}=\frac{{\tilde{\alpha}}_{1}}{{\tilde{\beta}}_{1}}$ |

$q({\theta}_{2})=\delta ({\theta}_{2}-{\tilde{\theta}}_{2})$ | $q({\theta}_{2})=\delta ({\theta}_{2}-{\tilde{\theta}}_{2})$ | $q({\theta}_{2})=\mathcal{G}({\theta}_{2}|{\tilde{\alpha}}_{2},{\tilde{\beta}}_{2})$ |

${\tilde{\alpha}}_{2}=({\alpha}_{20}-1)+\frac{\mathrm{M}}{2}$ | ${\tilde{\alpha}}_{2}=({\alpha}_{20}-1)+\frac{\mathrm{M}}{2}$ | ${\tilde{\alpha}}_{2}=({\alpha}_{20}-1)+\frac{\mathcal{N}}{2}$ |

${\tilde{\beta}}_{2}={\beta}_{10}+\frac{1}{2}\Vert \mathit{f}\Vert {}_{2}^{2}$ | ${\tilde{\beta}}_{2}={\beta}_{20}+\frac{1}{2}<\Vert \mathit{f}\Vert {}_{2}^{2}>$ | ${\tilde{\beta}}_{2}={\beta}_{20}+\frac{1}{2}<\Vert \mathit{f}\Vert {}_{2}^{2}>$ |

${\tilde{\theta}}_{2}=\frac{{\tilde{\alpha}}_{2}}{{\tilde{\beta}}_{2}}$ | ${\tilde{\theta}}_{2}=\frac{{\tilde{\alpha}}_{2}}{{\tilde{\beta}}_{2}}$ | ${\tilde{\theta}}_{2}=\frac{{\tilde{\alpha}}_{2}}{{\tilde{\beta}}_{2}}$ |

$\tilde{\lambda}=\frac{{\tilde{\theta}}_{2}}{{\tilde{\theta}}_{1}}$ | $\tilde{\lambda}=\frac{{\tilde{\theta}}_{2}}{{\tilde{\theta}}_{1}}$ | $\tilde{\lambda}=\frac{{\tilde{\theta}}_{2}}{{\tilde{\theta}}_{1}}$ |

© 2015 by the authors; licensee MDPI, Basel, Switzerland This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution license (http://creativecommons.org/licenses/by/4.0/).