Enhanced sampling schemes for MCMC based blind Bernoulli–Gaussian deconvolution

Abstract

This paper proposes and compares two new sampling schemes for sparse deconvolution using a Bernoulli–Gaussian model. To tackle such a deconvolution problem in a blind and unsupervised context, the Markov Chain Monte Carlo (MCMC) framework is usually adopted, and the chosen sampling scheme is most often the Gibbs sampler. However, such a sampling scheme fails to explore the state space efficiently. Our first alternative, the K-tuple Gibbs sampler, is simply a grouped Gibbs sampler. The second one, called partially marginalized sampler, is obtained by integrating the Gaussian amplitudes out of the target distribution. While the mathematical validity of the first scheme is obvious as a particular instance of the Gibbs sampler, a more detailed analysis is provided to prove the validity of the second scheme.

For both methods, optimized implementations are proposed in terms of computation and storage cost. Finally, simulation results validate both schemes as more efficient in terms of convergence time compared with the plain Gibbs sampler. Benchmark sequence simulations show that the partially marginalized sampler takes fewer iterations to converge than the K-tuple Gibbs sampler. However, its computation load per iteration grows almost quadratically with respect to the data length, while it only grows linearly for the K-tuple Gibbs sampler.

Introduction

The problem of the restoration of a sparse spike train $\mathit{x}$ distorted by a linear system $\mathit{h}$ and corrupted by noise $\epsilon$ such that $\mathit{z}=\mathit{x}*\mathit{h}+\epsilon$, arises in many fields such as seismic exploration [1], [2] and astronomy [3].

In this paper, we adopt a Bernoulli–Gaussian (BG) model for the spike train $\mathit{x}$, following [2] and many subsequent contributions such as [1], [4]. A BG signal is an independent, identically distributed (iid) process defined in two stages. Firstly, the sparse nature of the spikes is governed by the Bernoulli law:

$$P(\mathit{q})={\lambda }^L(1-\lambda {)}^{M-L}$$

with the Bernoulli sequence $\mathit{q}=[q_1,\dots ,q_M{]}^{\mathrm{t}}$ (a binary sequence of length M) and $L={\sum }_{m=1}^M{q}_m$, the number of non-zero entries of $\mathit{q}$. Secondly, amplitudes $\mathit{x}=[x_1,\dots ,x_M{]}^{\mathrm{t}}$ are assumed iid zero-mean Gaussian conditionally to $\mathit{q}$:

$$\mathit{x}|\mathit{q}\sim \mathcal{N}(\mathbf{0},{\sigma }_x^2\mathrm{diag}(\mathit{q}))\text{,}$$

where $\mathrm{diag}(\mathit{q})$ denotes a diagonal matrix whose diagonal is $\mathit{q}$.

The MCMC approach [5], [6] is a powerful numerical tool, appropriate to solve complex inference problems such as blind deconvolution. In the field of blind BG deconvolution, Cheng et al. pioneered the introduction of MCMC methods [1]. They proposed to rely on a plain Gibbs sampler, i.e., with a site-by-site updating scheme for the spike train, for which their algorithm constitutes a simple and canonical example. However, simulation results indicate that it lacks reliability: from different initial conditions, significantly different estimations are obtained, even after a considerable number of iterations.

The recent contribution of [7] already identified a convergence issue linked to time-shift ambiguities, and proposed an efficient way to solve it. In addition, the scale ambiguities are treated by [8] by proposing a scale re-sampling step to accelerate the convergence rate of the Markov chain. In this study, we point out another source of inefficiency, unrelated to the above-mentioned ambiguities: instead of exploring the 2_M configurations of $\mathit{q}$ at an acceptable speed, the Gibbs sampler tends to get stuck for many iterations around some particular configurations of $\mathit{q}$, often corresponding to local optimal configurations of $\mathit{x}$ of the posterior distribution as illustrated by the example in Section 2.2. This conclusion meets Bourguignon and Carfantan's analysis [3]: the Markov chain equilibrates rapidly around a mode (i.e., a local optimal configuration), but takes a long time to move from mode to mode.

In order to make up for this deficiency, our first proposition is to adopt a grouped Gibbs sampler [5], called a K-tuple Gibbs sampler, where blocks of K adjacent BG variables (q_i,…,q_i+K−1) and their associated amplitudes (x_i,…,x_i+K−1) are jointly sampled. A comparable idea first appeared in [9], and more recently in [3], under the form of a deterministic iteration aiming at producing an increase of the posterior probability. We then propose a second solution based on sampling the posterior distribution marginally to the amplitudes $\mathit{x}$ [10]. A comparable idea is found in statistical signal segmentation [11] where some hyper-parameters are partially marginalized. In fact, as Liu pointed out in [5], completely integrating out some components (the Gaussian amplitudes $\mathit{x}$ in our case), leads to a more efficient sampling scheme called collapsed Gibbs sampling. However, a plain collapsed Gibbs sampling on the marginal posterior distribution involves hardly tractable sampling steps. In particular, it is all but simple to sample $\mathit{h}$ conditionally on $(\mathit{z},\mathit{q})$ and marginally with respect to $\mathit{x}$. Our scheme solves this problem by combining a step that samples $\mathit{q}$ marginally with respect to $\mathit{x}$ and other sampling steps involving $\mathit{x}$. Such a partially marginalized sampler is fully valid from the mathematical point of view. In this paper, we show that it can be interpreted as a plain Gibbs sampler with a particular scanning order of the variables.

Simulation tests on toy examples and on the so-called Mendel's sequence [2] confirm the efficiency of both methods in terms of computation time before convergence of the Markov chain. Further analysis shows the data length as a criterion to choose between the two proposed methods.

This paper is organized as follows. In Section 2, after a brief formulation of the blind BG deconvolution problem, the Gibbs sampler of the joint posterior distribution [1] is presented, and an example illustrates its inefficiency as regards the sampling of $\mathit{q}$. 3 Generalized Gibbs on K-tuple variables, 4 Partially marginalized Gibbs sampler, respectively, introduce the generalized K-tuple Gibbs sampler and the partially marginalized sampler. In both cases, a toy example is used to evaluate the capability of the sampler to escape from local optimal configurations, and implementation issues are carefully dealt with.

Finally, simulation results are presented in Section 5 to compare the efficiency of the sampling schemes according to Brooks and Gelman's convergence diagnostic [12], and conclusions are drawn in Section 6.