Wavelet-based hyperparameter estimation for solving inverse problems

We consider the problem of image restoration/reconstruction where the acquisition system is modeled by a linear operator with additive Gaussian noise. A variational approach is adopted for image inversion in order to compute a restored/reconstructed image, consisting in minimizing a convex criterion composed of two parts: a data fidelity term (e.g. quadratic) and a regularization term (e.g. ℓ1-norm) expressed in the wavelet domain. The purpose of this paper is to estimate the regularization hyperparameters (one per subband) based on a Maximum Likelihood (ML) estimator, only knowing the observed data. A difficult task in such estimation is to compute the expectation according to the a posteriori probability as there is no analytical form. This expectation must be approximated numerically by sampling the distribution. However, sampling the a posteriori distribution is a difficult task because of pixel interactions introduced by the linear operator (image acquisition) in the same time as the wavelet transform (regularization). Moreover, the possible different nature (ℓ2, ℓ1-norm ...) of the fidelity and regularization terms does not allow to easily process them simultaneously. We show that both operators can be separated using an auxiliary (hidden) variable and splitting the a posteriori probability in two parts which are sampled alternately using MCMC (Gibbs sampling and Metropolis-Hastings). We show the equivalence between both formulations of the a posteriori distribution. Then a gradient method is used to estimate the hyperparameters. Simulation results demonstrate the good performance and behavior of the proposed approach