Advances in sampling schemes for MCMC based blind Bernoulli-Gaussian deconvolution
D. Ge\corauthref{cor}
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. %Two more efficient alternatives are studied and compared in this paper: one is derived from Bourguignon and Carfantan's~\cite{Bourguignon05} idea and yields a generalized $K$-tuple Gibbs sampler; The second one, called partially marginalized sampler, is obtained by integrating the Gaussian amplitudes out of the target distribution. %Unfortunately, for the latter case some of the sampling steps then become intractable. Therefore, our solution mixes steps in which the amplitudes are sampled with others where they are integrated out. To justify the mathematical validity of the proposed partially marginal sampler, the invariant condition is shown to hold first, and then we provide an alternative demonstration using the equivalent standard Gibbs sampler. Another contribution of this article consists of an enhanced numerical implementation for the marginal distribution computation, in terms of both efficiency and storage. 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. %In order to release samples from suboptimal convergence, a general principle is proposed to produce less correlated samples, thus making the algorithm more stable and efficient in terms of convergence rate. Inspired from the classical SMLR scheme studied by Kormylo (1982), appropriate modifications on the basic MCMC methods are proposed and then validated by simulation results that exhibit satisfactory behaviors.
alphabet.tex
% Document LaTeX        GLB (25/05/94)
% Modification :        Gio (15/09/95)
%                       JI  (format Latex 2e, taille adaptative 22/03/99)
%                           (fonte blackboard 27/08/99)

% Objet : Definitions compactes d'alphabets pour ecriture en Math
% b pour bold italique : \Ab
% v pour bold droit : \Av
% D pour droit : \AD
% c pour cursif : \Ac
% bb pour blackboard : \Abb
% Capitale pour capitale : \Lambda

\typeout{Fichier 'alphabet.tex' -- LSS 1994/95/99}
\typeout{(necessite les packages xspace, amsmath, bbm)}
\typeout{}

\RequirePackage{amsmath}
\RequirePackage{xspace}
\RequirePackage{bbm}
%\RequirePackage{amssymb}
%\RequirePackage{bbold}
%\RequirePackage{mathbbol}
%\RequirePackage{doublestroke}

%-------------------------------------------

% DEFINITIONS PRELIMINAIRES 

\def\XS{\xspace}
\DeclareMathAlphabet{\mathb}{OML}{cmm}{b}{it}
\def\sbm#1{\ensuremath{\mathb{#1}}}                % Style gras italique (necessite amsmath)           
%
\def\sbmm#1{\ensuremath{\boldsymbol{#1}}}          % Style gras italique (necessite amsmath)           
\def\sdm#1{\ensuremath{\mathrm{#1}}}               % Style droit en math
\def\sbv#1{\ensuremath{\mathbf{#1}}}               % Style gras droit
\def\scu#1{\ensuremath{\mathcal{#1\XS}}}           % Style cursif
\def\scb#1{\ensuremath{\boldsymbol{\mathcal{#1}}}} % Style gras cursif
\def\sbl#1{\ensuremath{\mathbbm{#1}}}              % Style blackboard (necessite bbm)
%\def\sbl#1{\ensuremath{\mathbb{#1}}}

%-------------------------------------------

% ALPHABET GRAS ITALIQUE, taille adaptative

\def\Ab{{\sbm{A}}\XS}  \def\ab{{\sbm{a}}\XS}
\def\Bb{{\sbm{B}}\XS}  \def\bb{{\sbm{b}}\XS}
\def\Cb{{\sbm{C}}\XS}  \def\cb{{\sbm{c}}\XS}
\def\Db{{\sbm{D}}\XS}  \def\db{{\sbm{d}}\XS}
\def\Eb{{\sbm{E}}\XS}  \def\eb{{\sbm{e}}\XS}
\def\Fb{{\sbm{F}}\XS}  \def\fb{{\sbm{f}}\XS}
\def\Gb{{\sbm{G}}\XS}  \def\gb{{\sbm{g}}\XS}
\def\Hb{{\sbm{H}}\XS}  \def\hb{{\sbm{h}}\XS}
\def\Ib{{\sbm{I}}\XS}  \def\ib{{\sbm{i}}\XS}
\def\Jb{{\sbm{J}}\XS}  \def\jb{{\sbm{j}}\XS}
\def\Kb{{\sbm{K}}\XS}  \def\kb{{\sbm{k}}\XS}
\def\Lb{{\sbm{L}}\XS}  \def\lb{{\sbm{l}}\XS}
\def\Mb{{\sbm{M}}\XS}  \def\mb{{\sbm{m}}\XS}
\def\Nb{{\sbm{N}}\XS}  \def\nb{{\sbm{n}}\XS}
\def\Ob{{\sbm{O}}\XS}  \def\ob{{\sbm{o}}\XS}
\def\Pb{{\sbm{P}}\XS}  \def\pb{{\sbm{p}}\XS}
\def\Qb{{\sbm{Q}}\XS}  \def\qb{{\sbm{q}}\XS}
\def\Rb{{\sbm{R}}\XS}  \def\rb{{\sbm{r}}\XS}
\def\Sb{{\sbm{S}}\XS}  \def\sb{{\sbm{s}}\XS}
\def\Tb{{\sbm{T}}\XS}  \def\tb{{\sbm{t}}\XS}
\def\Ub{{\sbm{U}}\XS}  \def\ub{{\sbm{u}}\XS}
\def\Vb{{\sbm{V}}\XS}  \def\vb{{\sbm{v}}\XS}
\def\Wb{{\sbm{W}}\XS}  \def\wb{{\sbm{w}}\XS}
\def\Xb{{\sbm{X}}\XS}  \def\xb{{\sbm{x}}\XS}
\def\Yb{{\sbm{Y}}\XS}  \def\yb{{\sbm{y}}\XS}
\def\Zb{{\sbm{Z}}\XS}  \def\zb{{\sbm{z}}\XS}

% Compl�ment int�ressant...

\def\ellb{\sbmm{\ell}\XS}


% ALPHABET CURSIF MAJUSCULE, taille adaptative
% Cursif,               % Cursif gras

\def\Ac{{\scu{A}}\XS}   \def\Acb{{\scb{A}}\XS}
\def\Bc{{\scu{B}}\XS}   \def\Bcb{{\scb{B}}\XS}
\def\Cc{{\scu{C}}\XS}   \def\Ccb{{\scb{C}}\XS}
\def\Dc{{\scu{D}}\XS}   \def\Dcb{{\scb{D}}\XS}
\def\Ec{{\scu{E}}\XS}   \def\Ecb{{\scb{E}}\XS}
\def\Fc{{\scu{F}}\XS}   \def\Fcb{{\scb{F}}\XS}
\def\Gc{{\scu{G}}\XS}   \def\Gcb{{\scb{G}}\XS}
\def\Hc{{\scu{H}}\XS}   \def\Hcb{{\scb{H}}\XS}
\def\Ic{{\scu{I}}\XS}   \def\Icb{{\scb{I}}\XS}
\def\Jc{{\scu{J}}\XS}   \def\Jcb{{\scb{J}}\XS}
\def\Kc{{\scu{K}}\XS}   \def\Kcb{{\scb{K}}\XS}
\def\Lc{{\scu{L}}\XS}   \def\Lcb{{\scb{L}}\XS}
\def\Mc{{\scu{M}}\XS}   \def\Mcb{{\scb{M}}\XS}
\def\Nc{{\scu{N}}\XS}   \def\Ncb{{\scb{N}}\XS}
\def\Oc{{\scu{O}}\XS}   \def\Ocb{{\scb{O}}\XS}
\def\Pc{{\scu{P}}\XS}   \def\Pcb{{\scb{P}}\XS}
\def\Qc{{\scu{Q}}\XS}   \def\Qcb{{\scb{Q}}\XS}
\def\Rc{{\scu{R}}\XS}   \def\Rcb{{\scb{R}}\XS}
\def\Sc{{\scu{S}}\XS}   \def\Scb{{\scb{S}}\XS}
\def\Tc{{\scu{T}}\XS}   \def\Tcb{{\scb{T}}\XS}
\def\Uc{{\scu{U}}\XS}   \def\Ucb{{\scb{U}}\XS}
\def\Vc{{\scu{V}}\XS}   \def\Vcb{{\scb{V}}\XS}
\def\Wc{{\scu{W}}\XS}   \def\Wcb{{\scb{W}}\XS}
\def\Xc{{\scu{X}}\XS}   \def\Xcb{{\scb{X}}\XS}
\def\Yc{{\scu{Y}}\XS}   \def\Ycb{{\scb{Y}}\XS}
\def\Zc{{\scu{Z}}\XS}   \def\Zcb{{\scb{Z}}\XS}


% ALPHABET DROIT EN MATHS, taille adaptative

\def\AD{{\sdm{A}}\XS}  \def\aD{{\sdm{a}}\XS}
\def\BD{{\sdm{B}}\XS}  \def\bD{{\sdm{b}}\XS}
\def\CD{{\sdm{C}}\XS}  \def\cD{{\sdm{c}}\XS}
\def\DD{{\sdm{D}}\XS}  \def\dD{{\sdm{d}}\XS}
\def\ED{{\sdm{E}}\XS}  \def\eD{{\sdm{e}}\XS}
\def\FD{{\sdm{F}}\XS}  \def\fD{{\sdm{f}}\XS}
\def\GD{{\sdm{G}}\XS}  \def\gD{{\sdm{g}}\XS}
\def\HD{{\sdm{H}}\XS}  \def\hD{{\sdm{h}}\XS}
\def\ID{{\sdm{I}}\XS}  \def\iD{{\sdm{i}}\XS}
\def\JD{{\sdm{J}}\XS}  \def\jD{{\sdm{j}}\XS}
\def\KD{{\sdm{K}}\XS}  \def\kD{{\sdm{k}}\XS}
\def\LD{{\sdm{L}}\XS}  \def\lD{{\sdm{l}}\XS}
\def\MD{{\sdm{M}}\XS}  \def\mD{{\sdm{m}}\XS}
\def\ND{{\sdm{N}}\XS}  \def\nD{{\sdm{n}}\XS}
\def\OD{{\sdm{O}}\XS}  \def\oD{{\sdm{o}}\XS}
\def\PD{{\sdm{P}}\XS}  \def\pD{{\sdm{p}}\XS}
\def\QD{{\sdm{Q}}\XS}  \def\qD{{\sdm{q}}\XS}
\def\RD{{\sdm{R}}\XS}  \def\rD{{\sdm{r}}\XS}
\def\SD{{\sdm{S}}\XS}  \def\sD{{\sdm{s}}\XS}
\def\TD{{\sdm{T}}\XS}  \def\tD{{\sdm{t}}\XS}
\def\UD{{\sdm{U}}\XS}  \def\uD{{\sdm{u}}\XS}
\def\VD{{\sdm{V}}\XS}  \def\vD{{\sdm{D}}\XS}
\def\WD{{\sdm{W}}\XS}  \def\wD{{\sdm{w}}\XS}
\def\XD{{\sdm{X}}\XS}  \def\xD{{\sdm{x}}\XS}
\def\YD{{\sdm{Y}}\XS}  \def\yD{{\sdm{y}}\XS}
\def\ZD{{\sdm{Z}}\XS}  \def\zD{{\sdm{z}}\XS}


% ALPHABET GRAS DROIT, taille adaptative

\def\Av{{\sbv{A}}\XS}  \def\av{{\sbv{a}}\XS}
\def\Bv{{\sbv{B}}\XS}  \def\bv{{\sbv{b}}\XS}
\def\Cv{{\sbv{C}}\XS}  \def\cv{{\sbv{c}}\XS}
\def\Dv{{\sbv{D}}\XS}  \def\dv{{\sbv{d}}\XS}
\def\Ev{{\sbv{E}}\XS}  \def\ev{{\sbv{e}}\XS}
\def\Fv{{\sbv{F}}\XS}  \def\fv{{\sbv{f}}\XS}
\def\Gv{{\sbv{G}}\XS}  \def\gv{{\sbv{g}}\XS}
\def\Hv{{\sbv{H}}\XS}  \def\hv{{\sbv{h}}\XS}
\def\Iv{{\sbv{I}}\XS}  \def\iv{{\sbv{i}}\XS}
\def\Jv{{\sbv{J}}\XS}  \def\jv{{\sbv{j}}\XS}
\def\Kv{{\sbv{K}}\XS}  \def\kv{{\sbv{k}}\XS}
\def\Lv{{\sbv{L}}\XS}  \def\lv{{\sbv{l}}\XS}
\def\Mv{{\sbv{M}}\XS}  \def\mv{{\sbv{m}}\XS}
\def\Nv{{\sbv{N}}\XS}  \def\nv{{\sbv{n}}\XS}
\def\Ov{{\sbv{O}}\XS}  \def\ov{{\sbv{o}}\XS}
\def\Pv{{\sbv{P}}\XS}  \def\pv{{\sbv{p}}\XS}
\def\Qv{{\sbv{Q}}\XS}  \def\qv{{\sbv{q}}\XS}
\def\Rv{{\sbv{R}}\XS}  \def\rv{{\sbv{r}}\XS}
\def\Sv{{\sbv{S}}\XS}  \def\sv{{\sbv{s}}\XS}
\def\Tv{{\sbv{T}}\XS}  \def\tv{{\sbv{t}}\XS}
\def\Uv{{\sbv{U}}\XS}  \def\uv{{\sbv{u}}\XS}
\def\Vv{{\sbv{V}}\XS}  \def\vv{{\sbv{v}}\XS}
\def\Wv{{\sbv{W}}\XS}  \def\wv{{\sbv{w}}\XS}
\def\Xv{{\sbv{X}}\XS}  \def\xv{{\sbv{x}}\XS}
\def\Yv{{\sbv{Y}}\XS}  \def\yv{{\sbv{y}}\XS}
\def\Zv{{\sbv{Z}}\XS}  \def\zv{{\sbv{z}}\XS}


% Nouveau (08/99) : ALPHABET BLACKBOARD, taille adaptative

\def\Abb{{\sbl{A}}\XS}  \def\abb{{\sbl{a}}\XS}
\def\Bbb{{\sbl{B}}\XS}  \def\bbb{{\sbl{b}}\XS}
\def\Cbb{{\sbl{C}}\XS}  \def\cbb{{\sbl{c}}\XS}
\def\Dbb{{\sbl{D}}\XS}  \def\dbb{{\sbl{d}}\XS}
\def\Ebb{{\sbl{E}}\XS}  \def\ebb{{\sbl{e}}\XS}
\def\Fbb{{\sbl{F}}\XS}  \def\fbb{{\sbl{f}}\XS}
\def\Gbb{{\sbl{G}}\XS}  \def\gbb{{\sbl{g}}\XS}
\def\Hbb{{\sbl{H}}\XS}  \def\hbb{{\sbl{h}}\XS}
\def\Ibb{{\sbl{I}}\XS}  \def\ibb{{\sbl{i}}\XS}
\def\Jbb{{\sbl{J}}\XS}  \def\jbb{{\sbl{j}}\XS}
\def\Kbb{{\sbl{K}}\XS}  \def\kbb{{\sbl{k}}\XS}
\def\Lbb{{\sbl{L}}\XS}  \def\lbb{{\sbl{l}}\XS}
\def\Mbb{{\sbl{M}}\XS}  \def\mbb{{\sbl{m}}\XS}
\def\Nbb{{\sbl{N}}\XS}  \def\nbb{{\sbl{n}}\XS}
\def\Obb{{\sbl{O}}\XS}  \def\obb{{\sbl{o}}\XS}
\def\Pbb{{\sbl{P}}\XS}  \def\pbb{{\sbl{p}}\XS}
\def\Qbb{{\sbl{Q}}\XS}  \def\qbb{{\sbl{q}}\XS}
\def\Rbb{{\sbl{R}}\XS}  \def\rbb{{\sbl{r}}\XS}
\def\Sbb{{\sbl{S}}\XS}  \def\sbb{{\sbl{s}}\XS}
\def\Tbb{{\sbl{T}}\XS}  \def\tbb{{\sbl{t}}\XS}
\def\Ubb{{\sbl{U}}\XS}  \def\ubb{{\sbl{u}}\XS}
\def\Vbb{{\sbl{V}}\XS}  \def\vbb{{\sbl{v}}\XS}
\def\Wbb{{\sbl{W}}\XS}  \def\wbb{{\sbl{w}}\XS}
\def\Xbb{{\sbl{X}}\XS}  \def\xbb{{\sbl{x}}\XS}
\def\Ybb{{\sbl{Y}}\XS}  \def\ybb{{\sbl{y}}\XS}
\def\Zbb{{\sbl{Z}}\XS}  \def\zbb{{\sbl{z}}\XS}


% ALPHABET GREC GRAS italique, taille adaptative

\def\alphab      {{\sbmm{\alpha}}\XS}
\def\betab       {{\sbmm{\beta}}\XS}
\def\gammab      {{\sbmm{\gamma}}\XS}      \def\Gammab    {{\sbmm{\Gamma}}\XS}
\def\deltab      {{\sbmm{\delta}}\XS}      \def\Deltab    {{\sbmm{\Delta}}\XS}
\def\epsilonb    {{\sbmm{\epsilon}}\XS}
\def\varepsilonb {{\sbmm{\varepsilon}}\XS}
\def\etab        {{\sbmm{\eta}}\XS}
\def\thetab      {{\sbmm{\theta}}\XS}      \def\Thetab    {{\sbmm{\Theta}}\XS}
\def\varthetab   {{\sbmm{\vartheta}}\XS}
\def\iotab       {{\sbmm{\iota}}\XS}
\def\kappab      {{\sbmm{\kappa}}\XS}
\def\varkappab   {{\sbmm{\varkappa}}\XS}
\def\lambdab     {{\sbmm{\lambda}}\XS}     \def\Lambdab   {{\sbmm{\Lambda}}\XS}
\def\mub         {{\sbmm{\mu}}\XS}
\def\nub         {{\sbmm{\nu}}\XS}
\def\xib         {{\sbmm{\xi}}\XS}         \def\Xib       {{\sbmm{\Xi}}\XS}
\def\pib         {{\sbmm{\pi}}\XS}         \def\Pib       {{\sbmm{\Pi}}\XS}
\def\varpib      {{\sbmm{\varpi}}\XS}
\def\rhob        {{\sbmm{\rho}}\XS}
\def\varrhob     {{\sbmm{\varrho}}\XS}
\def\sigmab      {{\sbmm{\sigma}}\XS}      \def\Sigmab    {{\sbmm{\Sigma}}\XS}
\def\varsigmab   {{\sbmm{\varsigma}}\XS}   \def\Varsigmab {{\sbmm{\Varsigma}}\XS}
\def\phib        {{\sbmm{\phi}}\XS}        \def\Phib      {{\sbmm{\Phi}}\XS}
\def\varphib     {{\sbmm{\varphi}}\XS}
\def\chib        {{\sbmm{\chi}}\XS}
\def\psib        {{\sbmm{\psi}}\XS}        \def\Psib      {{\sbmm{\Psi}}\XS}
\def\omegab      {{\sbmm{\omega}}\XS}      \def\Omegab    {{\sbmm{\Omega}}\XS}
\def\taub        {{\sbmm{\tau}}\XS}
\def\upsilonb    {{\sbmm{\upsilon}}\XS}    \def\Upsilonb  {{\sbmm{\Upsilon}}\XS}
\def\zetab       {{\sbmm{\zeta}}\XS}


% CHIFFRES
% Gras et Blackboard, taille adaptative

\def\zerob   {{\sbv{0}}\XS}        \def\zerobb   {{\sbl{O}}\XS}
\def\unb     {{\sbv{1}}\XS}        \def\unbb     {{\sbl{1}}\XS}
\def\deuxb   {{\sbv{2}}\XS}        \def\deuxbb   {{\sbl{2}}\XS}
\def\troisb  {{\sbv{3}}\XS}        \def\troisbb  {{\sbl{3}}\XS}
\def\quatreb {{\sbv{4}}\XS}        \def\quatrebb {{\sbl{4}}\XS}
\def\cinqb   {{\sbv{5}}\XS}        \def\cinqbb   {{\sbl{5}}\XS}
\def\sixb    {{\sbv{6}}\XS}        \def\sixbb    {{\sbl{6}}\XS}
\def\septb   {{\sbv{7}}\XS}        \def\septbb   {{\sbl{7}}\XS}
\def\huitb   {{\sbv{8}}\XS}        \def\huitbb   {{\sbl{8}}\XS}
\def\neufb   {{\sbv{9}}\XS}        \def\neufbb   {{\sbl{9}}\XS}


% QUELQUES ENSEMBLES

\def\eC{\Cbb}
\def\eE{\Ebb}
\def\eN{\Nbb}
\def\eQ{\Qbb}
\def\eR{\Rbb}
\def\eZ{\Zbb}


beginend.tex
% TEX 7(ascii) bits
% BEGINEND.tex                          LaTeX document
% Author: J.I.          Date : 07/96    Derniere modif le 10/07/96
% Raccourcis des \begin et \end, d'apr\`es les Grands Anciens
% Modification provisoire le 12/02/97 sur \beq#1\eeq --->
% ---> Compatibilit\'e avec Amsmath.sty

\def\barr{\begin{array}}                \def\earr{\end{array}}
\def\bcc{\begin{center}}                \def\ecc{\end{center}}
\def\cl#1{\centerline{#1}}
\def\bdoc{\begin{document}}             \def\edoc{\end{document}}
\def\ben{\begin{enumerate}}             \def\een{\end{enumerate}}
\def\beqn{\begin{eqnarray}}             \def\eeqn{\end{eqnarray}}
\def\beqnl#1{\beqn\label{#1}}           \def\eeqnl#1{\label{#1}\eeqn}
\def\beqnx{\begin{eqnarray*}}           \def\eeqnx{\end{eqnarray*}}
\def\bseqn{\begin{subeqnarray}}         \def\eseqn{\end{subeqnarray}}
%
\def\beq#1\eeq{\begin{equation}#1\end{equation}}
%\def\beq{\begin{equation}}             \def\eeq{\end{equation}}
%\def\beql#1{\beq\label{#1}}            \def\eeql#1{\label{#1}\eeq}
%
\def\bal#1\eal{\begin{align}#1\end{align}}
\def\balx#1\ealx{\begin{align*}#1\end{align*}}
%
\def\beqx{$$}                           \def\eeqx{$$}
\def\bfig{\protect\begin{figure}}       \def\efig{\protect\end{figure}}
\def\bfigx{\protect\begin{figure*}}     \def\efigx{\protect\end{figure*}}
\def\bfigt{\protect\begin{figurette}}   \def\efigt{\protect\end{figurette}}
\def\bfl{\begin{flushleft}}             \def\efl{\end{flushleft}}
\def\bfr{\begin{flushright}}            \def\efr{\end{flushright}}
\def\bit{\begin{itemize}}               \def\eit{\end{itemize}}
\def\bmi{\begin{minipage}}              \def\emi{\end{minipage}}
%\def\bfmi{\begin{fminipage}}            \def\efmi{\end{fminipage}}
\newcommand{\bfmi}{\begin{fminipage}}   \newcommand{\efmi}{\end{fminipage}}
\def\bpic{\begin{picture}}              \def\epic{\end{picture}}
\def\btabl{\begin{table}}               \def\etabl{\end{table}}
\def\btab{\begin{tabular}} %\def\etab{\end{tabular}} Conflit avec eta gras... Elle est pas grasse, ma soeur !
\def\btabu{\begin{tabular}}             \def\etabu{\end{tabular}}
\def\btabx{\begin{tabular*}}            \def\etabx{\end{tabular*}}

jidef.tex

%\tiny\scriptsize\footnotesize\small\normalsize\large\Large\LARGE\huge\Huge

%\setlength{\abovedisplayskip}{6pt}         % Pour bricoler l'espace,
%\setlength{\belowdisplayskip}{6pt}         % Autour des equations...
%\setlength{\abovedisplayshortskip}{6pt}    % (Bon d'accord, je finasse)
%\setlength{\belowdisplayshortskip}{6pt}    % Pas d'improvisation!

% pour tasser la biblio: ajouter \setlength{\itemsep}{0pt} en deuxieme ligne du .bbl

\RequirePackage{accents} % e.g. \frown
\RequirePackage{amssymb} % e.g. \geqslant

\def\nnb{\nonumber}
\def\dd{\,\dD}                         % doit etre droit en francais
\def\leq{\leqslant}
\def\geq{\geqslant}
\renewcommand{\det}[1]{\bars{#1}}
\def\V{,\kern-.1em}  	% Virgule pour les nombres decimaux
%\def\V{.}              % Point pour les anglo-saxons
\def\e#1{.10^{#1}}      % Notation scientifique a la francaise
%\def\e#1{\rm{e}^{#1}}  % Notation scientifique anglo-saxonne

%\def\et{{\normalfont et }}
\def\et{{et}\ }
\def\numero{}
\def\editornames{(dir.)}
\def\editorname{(dir.)}
\def\eqname{�quation~}
\def\eqnames{�quations~}
\def\tabname{tableau~}
\def\tabnames{tableaux~}
\def\figname{figure~}
\def\fignames{figures~}
\def\chapname{chapitre~}
\def\chapnames{chapitres~}
\def\parname{pa\-ra\-gra\-phe~}
\def\parnames{pa\-ra\-gra\-phes~}
\def\sectname{section~}
\def\theoname{th�or�me~}
\def\gui#1{\<#1\>}
%\def\<{<<\,}\def\>{\,>>}
\def\<{\og}\def\>{\fg}
\def\ier{-ier\XS}
\def\iere{-i�re\XS}
\def\ieme{-i�me\XS}

\def\para{\cl{\btabu{c}$\ast$~~~$\ast$\\$\ast$\\[-.2cm]\etabu}}
\def\Separe{\medskip\cl{\rule{3cm}{.02cm}}}
\def\Fin{\cl{\rule{1.5cm}{.02cm}}\cl{\rule{.75cm}{.02cm}}}  % deux traits pour fin

%------------------------------------------------------------------
% Mise en commentaire rapide !
\def\rem#1{}                    % rem{bla bla bla}

\def\In{in }
\def\volume{vol.\,}
\def\PAGE{p.\,}
\def\PAGES{p.\,}
%\def\editionname{�dition}

\newcommand{\tabref}[2][]{\tabname\ref{#2}{\normalfont#1}\XS}
\newcommand{\figref}[2][]{\figname\ref{#2}{\normalfont#1}\XS}
\newcommand{\refeq}[1]{\eqname(\ref{#1})\XS}
\newcommand{\refeqs}[1]{\eqnames(\ref{#1})\XS}
\def\chapref#1{\chapname\ref{#1}\XS}
\def\parref#1{\parname\ref{#1}\XS}
\def\parsref#1{\parnames\ref{#1}\XS}
\def\sectref#1{\sectname\ref{#1}\XS}
\def\refpage#1{\page\pageref{#1}\XS}
\def\reff{\XS\eqref}

\newcommand{\page}{p.~}
\newcommand{\pages}{pp.~}
%\def\nompage{p.\,}
%\def\nompages{pp.\,}
\def\numero{n\ensuremath{^{\circ}}}
\def\primo{\ensuremath{1^{\circ}}\XS}
\def\secundo{\ensuremath{2^{\circ}}\XS}
\def\tertio{\ensuremath{3^{\circ}}\XS}
\def\quatro{\ensuremath{4^{\circ}}\XS}

\def\B{\text{\large$\bullet$}}
\def\incirc#1{\pscirclebox[framesep=1pt]{\scriptsize#1}}
\def\incircp#1{\pscirclebox[framesep=1pt]{\tiny#1}}
\def\outcirc#1{\raisebox{-3pt}{\huge$\times$}\hspace*{-.45cm}\text{\incirc{#1}}}

\def\traitplein{\text{$\kern-1pt\frac{\kern5pt}{\kern7pt}\kern-2pt$}}
\def\tirets{\text{\kern-1pt-{}-{}-\kern-1pt}}
\def\pointilles{\text{$\cdot$$\cdot$$\cdot$\kern-1pt}}

\def\para{~\\[-.75cm]\bcc $\ast$~~~$\ast$\\$\ast$\ecc~\\[-.75cm]}

%-- Accolades, parenth�ses, etc --------------------------------------
        
\def\pth#1{\left(#1\right)}                \def\stdpth#1{(#1)}
\def\acc#1{\left\{#1\right\}}              \def\stdacc#1{\{#1\}}
\def\cro#1{\left[#1\right]}                \def\stdcro#1{[#1]}
\def\croext#1{\left]#1\right[}             \def\stdcroext#1{]#1[}
\def\bars#1{\left|#1\right|}               \def\stdbars#1{|#1|}
\def\norm#1{\left\|#1\right\|}             \def\stdnorm#1{\|#1\|}
\def\scal#1{\left\langle#1\right\rangle}   \def\stdscal#1{\langle#1\rangle}
 
\def\bigpth#1{\bigl(#1\bigr)}              \def\biggpth#1{\biggl(#1\biggr)}
\def\bigacc#1{\bigl\{#1\bigr\}}            \def\biggacc#1{\biggl\{#1\biggr\}}
\def\bigcro#1{\bigl[#1\bigr]}              \def\biggcro#1{\biggl[#1\biggr]}
\def\bigcroext#1{\bigl]#1\bigr[}           \def\biggcroext#1{\biggl]#1\biggr[}
\def\bigbars#1{\bigl|#1\bigr|}             \def\biggbars#1{\biggl|#1\biggr|}
\def\bignorm#1{\bigl\|#1\bigr\|}           \def\biggnorm#1{\biggl\|#1\biggr\|}
\def\bigscal#1{\bigl\langle#1\bigr\rangle} \def\biggscal#1{\biggl\langle#1\biggr\rangle}

\def\Bigpth#1{\Bigl(#1\Bigr)}              \def\Biggpth#1{\Biggl(#1\Biggr)}
\def\Bigacc#1{\Bigl\{#1\Bigr\}}            \def\Biggacc#1{\Biggl\{#1\Biggr\}}
\def\Bigcro#1{\Bigl[#1\Bigr]}              \def\Biggcro#1{\Biggl[#1\Biggr]}
\def\Bigcroext#1{\Bigl]#1\Bigr[}           \def\Biggcroext#1{\Biggl]#1\Biggr[}
\def\Bigbars#1{\Bigl|#1\Bigr|}             \def\Biggbars#1{\Biggl|#1\Biggr|}
\def\Bignorm#1{\Bigl\|#1\Bigr\|}           \def\Biggnorm#1{\Biggl\|#1\Biggr\|}
\def\Bigscal#1{\Bigl\langle#1\Bigr\rangle} \def\Biggscal#1{\Biggl\langle#1\Biggr\rangle}

%------------------------------------------------------------------
\def\cl#1{\centerline{#1}} \def\eqcl#1{\cl{$\displaystyle#1$}} \def\eq#1{{$\displaystyle#1$}}

\def\nablab{\text{\boldmath$\nabla$}}
\def\wh#1{\widehat{#1}}    \def\wt#1{\widetilde{#1}}
\def\foo#1{{\text{\def\scad{-3}\tiny #1}}}
\def\T{^\tD}
%\def\T{^T}
\def\M{^{-1}}
\def\MT{^{-\tD}}
\def\+{^\dagger}
\def\I{\,|\,}           % "sachant" bien espac� pour les formules
\def\J{\mathop{\,;\,}}  % "point virgule" bien espac� pour les formules.
\def\w{,\thinspace}
\def\W{,\thickspace}
\def\ldotsv{,\,\ldots,\,}
\def\Exp#1{\exp\acc{#1}}
\def\Expp#1{e^{#1}}
\def\argmax{\mathop{\mathrm{arg\,max}}} % Mieux que \def\argmax{\arg\max}
\def\argmin{\mathop{\mathrm{arg\,min}}} % car l'indice est reparti
\def\Vect{\mathop{\text{Vect}}}
\def\Card{\mathrm{Card}}
\def\Circul#1{\mathrm{Circul}\{#1\}}
\def\toeplitz#1{\mathrm{T{\oe}plitz}\{#1\}}
\def\Id{\mathit{Id}}
\def\cste{\mathrm{cste}}
\def\trace{\mathrm{trace}}
\def\rang{{\mathrm{rang}}}
\def\Cov{\mathrm{Cov}}
\def\diag{{\mathrm{diag}}}
\def\Diag{\mathrm{Diag}}
\def\DIAG#1{\mathrm{Diag}\acc{#1}}
\def\Ker{\mathrm{Ker\,}}
\def\Img{\mathrm{Im\,}}
\def\Div{{\mathrm{div}}}
\def\FT#1{\protect\accentset{\frown}{#1}}
% Pour mettre un \No{} au-dessus d'une lettre (\TM, 02/97) (package{amsmath.sty})
%\def\rond#1{\overset{\kern-0.33em~_\circ}{#1}}
%\def\rondit[#1]#2{\overset{\kern#1~_\circ}{#2}}
\def\EL{\mathrm{EL}}
\def\esp{\ED\,}
\newcommand{\Esp}[2][]{\ED{_{#1}}\kern-2pt\cro{#2}}
%\def\Esp#1{\ED\kern-1pt\pth{#1}}
\def\stdEsp#1{\ED\stdcro{#1}}
%\def\EMQ{\mathit{EMQ}}
\def\var{\mathrm{var}\,}	\def\Var#1{\mathrm{var}\pth{#1}}	
\def\cov{\mathrm{cov}}		\def\Cov#1{\mathrm{cov}\kern-2pt\pth{#1}}
\def\egal{\kern-.5em=\kern-.5em}
\def\egdef{\stackrel{\Delta}{=}}
\newcommand{\imag}[1]{\mathrm{Im}\acc{#1}}
\newcommand{\reel}[1]{\mathrm{Re}\acc{#1}}
\newcommand{\impur}{\ensuremath{j}}

\def\sump{\mathop{\raisebox{-.3ex}{\text{\Large$\Sigma$}}}}
\def\intp{{\text{\small$\int$}}}
\def\prodp{\mathop{\raisebox{-.3ex}{\text{\Large$\Pi$}}}}
\newcommand{\fracp}[2]{\frac{\raisebox{-1pt}{$#1$}}{#2}}
\newcommand{\UNsur}[1]{\fracp{1}{#1}}
\def\demi{\UNsur{2}}
\def\ds{\displaystyle}
\def\ss{\scriptstyle}
\def\fracds#1#2{{\ds\frac{#1}{#2}}}
\def\fracts#1#2{{\ss\frac{#1}{#2}}}
\def\derpar#1#2{{\frac{\partial #1}{\partial #2}}}
\def\froc#1#2{{#1/#2}}     \def\fric#1#2{\frac1{#2}#1}
\def\derpor#1#2{{\froc{\partial #1}{\partial #2}}}
\def\parsec#1#2#3{{\frac{\partial^2 #1}{\partial #2\,\partial #3}}}
\def\parsecd#1#2{{\frac{\partial^2 #1}{{\partial #2}^2}}}
\def\porsec#1#2#3{{\froc{\partial^2 #1}{\partial #2\,\partial #3}}}
\def\porsecd#1#2{{\froc{\partial^2 #1}{{\partial #2}^2}}}
\def\pornd#1#2#3{{\froc{\partial^{#3} #1}{{\partial #2}^{#3}}}}
%\def\sur{\ensuremath{\kern-1pt/\kern-1pt}}
\def\diff#1#2{{\frac{\dd#1}{\dd#2}}}

%------------------------------------------------------------

\newsavebox{\fminibox}
\newlength{\fminilength}
\newenvironment{fminipage}[1][\linewidth]
  {\setlength{\fminilength}{#1}%-2\fboxsep-2\fboxrule%
   \begin{lrbox}{\fminibox}\begin{minipage}{\fminilength}}
  {\end{minipage}\end{lrbox}\noindent\fbox{\usebox{\fminibox}}}

\def\arrayp{\renewcommand{\arraystretch}{.7}\setlength{\arraycolsep}{2pt}}
\def\tabp{\renewcommand{\arraystretch}{.7}\setlength{\tabcolsep}{2pt}}

% Centre verticalement la colonne :
\def\cc#1{\setlength{\tabcolsep}{0pt}\btabu{c}#1\etabu}
\newcommand{\figc}[2][]
   {\setlength{\tabcolsep}{0pt}\btabu{c}\includegraphics[#1]{#2}\etabu}
   %{\setlength{\tabcolsep}{0pt}\btabu{c}\includegraphics[draft=true,#1]{#2}\etabu}
   %{\btabu{c}\small\texttt{\btabu{c}includegraphics[#1]\\\{\url{#2}\}\etabu}\etabu}

%\def\yaxis#1{\cc{\rotateleft{#1}}}
\def\ylabel#1{\cc{\rotateleft{#1}}}

\newlength{\vshift} \newlength{\hshift} % Pour hspace et vspace en pstricks

%\boite{length}{blabla} utilisable dans \caption... marche pas !
%\def\boite#1#2{#2\settowidth{\hshift}{#2}\addtolength{\hshift}{-#1}\hspace*{-\hshift}}
%\boite{blabla1}{blabla2} utilisable dans \caption
\def\boite#1#2{#2\settowidth{\hshift}{#1}\hspace*{\hshift}\settowidth{\hshift}{#2}\hspace*{-\hshift}}

%------------------------------------------------------------------

\def\UP#1{\uppercase{#1}}

\RequirePackage{ifthen}
\RequirePackage{calc}

\newcommand{\taille}[1][\scad]{%
\ifthenelse{#1 = -5}{}{}%
\ifthenelse{#1 = -4}{\tiny}{}%
\ifthenelse{#1 = -3}{\scriptsize}{}%
\ifthenelse{#1 = -2}{\footnotesize}{}%
\ifthenelse{#1 = -1}{\small}{}%
\ifthenelse{#1 = 0}{\normalsize}{}%
\ifthenelse{#1 = 1}{\large}{}%
\ifthenelse{#1 = 2}{\Large}{}%
\ifthenelse{#1 = 3}{\LARGE}{}%
\ifthenelse{#1 = 4}{\huge}{}%
\ifthenelse{#1 = 5}{\Huge}{}}

\newboolean{@serif}
\setboolean{@serif}{true}
\def\scad{-5} % scadefaultsize

\newcounter{taille}
\newcommand{\sca}[2][\scad]{\setcounter{taille}{#1}%
  \ifthenelse{\boolean{@serif}}
  {{\taille[\thetaille]\textsc{#2}}}
  {\setcounter{taille}{\value{taille}-1}{\uppercase{\taille[\thetaille]#2}}}}

%------------------------------------------------------------

\def\xbwh{\wh{\xb}{}}
%\def\qbwh{\kern1pt\wh{\kern-1pt\qb}{}}
\def\Xbwh{\kern1pt\wh{\kern-1pt\Xb}}
\def\Xvwh{\wh{\Xv}{}}
\def\Xvh#1{\Xvwh{}^{\foo{#1}}}
\def\xbh#1{\xbwh^{\foo{\normalfont#1}}}
\def\xh#1{\xwh^{\foo{\normalfont#1}}}
\def\xwh{\wh{x}{}}
\def\whtfx{{\stackrel{\raisebox{-5pt}{$\scriptscriptstyle\smash{\frown}$}}{\smash{\kern-.4pt \xwh}}}{}}
%\def\bvwh{\wh{\bv}{}}
%\def\hvwh{\wh{\hv}{}}
\def\bbwh{\wh{\bb}{}}
\def\hbwh{\wh{\hb}{}}
\def\thetabh#1{\wh{\thetab}{}^{\foo{\normalfont#1}}}

%------------------------------------------------------------------

\newcommand\unit[1]{\,\text{#1}\XS}
\def\Mo{\unit{Mo}\XS}
\def\MHz{\unit{MHz}\XS}
\def\kHz{\unit{kHz}\XS}
\def\Hz{\unit{Hz}\XS}
\def\cm{\unit{cm}\XS}
\def\mm{\unit{{mm}}\XS}
\def\ms{\unit{ms}\XS}
\def\dB{\unit{dB}\XS}
\def\sec{\unit{s}\XS}
\def\arcsec{\unit{arcsec}\XS}

\def\const{\textrm{const.}}
\def\SI{\text{si\:}}       \def\Si{\text{Si\:}}
\def\ET{\text{et\:}}       \def\OU{\text{ou\:}}
\def\ALORS{\text{alors\:}}
\def\DOU{\text{d'o\`u\:}}  \def\Ou{\text{o\`u\:}}
\def\QUAND{\text{quand\:}}
\def\POUR{\text{pour\:}}   \def\POURTOUT{\text{pour tout\:}}
\def\SC{\text{s.\,c.\:}}   \def\SOUSC{\text{sous contraintes\:}}
\def\SINON{\text{sinon}}
\def\AVEC{\text{avec\:}}
\def\DANS{\text{dans\:}}
\def\DONC{\text{donc}}

\def\adhoc{\textit{ad hoc}\XS}
\def\apost{\textit{a posteriori}\XS}
\def\Apost{\textit{A posteriori}\XS}
\def\APost{\textit{A Posteriori}\XS}
\def\post{\textit{posterior}\XS}
\def\Post{\textit{Posterior}\XS}
\def\afortio{\textit{a fortiori}\XS}
\def\Afortio{\textit{A fortiori}\XS}
\def\aprio{\textit{a priori}\XS}
\def\Aprio{\textit{A priori}\XS}
\def\Cad{C'est-\`a-dire\XS}
\def\cad{c'est-\`a-dire\XS}
\def\coeur{c{\oe}ur\XS}
\def\cf{\textit{cf.}\XS}
\def\cv{\textit{curriculum vit\ae}\XS}
\def\Cv{\textit{Curriculum vit\ae}\XS}
\def\CV{\textit{Curriculum Vit\ae}\XS}
\def\eg{\textit{e.g.,}\XS}
\def\andname{et\xspace}
\def\etal{{\itshape et coll.}}
\def\etc{{etc.}}
\def\etcoll{{\itshape et coll.}}
\def\idem{\textit{idem}\XS}
\def\ie{\textit{i.e.},\XS}
\def\IFF{\textit{if and only if}\XS}
\def\iid{i.i.d.\XS}
\def\infine{\textit{in fine}\XS}
\def\insit{\textit{in situ}\XS}
\def\inviv{\textit{in vivo}\XS}
\def\oeuvre{{\oe}uvre\XS}
\def\oeuvres{{\oe}uvres\XS}
\def\recto{\textit{recto}\XS}
\def\Toeplitz{T{\oe}plitz\XS}
\def\vav{vis-\`a-vis\XS}
\def\via{\textit{via}\XS}
\def\vic{\textit{vice-versa}\XS}
\def\verso{\textit{verso}\XS}
%\def\vitae{\textit{vit\ae}\XS}
\def\vs{\textit{vs}\XS}
\def\versus{\textit{versus}\XS}
\def\ssi{\textit{si et seulement si}\XS}
\def\strsen{\textit{stricto sensu}\XS}

%------------------------------------------------------------

\newcommand{\aic}[1][\scad]{\sca[#1]{aic}\XS}
\newcommand{\ar}[1][\scad]{\sca[#1]{ar}\XS}
\newcommand{\art}[1][\scad]{\sca[#1]{art}\XS}
\newcommand{\arma}[1][\scad]{\sca[#1]{arma}\XS}
\newcommand{\ma}[1][\scad]{\sca[#1]{ma}\XS}
%\newcommand{\bttb}[1][\scad]{\sca[#1]{bttb}\XS}
\newcommand{\bccb}[1][\scad]{\sca[#1]{bccb}\XS}
\newcommand{\bg}[1][\scad]{\sca[#1]{bg}\XS}
\newcommand{\dbg}[1][\scad]{\sca[#1]{dbg}\XS}
\newcommand{\cat}[1][\scad]{\sca[#1]{cat}\XS}
\newcommand{\ccd}[1][\scad]{\sca[#1]{ccd}\XS}
\newcommand{\cgm}[1][\scad]{\sca[#1]{cgm}\XS}
\newcommand{\cgmgg}[1][\scad]{\sca[#1]{cgmgg}\XS}
\newcommand{\clean}[1][\scad]{\sca[#1]{clean}\XS}
\newcommand{\clv}[1][\scad]{\sca[#1]{clv}\XS}
\newcommand{\clvh}[1][\scad]{\sca[#1]{clvh}}
\newcommand{\clvm}[1][\scad]{\sca[#1]{clvm}\XS}
\newcommand{\clpv}[1][\scad]{\sca[#1]{clpv}\XS}
\newcommand{\cnd}[1][\scad]{\sca[#1]{cnd}\XS}
\newcommand{\ct}[1][\scad]{\sca[#1]{ct}\XS}
\newcommand{\CP}[1][\scad]{\sca[#1]{cp}\XS}
\newcommand{\dse}[1][\scad]{\sca[#1]{dse}\XS}
\newcommand{\dsp}[1][\scad]{\sca[#1]{dsp}\XS}
\newcommand{\edp}[1][\scad]{\sca[#1]{edp}\XS}
\newcommand{\EM}[1][\scad]{\sca[#1]{em}\XS}
\newcommand{\EQM}{\sca{eqm}\XS}
\newcommand{\sem}[1][\scad]{\sca[#1]{sem}\XS}
\newcommand{\sirt}[1][\scad]{\sca[#1]{sirt}\XS}
\newcommand{\mvg}[1][\scad]{\sca[#1]{mvg}\XS}
\newcommand{\fft}[1][\scad]{\sca[#1]{fft}\XS}
\newcommand{\fpe}[1][\scad]{\sca[#1]{fpe}\XS}
\newcommand{\fk}[1][\scad]{\sca[#1]{fk}\XS}
\newcommand{\gc}[1][\scad]{\sca[#1]{gc}\XS}
\newcommand{\gcp}[1][\scad]{\sca[#1]{gcp}\XS}
\newcommand{\gcpc}[1][\scad]{\sca[#1]{gcpc}\XS}
\newcommand{\gnc}[1][\scad]{\sca[#1]{gnc}\XS}
\newcommand{\icd}[1][\scad]{\sca[#1]{icd}\XS}
\newcommand{\ig}[1][\scad]{\sca[#1]{ig}\XS}
\newcommand{\icm}[1][\scad]{\sca[#1]{icm}\XS}
\newcommand{\iwm}[1][\scad]{\sca[#1]{iwm}\XS}
\newcommand{\irm}[1][\scad]{\sca[#1]{irm}\XS}
\newcommand{\ldhy}[1][\scad]{\sca[#1]{l2h}y\XS}
\newcommand{\dldhy}[1][\scad]{\sca[#1]{dl2h}y\XS}
\newcommand{\ldlp}[1][\scad]{\sca[#1]{l2lp}\XS}
\newcommand{\lk}[1][\scad]{\sca[#1]{lk}\XS}
\newcommand{\map}[1][\scad]{\sca[#1]{map}\XS}
\newcommand{\mc}[1][\scad]{\sca[#1]{mc}\XS}
\newcommand{\mca}[1][\scad]{\sca[#1]{mca}\XS}
\newcommand{\mcmc}[1][\scad]{\sca[#1]{mcmc}\XS}
\newcommand{\mcreg}[1][\scad]{\sca[#1]{mcr}eg\XS}
\newcommand{\mdl}[1][\scad]{\sca[#1]{mdl}\XS}
\newcommand{\mem}[1][\scad]{\sca[#1]{mem}\XS}
\newcommand{\mmap}[1][\scad]{\sca[#1]{mmap}\XS}
\newcommand{\mpm}[1][\scad]{\sca[#1]{mpm}\XS}
\newcommand{\MP}[1][\scad]{\sca[#1]{mp}\XS}
\newcommand{\mV}[1][\scad]{\sca[#1]{mv}\XS}
\newcommand{\os}[1][\scad]{\sca[#1]{os}\XS}
\newcommand{\pC}[1][\scad]{\sca[#1]{pc}\XS}
\newcommand{\pocs}[1][\scad]{\sca[#1]{pocs}\XS}
\newcommand{\psf}[1][\scad]{\sca[#1]{psf}\XS}
\newcommand{\ri}[1][\scad]{\sca[#1]{ri}\XS}
\newcommand{\rif}[1][\scad]{\sca[#1]{rif}\XS}
\newcommand{\rp}[1][\scad]{\sca[#1]{rp}\XS}
\newcommand{\rls}[1][\scad]{\sca[#1]{rls}\XS}
\newcommand{\rsb}[1][\scad]{\sca[#1]{rsb}\XS}
\newcommand{\smlr}[1][\scad]{\sca[#1]{smlr}\XS}
\newcommand{\mmlr}[1][\scad]{\sca[#1]{mmlr}\XS}
\newcommand{\sage}[1][\scad]{\sca[#1]{sage}\XS}
\newcommand{\sar}[1][\scad]{\sca[#1]{sar}\XS}
\newcommand{\spect}[1][\scad]{\sca[#1]{spect}\XS}
\newcommand{\SF}[1][\scad]{\sca[#1]{sf}\XS}
\newcommand{\sq}[1][\scad]{\sca[#1]{sq}\XS}
\newcommand{\svd}[1][\scad]{\sca[#1]{svd}\XS}
\newcommand{\tcd}[1][\scad]{\sca[#1]{tcd}\XS}
\newcommand{\tep}[1][\scad]{\sca[#1]{tep}\XS}
\newcommand{\Th}[1][\scad]{\sca[#1]{th}\XS}
\newcommand{\tf}[1][\scad]{\sca[#1]{tf}\XS}
\newcommand{\tfd}[1][\scad]{\sca[#1]{tfd}\XS}
\newcommand{\tfr}[1][\scad]{\sca[#1]{tfr}\XS}
\newcommand{\tr}[1][\scad]{\sca[#1]{tr}\XS}
\newcommand{\tsvd}[1][\scad]{\sca[#1]{tsvd}\XS}
\newcommand{\us}[1][\scad]{\sca[#1]{us}\XS}
\newcommand{\vcg}[1][\scad]{\sca[#1]{vcg}\XS}

%\def\TF{\textrm{TF}\XS} \def\TFD{\foo{TFD}} \def\TFDD{\textrm{TFD}\XS}
%\def\tfd{TFD\XS}
%\def\fft{\foo{FFT}}
%\def\RSB{\textrm{RSB}\XS}

%------------------------------------------------------------

% Famous people - VIP

%\newcommand{\beaunom}[4][0]{{\taille[#1]#2 \uppercase{#3}}{\sca[#1]{#4}}\XS}
%\newcommand{\beaunom}[3][0]{#2 \sca[#1]{#3}\XS}

\newcommand{\MAL}[1][\scad]{Marc A\sca[#1]{llain}\XS}
\newcommand{\MBA}[1][\scad]{Michel B\sca[#1]{arlaud}\XS}
\newcommand{\LBE}[1][\scad]{Luigi B\sca[#1]{edini}\XS}
\newcommand{\ABL}[1][\scad]{Amandine B\sca[#1]{lanc}\XS}
\newcommand{\LBF}[1][\scad]{Laure B\sca[#1]{lanc}-F\sca[#1]{�raud}\XS}
\newcommand{\CBO}[1][\scad]{Charles A. B\sca[#1]{ouman}\XS}
\newcommand{\YBR}[1][\scad]{Yoram B\sca[#1]{resler}\XS}
\newcommand{\SBR}[1][\scad]{St�phane B\sca[#1]{rette}\XS}
\newcommand{\DBR}[1][\scad]{David B\sca[#1]{rie}\XS}
\newcommand{\JFC}[1][\scad]{Jean--Fran\c{c}ois C\sca[#1]{ardoso}\XS}
\newcommand{\HCA}[1][\scad]{Herv� C\sca[#1]{arfantan}\XS}
\newcommand{\BCH}[1][\scad]{Bernard C\sca[#1]{halmond}\XS}
\newcommand{\FCH}[1][\scad]{Fr�d�ric C\sca[#1]{hampagnat}\XS}
\newcommand{\LCH}[1][\scad]{Laurence C\sca[#1]{hatellier}\XS}
\newcommand{\PCI}[1][\scad]{Philippe C\sca[#1]{iuciu}\XS}
\newcommand{\GCI}[1][\scad]{Gabriela C\sca[#1]{iuperca}\XS}
\newcommand{\GCO}[1][\scad]{Guylaine C\sca[#1]{ollewet}\XS}
\newcommand{\PCO}[1][\scad]{Patrick C\sca[#1]{ombettes}\XS}
\newcommand{\CCO}[1][\scad]{Claude C\sca[#1]{omtat}\XS}
\newcommand{\MDA}[1][\scad]{Manuel D\sca[#1]{avy}\XS}
\newcommand{\ADE}[1][\scad]{Alexander H. D\sca[#1]{elaney}\XS}
\newcommand{\JDL}[1][\scad]{Jan De L\sca[#1]{eeuw}\XS}
\newcommand{\GDE}[1][\scad]{Guy D\sca[#1]{emoment}\XS}
\newcommand{\CDO}[1][\scad]{Christian D\sca[#1]{oncarli}\XS}
\newcommand{\FDU}[1][\scad]{Fr�d�ric D\sca[#1]{ublanchet}\XS}
\newcommand{\SDU}[1][\scad]{St�phanie D\sca[#1]{ubost}\XS}
\newcommand{\PDU}[1][\scad]{Patrick D\sca[#1]{uvaut}\XS}
\newcommand{\MFA}[1][\scad]{Marc F\sca[#1]{ayolle}\XS}
\newcommand{\DFA}[1][\scad]{Dario F\sca[#1]{arina}\XS}
\newcommand{\JFE}[1][\scad]{Jeffrey A. F\sca[#1]{essler}\XS}
\newcommand{\CFE}[1][\scad]{C�dric F\sca[#1]{�votte}\XS}
\newcommand{\MFI}[1][\scad]{Myriam F\sca[#1]{iani}\XS}
\newcommand{\IFI}[1][\scad]{Inbar F\sca[#1]{ijalkow}\XS}
\newcommand{\JJF}[1][\scad]{Jean-Jacques F\sca[#1]{uchs}\XS}
\newcommand{\LGA}[1][\scad]{Line G\sca[#1]{arnero}\XS}
\newcommand{\SGA}[1][\scad]{St�phane G\sca[#1]{autier}\XS}
\newcommand{\DGE}[1][\scad]{Di G\sca[#1]{e}\XS}
\newcommand{\JFG}[1][\scad]{Jean--Fran\c{c}ois G\sca[#1]{iovannelli}\XS}
\newcommand{\SGO}[1][\scad]{Simon G\sca[#1]{odsill}\XS}
\newcommand{\EGO}[1][\scad]{Evans G\sca[#1]{ouno}\XS}
\newcommand{\YGO}[1][\scad]{Yves G\sca[#1]{oussard}\XS}
\newcommand{\JPG}[1][\scad]{Jean-Pierre G\sca[#1]{u�don}\XS}
\newcommand{\AHE}[1][\scad]{Alain H\sca[#1]{erment}\XS}
\newcommand{\DHU}[1][\scad]{David H\sca[#1]{unter}\XS}
\newcommand{\JID}[1][\scad]{J�r\^ome I\sca[#1]{dier}\XS}
\newcommand{\CLA}[1][\scad]{Christian L\sca[#1]{abat}\XS}
\newcommand{\JFL}[1][\scad]{Jean--Fran\c{c}ois L\sca[#1]{afay}\XS}
\newcommand{\JJL}[1][\scad]{Jean--Jacques L\sca[#1]{aurin}\XS}
\newcommand{\KLA}[1][\scad]{Kenneth L\sca[#1]{ange}\XS}
\newcommand{\RLE}[1][\scad]{Richard M. L\sca[#1]{eahy}\XS}
\newcommand{\GLB}[1][\scad]{Guy Le B\sca[#1]{esnerais}\XS}
\newcommand{\ELC}[1][\scad]{\'Eric Le C\sca[#1]{arpentier}\XS}
\newcommand{\MFL}[1][\scad]{Marie--Fran\c{c}oise L\sca[#1]{ucas}\XS}
\newcommand{\SMA}[1][\scad]{Salima M\sca[#1]{akni}\XS}
\newcommand{\TMA}[1][\scad]{Thierry M\sca[#1]{artin}\XS}
\newcommand{\VMA}[1][\scad]{Vincent M\sca[#1]{azet}\XS}
\newcommand{\OMI}[1][\scad]{Olivier M\sca[#1]{ichel}\XS}
\newcommand{\AMD}[1][\scad]{Ali M\sca[#1]{ohammad--}D\sca[#1]{jafari}\XS}
\newcommand{\SMO}[1][\scad]{Said M\sca[#1]{oussaoui}\XS}
\newcommand{\LMU}[1][\scad]{Laurent M\sca[#1]{ugnier}\XS}
\newcommand{\MNI}[1][\scad]{Mila N\sca[#1]{ikolova}\XS}
\newcommand{\NNO}[1][\scad]{Nicolas N\sca[#1]{ormand}\XS}
\newcommand{\JCP}[1][\scad]{Jean-Christophe P\sca[#1]{esquet}\XS}
\newcommand{\HPL}[1][\scad]{H�l�ne P\sca[#1]{iet}-L\sca[#1]{ahanier}\XS}
\newcommand{\NOU}[1][\scad]{Nadia O\sca[#1]{udjane}\XS}
\newcommand{\JBP}[1][\scad]{Jean-Baptiste P\sca[#1]{oline}\XS}
\newcommand{\ARI}[1][\scad]{Andrea R\sca[#1]{idolfi}\XS}
\newcommand{\KSA}[1][\scad]{Ken D. S\sca[#1]{auer}\XS}
\newcommand{\MSE}[1][\scad]{Myriam S\sca[#1]{ervi�res}\XS}
\newcommand{\HSN}[1][\scad]{Hichem S\sca[#1]{noussi}\XS}
\newcommand{\CSO}[1][\scad]{Charles S\sca[#1]{oussen}\XS}
\newcommand{\ITA}[1][\scad]{Ina T\sca[#1]{aralova}\XS}
\newcommand{\ETH}[1][\scad]{\'Eric T\sca[#1]{hi�baut}\XS}
\newcommand{\ATO}[1][\scad]{Anna T\sca[#1]{onazzini}\XS}
\newcommand{\TVE}[1][\scad]{Thomas V\sca[#1]{eit}\XS}
\newcommand{\CVO}[1][\scad]{Curtis R. V\sca[#1]{ogel}\XS}
\newcommand{\JZE}[1][\scad]{Josiane Z\sca[#1]{erubia}\XS}

\def\CNRS{Centre National de la Recherche Scientifique\XS}
%\newcommand{\cnrs}[1][\scad]{\sca[#1]{cnrs}\XS}

\def\UPS{Universit� de Paris-Sud\XS}
%\newcommand{\ups}[1][\scad]{\sca[#1]{ups}\XS}
%\newcommand{\ESE}[1][\scad]{\sca[#1]{sup�lec}\XS}
\def\LSS{Laboratoire des Signaux et Syst\`emes\XS}
\newcommand{\lss}[1][\scad]{\sca[#1]{l{\footnotesize 2}s}\XS}

%\newcommand{\tutelle}[1][\scad]{\cnrs[#1]~--~\ESE[#1]~--~\ups[#1]}
%\def\ALSS{\LSS~(\tutelle)\\\adresse}
%\def\ALSSun{\LSS (\tutelle) \adresse}
%\newcommand{\adresse}[1][\scad]{\ESE[#1], Plateau de Moulon, 91192 Gif--sur--Yvette Cedex, France\XS}
%\def\adrshort{\ESE, Plateau de Moulon, 91192 Gif--sur--Yvette Cedex\XS}

\def\irccyn{Institut de recherche en communications et cybern�tique\XS}
%\newcommand{\IRCCyN}[1][\scad]{\sca[#1]{ircc}y\sca[#1]{n}\XS}
\def\adts{Analyse et d�cision en traitement du signal\XS}
\def\adtsi{Analyse et d�cision en traitement du signal et de l'image\XS}
%\newcommand{\ADTS}[1][\scad]{\sca[#1]{adts}\XS}
%\newcommand{\ECN}[1][\scad]{\sca[#1]{ecn}\XS}
%\newcommand{\EMN}[1][\scad]{\sca[#1]{emn}\XS}
\def\UDN{Universit� de Nantes\XS}
%\newcommand{\udn}[1][\scad]{\sca[#1]{udn}\XS}
%\newcommand{\tutelle}[1][\scad]{\cnrs[#1]/\ECN[#1]/\EMN[#1]/\UDN}

%\RequirePackage{color}
\RequirePackage{pifont}
%\RequirePackage{/home/idier/tex/fonts/marvosym}	% ATTENTION : interferences bizarres (eg., \implies)
\RequirePackage[alpine]{ifsym}
\def\mail{}
\def\phone{~{\large \Telefon}~}%\ding{37}
\def\fax{~\smash{\normalsize\Faxmachine}~}
%\def\helvnormal#1{{\fontfamily{phv}\fontsize{10}{14pt}\selectfont #1}}
%\def\helvfoot  #1{{\fontfamily{phv}\fontsize{9}{14pt} \selectfont #1}}
%\def\helvscript#1{{\fontfamily{phv}\fontsize{8}{14pt} \selectfont #1}}

\newcommand{\edf}[1][\scad]{\sca[#1]{�df}\XS}
\newcommand{\cea}[1][\scad]{\sca[#1]{cea}\XS}
\newcommand{\hermes}{Hermes}
\newcommand{\onera}[1][\scad]{\sca[#1]{onera}\XS}
\def\EPM{\'Ecole polytechnique de Montr�al\XS}
\def\SHFJ{Service Hospitalier Fr�d�ric Joliot\XS}
\newcommand{\shfj}[1][\scad]{\sca[#1]{shfj}\XS}
\def\LENA{Laboratoire d'Enc�phalographie et de Neurophysiologie Appliqu�e\XS}
%\newcommand{\lena}[1][\scad]{\sca[#1]{lena}\XS}
%\newcommand{\cnes}[1][\scad]{\sca[#1]{cnes}\XS}
\newcommand{\omp}[1][\scad]{\sca[#1]{omp}\XS}

\newcommand{\dea}[1][\scad]{\sca[#1]{dea}\XS}
\newcommand{\deaats}[1][\scad]{\sca[#1]{dea\,ats}\XS}
\newcommand{\deatis}[1][\scad]{\sca[#1]{dea\,tis}\XS}

%\newcommand{\ieee}[1][\scad]{\sca[#1]{ieee}\XS}
\newcommand{\ieee}{IEEE\XS}
\def\ieeeIP{\ieee Trans. on Image Processing}
\def\ieeeSP{\ieee Trans. on Signal Processing}
\def\ieeeIT{\ieee Trans. on Information Theory}
\def\ieeeSAP{\ieee Trans. on Speech and Audio Processing}
\def\ieeePAMI{\ieee Trans. on Pattern Analysis and Machine Intelligence}
\def\ICIP{IEEE Int. Conf. on Image Processing}
\def\GMIP{Graphical models and Image Processing}
\def\CRAS{Compte-rendus de l'Acad�mie des sciences}
\def\TS{Traitement du signal}

jidef_en.tex
\def\Card{\#}
\def\dd{\,d}
\def\leq{\le}
\def\geq{\ge}
\renewcommand{\det}[1]{\bars{#1}}
\def\rang{{\mathrm{rank}}}
%\def\V{,\kern-.1em}  	 % Virgule pour les nombres decimaux
\def\V{.}              % Point pour les anglo-saxons
%\def\e#1{.10^{#1}}      % Notation scientifique a la francaise
\def\e#1{\rm{e}^{#1}}  % Notation scientifique anglo-saxonne
\def\numero{n\ensuremath{^{\circ}}}

\def\et{{\normalfont and }}
\def\editorname{(ed.)}
\def\editornames{(eds.)}
\def\eqname{Eq.~}
\def\eqnames{Eqs~}
\def\tabname{Table~}
\def\tabnames{Tables~}
\def\figname{Figure~}
\def\fignames{Figures~}
\def\chapname{Chapter~}
\def\chapnames{Chapters~}
\def\parname{\S~}
\def\parnames{\S~}
\def\sectname{Section~}
\def\sectnames{Sections~}
\def\theoname{Theorem~}
\newtheorem{theorem}{Theorem}[section]
\def\gui#1{``#1''}
\def\<{``}\def\>{''}
\def\ier{st\XS}
\def\iere{st\XS}
\def\ieme{th\XS}

\def\SI{\text{if\:}}       \def\Si{\text{If\:}}
\def\ET{\text{and\:}}       \def\OU{\text{or\:}}
\def\ALORS{\text{then\:}}
\def\DOU{\text{hence\:}}  \def\Ou{\text{where\:}}
\def\QUAND{\text{when\:}}
\def\POUR{\text{for\:}}   \def\POURTOUT{\text{for all\:}}
\def\SC{\text{u.\,c.\:}}   \def\SOUSC{\text{under constraints\:}}
\def\SINON{\text{otherwise}}
\def\AVEC{\text{with\:}}
\def\DANS{\text{in\:}}
\def\IFF{\textit{if and only if}\XS}
\def\ssi{\textit{si et seulement si}\XS}

\def\post{\textit{posterior}\XS}
\def\Post{\textit{Posterior}\XS}
\def\etal{{\itshape et al.}}
\def\etcoll{\textit{et al.}}
\def\andname{and\xspace}
\def\Toeplitz{Toeplitz\XS}
\def\wrt{w.r.t.\ }

% Pour les intervalles ouverts, on utilise des parentheses
\def\croext#1{\left(#1\right)}             \def\stdcroext#1{(#1)}
\def\bigcroext#1{\bigl(#1\bigr)}           \def\biggcroext#1{\biggl(#1\biggr)}
\def\Bigcroext#1{\Bigl(#1\Bigr)}           \def\Biggcroext#1{\Biggl(#1\Biggr)}

sigproc09.tex
\documentclass[onecolumn]{elsart3p}
\usepackage[T1]{fontenc}
\usepackage{amsmath, epsfig,amssymb, graphicx}
\usepackage{setspace}
%\usepackage{psfrag}
\usepackage{algorithmic, appendix, enumitem}
\graphicspath{{figures/}}
%%   INCLUDE THE EUSIPCO STYLE BY UNCOMMENTING ONE OF THE FOLLOWING FOUR
%%%  STATEMENTS
%\usepackage{eusipco}          % 10pt Times Roman  Fonts
%\usepackage[ninepoint]{eusipco}          % 9pt Times Roman  Fonts
% \usepackage[cmfonts]{eusipco}            % 10pt Computer Modern Fonts
% \usepackage[cmfonts,ninepoint]{eusipco}  % 9pt Computer Modern Fonts
%-------------pstricks--------------------------
 %\usepackage[usenames,dvipsnames]{pstricks}
 %\usepackage{epsfig}
% \usepackage{pst-grad} % For gradients
% \usepackage{pst-plot} % For axes
 %\usepackage{rotating}

\input jidef
\input jidef_en
\input alphabet         % Definitions de macros pour les alphabets
\input beginend   % Abbreviations des \begin{} et des \end{}
\def\MT{^{-\tD}}
% def JI
%\renewcommand{\includegraphics}[1]{#1}
\usepackage{color}
%\definecolor{pink}{rgb}{1,0.2,.8}
%\definecolor{pink}{rgb}{1,0.7,.8}
\definecolor{col1}{rgb}{.7,0,0}
\def\suppJI#1{{\footnotesize \setstcolor{blue}\st{#1}}}
\def\remJI#1{{\noindent\color{blue}{{\footnotesize JI: #1}}}}
\def\addJI#1{{\noindent\color{blue}{#1}}}
%\usepackage{soul,ulem}


%\def\stepun{\textcircled{\tiny1}}
%\def\stepdeux{\textcircled{\tiny2}}
%\def\steptrois{\textcircled{\tiny3}}
%\def\stepun{1}
%\def\stepdeux{2}
%\def\steptrois{3}
%\def\tinysetminus{\text{\small$\setminus$}}
% 
\renewcommand{\algorithmiccomment}[1]{\!\texttt{\small{\% #1}}}

\newcommand{\mongraphe}[4]{
{\centering #1} \\
\begin{picture}(0,0)
\put(0,0){\makebox(0,0)[b]{\rotatebox[origin=bl]{90}{#3}}}
\put(0,0){\makebox(0,0)[bl]{#4}}
\end{picture}{\centering #2}
}
% Les interlignes
%\newcommand\singlespace{\renewcommand\baselinestretch{1.2}}
%\newcommand\doublespace{\renewcommand\baselinestretch{1.8}}
%
% Single address.
% ---------------
%\name{Di Ge, J�r�me Idier, Eric Le Carpentier}
%\address{IRCCyN (CNRS UMR 6597) 
%\\ 1, rue de la No�, BP 92101, F44321 Nantes Cedex 3, France\\
%{\small
%\url{(Di.Ge|Jerome.Idier|Eric.Le-Carpentier)@irccyn.ec-nantes.fr}}
%}
%\pagestyle{plain}
\begin{document}

\doublespace
%\ninept
%
%\maketitle
\begin{frontmatter}
%\title{Advances in sampling schemes for MCMC based blind Bernoulli-Gaussian deconvolution}
\title{Enhanced sampling schemes for MCMC based blind Bernoulli-Gaussian deconvolution}

\author[ECN]{D. Ge\corauthref{cor}}
\corauth[cor]{Corresponding author. Tel.:+33 2 40 37 69 24; fax: +33 2 40 37 69 30}
\ead{Di.Ge@irccyn.ec-nantes.fr}
\author[ECN]{J. Idier}
\ead{Jerome.Idier@irccyn.ec-nantes.fr}
\author[ECN]{E. Le Carpentier}
\ead{Eric.Le-Carpentier@irccyn.ec-nantes.fr}

\address[ECN]{IRCCyN (CNRS UMR 6597), 1 rue de la No�, BP 92101, 44321 Nantes Cedex 3, France}

\begin{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. 
%Two more efficient alternatives are studied and compared in this paper: one is derived from Bourguignon and Carfantan's~\cite{Bourguignon05} idea and  yields a generalized $K$-tuple Gibbs sampler; 
The second one, called partially marginalized sampler, is obtained by integrating the Gaussian amplitudes out of the target distribution.
%Unfortunately, for the latter case some of the sampling steps then become intractable. Therefore, our solution mixes steps in which the amplitudes are sampled with others where they are integrated out. To justify the mathematical validity of the proposed partially marginal sampler, the \emph{invariant condition} is shown to hold first, and then we provide an alternative demonstration using the equivalent standard Gibbs sampler. Another contribution of this article consists of an enhanced numerical implementation for the marginal distribution computation, in terms of both efficiency and storage. 
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.  
%In order to release samples from suboptimal convergence, a general principle is proposed to produce less correlated samples, thus making the algorithm more stable and efficient in terms of convergence rate. Inspired from the classical SMLR scheme studied by Kormylo (1982), appropriate modifications on the basic MCMC methods are proposed and then validated by simulation results that exhibit satisfactory behaviors. 
\end{abstract}

\begin{keyword}%
Blind deconvolution, Bernoulli-Gaussian model, Markov chain Monte Carlo methods 
\end{keyword}
\end{frontmatter}

%\remJI{Et le reechantillonnage de l'echelle?}

\section{Introduction}
\label{sec:intro}
The problem of the restoration of a sparse spike train \xb distorted by a linear system \hb and corrupted by noise $\epsilon$ such that $\zb = \xb \ast \hb + \epsilon$, arises in many fields such as seismic exploration~\cite{Cheng96,Kormylo83a} and astronomy~\cite{Bourguignon05}.
%It is classically dealt with using a discrete-time noisy convolution model for the observed vector $\zb=\cro{z_1,\ldots,z_N}\T$:

In this paper, we adopt a Bernoulli-Gaussian (BG) model for the spike train \xb, following~\cite{Kormylo83a} and many subsequent contributions such as~\cite{Cheng96,Champagnat96a}.
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:
\beq
\label{eqn:BG1}
P(\qb)=\lambda^L(1-\lambda)^{M-L}
\eeq
with the Bernoulli sequence $\qb=[q_1,\ldots,q_M]\T$ and $L=\sum_{m=1}^M q_m$, the number of non-zero realizations. Secondly, amplitudes $\xb=[x_1,\ldots,x_M]\T$ are assumed iid zero-mean Gaussian conditionally to \qb:
\beq
\label{eqn:BG2}
\xb\I\qb\ \sim \ \mathcal{N}(\zerob,\sigma_x^2\,\diag(\qb)),
\eeq
where $\diag(\qb)$ denotes a diagonal matrix whose diagonal is $\qb$. 
%admitting that a zero variance Gaussian distribution $\mathcal{N}(0,0)$ represents the Dirac function.
%The sparse nature of the spikes is governed by the Bernoulli law, while nonzero amplitudes are assumed iid zero-mean Gaussian.

The MCMC approach~\cite{Liu01,Robert04} is a powerful numerical tool, appropriate to solve complex inference problems such as blind deconvolution. In the field of blind BG deconvolution, Cheng \etal\ pioneered the introduction of MCMC methods~\cite{Cheng96}. They proposed to rely on a \emph{plain} Gibbs sampler, \ie with a site-by-site updating scheme for the spike train, for which their algorithm constitutes a simple and canonical example.
%that constitutes a simple and canonical example of the BG deconvolution.
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~\cite{Labat06b} 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~\cite{Veit08} 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 \qb at an acceptable speed, the Gibbs sampler tends to get stuck for many iterations around some particular configurations of \qb, often corresponding to local optimal configurations of \xb of the posterior distribution as illustrated by the example in Section~\ref{subsec:Classical Gibbs Sampling}. 
This conclusion meets Bourguignon and Carfantan's analysis~\cite{Bourguignon05}: the Markov chain equilibrates rapidly around a mode (\ie a local optimal configuration), but takes a long time to move from mode to mode.
%Guan~\cite{Guan07} studied a comparable slow convergence phenomenon of a \emph{local} chain in comparison with his \emph{small world} MCMC approach.  %Though globally more efficacious, the heavy-tailed proposals have the price to pay in terms of local mode convergence.  

%In spite of the simplicity of the above-mentioned MCMC scheme, simulation results show that different initial conditions might still yield inconsistent solutions of \xb, and this particularly in the case of standard MCMC approach as described by~\cite{Cheng96}.
In order to make up for this deficiency, our first proposition is to adopt a grouped Gibbs sampler~\cite{Liu01}, called a
%A similar idea is first appeared in Bourguignon and Carfantan~\cite{Bourguignon05}, in which spike shift is proposed and accepted if the posterior distribution is increased.
%We generalize it to a so-called
\emph{K-tuple} Gibbs sampler, where blocks of $K$ adjacent BG variables $(q_i,\ldots,q_{i+K-1})$ and their associated amplitudes $(x_i,\ldots,x_{i+K-1})$ are jointly sampled.
A comparable idea first appeared in~\cite{Chi84}, and more recently in~\cite{Bourguignon05}, under the form of a deterministic iteration aimed at producing an increase of the posterior probability, whereas our goal is to sample the latter.
%Joint updates allow more efficient jumps of the whole block as shown by simulation tests. 

We then propose a second solution based on sampling the posterior distribution \emph{marginally} to the amplitudes \xb~\cite{Ge08}. A comparable idea is found in statistical signal segmentation~\cite{Dobigeon07} where some hyper-parameters are partially marginalized.
In fact, as Liu pointed out in~\cite{Liu01}, completely integrating out some components (the Gaussian amplitudes \xb in our case), leads to a more efficient sampling scheme called \emph{collapsed Gibbs sampling}. 
%, though without mathematical justification.
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 \hb conditional on $(\zb,\qb)$ and marginally with respect to \xb. 
Our scheme solves this problem by combining a step that samples \qb marginally with respect to \xb and other sampling steps involving \xb.
Such a \emph{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.  
% since the invariant condition is shown to hold (see demonstration in~\cite{Ge08}). 
%Therefore, it is a \emph{partially marginal} sampler, though fully valid from the mathematical viewpoint, since the invariant condition is shown to hold. 
%Finally, for the sake of conciseness, it will remain implicit that the time-shift ambiguities are dealt with according to~\cite{Labat06b} in Cheng \etal's Gibbs sampler, in the generalized K-tuple Gibbs sampler as well as in the partially marginal sampler.

Simulation tests on toy examples and on the so-called Mendel's sequence~\cite{Kormylo83a} %(a benchmark sequence for BG deconvolution problems)
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.  %\remJI{(phrase provisoire)}.  

This paper is organized as follows.  
In Section~\ref{sec:formulation}, after a brief formulation of the blind BG deconvolution problem, the Gibbs sampler of the joint posterior distribution~\cite{Cheng96} is presented, and an example illustrates its inefficiency as regards the sampling of \qb.
Sections ~\ref{sec:Generalized Gibbs on K-tuple variables} and~\ref{sec:Partially marginal 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~\ref{sec:Simulation results} to compare the efficiency of the sampling schemes according to Brooks and Gelman's convergence diagnostic~\cite{Brooks98}, and conclusions are drawn in Section~\ref{sec:Conclusion}.%, and conclusions drawn in Section~\ref{sec:Conclusions}.

\section{Problem formulation}
\label{sec:formulation}
\subsection{Statistical model}
\label{subsec:Statistical model}
The mathematical model of convolution reads
\beq
\label{eqn:direct}
z_n=\sum_{k=0}^Ph_k\,x_{n-k}+\epsilon_n
\eeq
for all $n \in \{1,\ldots, N\}$, where $\zb=\cro{z_1,\ldots,z_N}\T$ denotes the observed vector, $\xb=\cro{x_1,\ldots,x_M}\T$ is the unknown spike train, $\hb=\cro{h_0,\ldots,h_P}\T$ the impulse response (IR) of the system (assumed finite here) and $\epsilonb=\cro{\epsilon_1,\ldots,\epsilon_M}\T$ an noise vector, often assumed white stationary Gaussian. 
The deconvolution problem is said \emph{blind} when $\hb$ is unknown, which is the studied case here.
Akin to~\cite{Cheng96}, the following assumptions are made:
\bit
\item $\epsilonb\sim\mathcal{N}(0,\sigma_{\epsilon}^2\Iv)$ is independent of \xb and \hb;
\item $\xb$ is a BG process defined by Eqs.\eqref{eqn:BG1}-\eqref{eqn:BG2} with $\sigma_{\xb}^2=1$;
\item $\hb \sim \mathcal{N}(0,\sigma_h^2\Iv_{P+1})$.
\eit
Let us remark that by setting $\sigma_{\xb}^2$ arbitrarily to one while imposing the Gaussian prior law \hb removes the scale ambiguity inherent to the blind deconvolution problem~\cite{Labat06b}. 
By adopting the ``zero boundary'' condition so that $M=N-P$, \zb can be rewritten in the form of a matrix multiplication $\zb = \Hv\xb + \epsilonb$, where \Hv denote the $N\times M$ Toeplitz matrix of convolution. 
It is also useful to introduce the Toeplitz matrix \Xv of size $N\times (P+1)$ such that $\Hv\xb = \Xv\hb$. 
%The remaining variables to estimate are $\Theta = \{\qb,\xb,\hb,\lambda,\sigma_{\epsilon},\sigma_{\hb}\}$, given $\zb$.
According to the Monte Carlo principle, a posterior mean estimator of $\Theta = \{\qb,\xb,\hb,\lambda,\sigma^2_{\epsilon},\sigma^2_{\hb}\}$ given \zb can be approximated by:
\begin{equation}\label{eqn:MonteCarlo}
\widehat{\Theta} = \frac1{I-J}\sum_{k=J+1}^I\Theta^{(k)},
\end{equation}
where the sum extends over the last $I-J$ samples. In the MCMC framework, the samples are generated iteratively, so that asymptotically  $\Theta^{(k)}$ follows the joint posterior distribution~\cite{Robert04}:
\bal
P(\Theta\,|\,\zb) \propto~
&g(\zb-\Hv\xb;\sigma_{\epsilon}^2\Iv_N)\, 
g(\xb;\sigma_{\xb}^2\text{diag}\{\qb\})\,%\notag\\
%&
g(\hb;\sigma_h^2\Iv_{P+1})\,P(\qb;\lambda)\,P(\sigma^2_{\hb})\,P(\lambda)\,P(\sigma^2_{\epsilon})
\label{eqn:JointPostdistribution}
\eal
where $g(\cdot;\Rv)$ denotes the centered Gaussian density of covariance \Rv. Conjugate prior laws are adopted for the last three terms $P(\sigma^2_{\hb})$,$P(\lambda)$ and $P(\sigma^2_{\epsilon})$:  
\balx
\sigma_{\hb}^2 \sim \mathit{IG}(1,1),\quad \lambda \sim \mathit{Be}(1,1), \quad \sigma_{\epsilon}^2 \sim \mathit{IG}(1,1),
\ealx
where $\mathit{IG}$ and $\mathit{Be}$ respectively represent the inverse Gamma and Beta distributions. The flat shapes of the two distributions $\mathit{IG}(1,1)$ and $\mathit{Be}(1,1)$ (a uniform distribution on $[0,1]$) can be interpreted as conveying no specific prior information of the parameters. 

\subsection{Classical Gibbs Sampling}
\label{subsec:Classical Gibbs Sampling}
A Gibbs sampler circumvents the difficulties of directly inferencing from the joint posterior distribution as in Eq.~\eqref{eqn:JointPostdistribution} and instead draws each parameter from its conditional law while all others are fixed. The principle of Cheng \etal's Gibbs sampler is given in Table~\ref{Cheng0} and a pseudocode can be found in Table~\ref{Cheng}.
%A Markov chain is generated to converge in distribution to the joint posterior distribution~\cite{Robert04}. 
\begin{table}[b]
\caption{Cheng \etal's Gibbs sampler.}
\label{Cheng0}
\centering
\medskip
\bfmi[.45\linewidth]
\btabu{ll}
%\text{\textcircled{\tiny1} }&\text{Let $\yb =(\qb,\xb)$.}\\
%&\text{for each } i=1\ldots,M,  $\\
\textcircled{\tiny1} &Let $\yb =(\qb,\xb)$. For each $i=1\ldots,M,$\\
& draw $y_i^{(k+1)}\I\yb_{1:i-1}^{(k+1)},\yb_{i+1:M}^{(k)},\hb^{(k)},\sigma_{\epsilon}^{(k)},\lambda^{(k)}, \zb$\\ %\label{AGibbs:q}
%\textcircled{\tiny2} & \xb^{(k+1)} \I\qb^{(k+1)},\hb^{(k)},\sigma_{\epsilon}^{(k)}, \zb$\\
\textcircled{\tiny2} & draw $\hb^{(k+1)} \I\xb^{(k+1)},\sigma_{\hb}^{(k)},\sigma_{\epsilon}^{(k)}, \zb$ \\%\label{Di}\\
\textcircled{\tiny3} & draw $\sigma_{\epsilon}^{(k+1)}\I\xb^{(k+1)},\hb^{(k+1)}, \zb$\\ %\label{AGibbs:eps}
\textcircled{\tiny4} & draw $\lambda^{(k+1)}\I\qb^{(k+1)}$\\
\textcircled{\tiny5} & draw $\sigma_{\hb}^{(k+1)}\I\hb^{(k+1)}$ %\label{AGibbs:end}
\etabu
\efmi
\end{table}
\begin{table}[b]
\caption{A pseudocode for Cheng \etal's Gibbs sampler (see~\cite{Cheng96} for implementation details).
%$\Iv_{P+1}$ denotes the $(P+1)$-by-$(P+1)$ identity matrix and $\hb_i$ the $i$th column of \Hv. $\eb_i=\zb-\Hv\xb+\hb_i x_i$. $\texttt{randn}(j,k)$ returns j-by-k matrix whose elements are iid centered Gaussians with variance $\sigma^2=1$.
}
\label{Cheng}
\centering
\medskip
\bfmi[.8\linewidth]
%~\\[-7mm]
\begin{algorithmic}
\STATE \COMMENT{Initialization}
\STATE $\hb\gets\zerob;~\hb(\texttt{round}(P/2))\gets1$ 
\STATE $\qb\gets\zerob;~\xb\gets\zerob$ 
\STATE $\sigma_{\xb}^2\gets1$
\STATE Sample $\sigma_{\epsilon}^2 \sim \mathit{IG}(1,1);~\sigma_{\hb}^2 \sim \mathit{IG}(1,1)$
\REPEAT 
%\STATE \texttt{\small{\% Step~1: Sample $\{q_i,x_i\}$---------------------}}
%\COMMENT{$\sigma_{\xb}^2=1$ is fixed}
\STATE \COMMENT{--------------------- Step~1: Sample $y_i = (q_i,x_i)$ ---------------------}
\STATE $\sigma_1^2\ \gets \sigma_{\epsilon}^2\sigma_{\xb}^2 / (\sigma_{\epsilon}^2 + \sigma_{\xb}^2\left\|\hb\right\|^2)$ 
\STATE $\eb\gets\zb-\Hv\xb$
\FOR{$i=1$ to $M$}
\STATE $\eb_i\gets\eb+\hb_i x_i$ \COMMENT{$\hb_i$ is the $i$-th column of \Hv}
\STATE $\mu_i \gets (\sigma_1^2 / \sigma_{\epsilon}^2) \hb_i\T\eb_i$
\STATE ${\nu_i} \gets \lambda (\sigma_1 / \sigma_{\xb})\exp\left(\mu_i^2 / (2\sigma_1^2)\right)$
\STATE $\lambda_i \gets {\nu_i}/({\nu_i}+1-\lambda)$
\STATE Sample $q_i \sim \textit{Bi}(\lambda_i)$ \COMMENT{Step 1a}
\STATE Sample $x_i \sim \Nc(\mu_i,\sigma_1^2)$ if $q_i=1$, $x_i=0$ otherwise \COMMENT{Step1b} 
\ENDFOR
\STATE \COMMENT{--------------------- Step~2: Sample $\hb$ ---------------------}
\STATE $\Sv \gets \sigma_{\epsilon}^{-2}\Xv\T\Xv+\sigma_{\hb}^{-2}\texttt{eye}(P+1)$ \COMMENT{Inverse of covariance matrix $\Rv$} 
\STATE $\Uv \gets \texttt{chol}(\Sv)$ \COMMENT{Cholesky factor of $\Sv=\Rv\M$, \ie $\Rv=\Uv\M\Uv\MT$}
\STATE $\hb \gets \Uv\backslash \left(\sigma_{\epsilon}^{-2} \Uv\T\backslash \Xv\T\zb + \texttt{randn}(P+1,1)\right)$ \quad \COMMENT{$\hb \sim \Nc(\mb,\Rv)$ with $\mb=\sigma_{\epsilon}^{-2}\Rv\Xv\T\zb$}
%\STATE Sample $\hb \sim P(\hb \I \yb, \zb,\sigma_{\hb}^2, \sigma_{\epsilon}^2) = \Nc(m,\Rv)$;
\STATE \COMMENT{--------------------- Step~3: Sample $\sigma_{\epsilon}^2$ ---------------------}
\STATE Sample $\sigma_{\epsilon}^2 \sim \text{IG}(N/2+1,\left\|\zb-\Hv\xb\right\|^2/2+1)$
\STATE \COMMENT{--------------------- Step~4: Sample $\lambda$ ------------------------}
\STATE Sample $\lambda \sim \text{Be}(1+L,1+M-L)$, with $L=\sum_m q_m$
\STATE \COMMENT{--------------------- Step~5: Sample $\sigma_{\hb}^2$ ---------------------}
\STATE Sample $\sigma_{\hb}^2 \sim \text{IG}(P/2+1,\left\|\hb\right\|^2/2+1)$
\UNTIL{Convergence}
\end{algorithmic}
%\end{algorithm}
\efmi\\
\end{table}
%
The five simulation steps are iterated until convergence toward the posterior distribution in Eq.~\eqref{eqn:JointPostdistribution}, and $\widehat{\Theta}$ is finally built according to (\ref{eqn:MonteCarlo}).
Let us remark than in Table~\ref{Cheng}, Step 1 is performed in two substeps: for each $i$, $q_i$ is first sampled conditional on $\left( \Theta \setminus (q_i, x_i), \zb \right)$, and then $x_i$ conditional on $\left( \Theta \setminus x_i, \zb\right)$. Taken jointly, these two substeps perform the sampling of $y_i=(q_i, x_i)$ conditional on $\left( \Theta \setminus y_i, \zb \right)$.

Despite its simplicity, this scheme presents several drawbacks.
In~\cite{Labat06b}, it is pointed out that time-shift ambiguities may lead to unreliable estimates, depending on the initialization of \hb.
To make up for this deficiency, an additional Metropolis-Hasting (MH) step is proposed between Steps~1 and~2 in order to allow moves to a time-shifted version of $(\yb,\hb)$. 
Another inefficiency in the Cheng \etal's Gibbs sampler has been discussed in~\cite{Veit08} with regard to the scale factor between \hb and \xb. A scale re-sampling step is proposed and tested to yield faster convergence. 
In Table~\ref{algo:TimeShift}, the modified version of Step~2 incorporates the time-shift operation and the scale re-sampling with an efficient implementation in which $f_{\text{GIG}}$ represents the \emph{Generalized Inverse Gaussian} law (see~\cite{Labat06b} and~\cite{Veit08} for more details). 
In what follows, the resulting global scheme is referred to as a hybrid sampler, since there is a Metropolis Hastings step within the Gibbs scheme and no more a pure Gibbs sampler. 
%It should be noted that given that \Rv and $\Rv'$ are Toeplitz, the time-shift compensation operation can be calculated with complexity $\Oc(P^2)$. 

%--------------------------------Time shift algorithm----------------------------------------------
\begin{table}[h]
\caption{Modified version of Step~2 to incorporate the time-shift compensation according to~\cite{Labat06b} and the scale re-sampling according to~\cite{Veit08} .
%$\texttt{circshift}(\cdot)$ defines the right circular shift operator.
}
\label{algo:TimeShift}
\centering
\medskip
\bfmi[.8\linewidth]
\begin{algorithmic}
\STATE \COMMENT{ --------------------- Time-shift compensation (Labat \etal~\cite{Labat06b})------------}
\STATE Propose $\yb'=(\qb',\xb')$ from $\yb=(\qb,\xb)$ with probability\\[1mm]
\centerline{$\pi(\yb'\I \yb)=\left\{
\barr{ll}
1-2\eta &\text{if}\; \yb' = \yb,\\
\eta &\text{if}\; \yb' = (\texttt{circshift}(\qb,1),\texttt{circshift}(\xb,1)),\\
\eta &\text{if}\; \yb' = (\texttt{circshift}(\qb,-1),\texttt{circshift}(\xb,-1)), \\
\earr
\right.
$}\smallskip
where $\eta\in(0,1/2)$, \eg $\eta=1/4$.
\STATE $\Sv \gets \sigma_{\epsilon}^{-2}\Xv\T\Xv+\sigma_{\hb}^{-2}\texttt{eye}(P+1)$ \COMMENT{Inverse of covariance matrix $\Rv$} 
\STATE $\Sv' \gets \sigma_{\epsilon}^{-2}(\Xv')\T\Xv'+\sigma_{\hb}^{-2}\texttt{eye}(P+1)$ \COMMENT{Inverse of covariance matrix $\Rv$} 
\STATE $\Uv \gets \texttt{chol}(\Sv)$
\STATE $\Uv' \gets \texttt{chol}(\Sv')$
\STATE $\rho \gets \sigma_{\epsilon}^{-4}\norm{(\Uv')\T\backslash (\Xv')\T\zb}^2 -\sigma_{\epsilon}^{-4}\norm{\Uv\T\backslash \Xv\T\zb}^2 + 2\texttt{log} |\Uv|/|\Uv'|$
\STATE \COMMENT{$\rho = (\mb')\T(\Rv')\M \mb' - \mb\T\Rv\M \mb + \log |\Rv'|/|\Rv|$}
\IF[with probability $\min\{1,\exp(\rho/2)\}$]{$2\texttt{log(rand)}<\rho$}
\STATE $\yb \gets \yb'$
\STATE $\hb \gets \Uv'\backslash \left(\sigma_{\epsilon}^{-2} (\Uv')\T\backslash (\Xv')\T\zb + \texttt{randn}(P+1,1)\right)$
\ELSE
\STATE $\hb \gets \Uv\backslash \left(\sigma_{\epsilon}^{-2} \Uv\T\backslash \Xv\T\zb + \texttt{randn}(P+1,1)\right)$
\ENDIF
\STATE \COMMENT{ --------------------- Scale re-sampling (Veit \etal~\cite{Veit08})------------}
\STATE $\lambda \leftarrow (\sum q_i -P-1)/2$
\STATE $\alpha \leftarrow |\xb|^2/\sigma_{\xb}^2, \quad \beta \leftarrow |\hb|^2/\sigma_{\hb}^2$
\STATE $s^2 \sim f_{\text{GIG}}(\lambda, \alpha, \beta)$ \quad \COMMENT{simulation by acceptation-rejection}
\STATE $\xb \leftarrow \xb\cdot s, \quad \hb \leftarrow \hb / s$
%\STATE \COMMENT{ --------------------- End of Time-shift compensation ------------}
%\STATE \COMMENT{ --------------------- Step~2-5 as in Table~\ref{Cheng} ---------------------}
%\UNTIL{Convergence}%\COMMENT{End of Main loop}
\end{algorithmic}
\efmi
\end{table}
%--------------------------------end algorithm-----------------------------------------------------

%Labat \etal's method has proven its success in the time-shift optimization of \hb within the given window length $P+1$.
Despite the improvement brought by the timeshift and scale ressampling operations in the hybrid sampler, we have identified another source of inefficiency, completely independent from the time-shift and/or scale ambiguity problem. Figure~\ref{fig:Configs} is an illustration based on simulated data generated from a single spike, convolved with an IR defined by 
%while the initial configuration in Figure~\ref{fig:Configs}(a) corresponds to a local maximum.
$$
\hb(i)=\cos\pth{(i-10)\frac{\pi}{4}}\exp\pth{-|0.225i-2|^{1.5}},\ i=0,\ldots,20,
$$
which is depicted by bullets on Figure~\ref{fig:MCMCTimeshifted}(b). A local optimal configuration of \xb (two neighboring spikes instead of a single one) is chosen as the initial state in Figure~\ref{fig:Configs}(a).
The hybrid sampler generates a Markov chain that typically takes several thousands of iterations (depending on the chosen random generator seed) to visit the optimal configuration for the first time.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[htb]
  \centering
\btabu{@{}c@{~~}c@{}}
  \figc[width=.47\linewidth]{Iter_TS0} & \figc[width=.47\linewidth]{Iter_TS1000} \\
{\footnotesize  (a) Initial configuration.} & {\footnotesize (b) 1000th iteration}\\[2mm]
  \figc[width=.47\linewidth]{Iter_TS2000} & \figc[width=.47\linewidth]{Iter_TS2415} \\
{\footnotesize  (c)  2000th iteration} & {\footnotesize (d) 2415th iteration}\\
\etabu
\caption{Sampled BG sequences obtained by the hybrid sampler on a simple example. From a local optimal configuration chosen as initial state, 
the Markov chain spends several hundreds of iterations before visiting the solution, \ie a unique spike in position 10 (marked as a bullet). 
}
\label{fig:Configs}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%----------on rajoute ici les simu de timeshift-compensated MCMC---------------------------------
%
\begin{figure}[htb]
  \centering
\setlength{\tabcolsep}{0pt}
\btabu{cc}
  \figc[height=6cm]{TSX20}%,height=3cm}
&  \figc[height=6cm]{TSIR20}%,height=3cm}
\\
  \figc[height=6cm]{TSX26}%,height=3cm}
&  \figc[height=6cm]{TSIR26}%,height=3cm}
\\
  \figc[height=6cm]{TSX32}%,height=3cm}
&
  \figc[height=6cm]{TSIR32}%,height=3cm}
\\
  {\footnotesize (a) estimation of \xb}
&  {\footnotesize (b) estimation of \hb}
%\end{minipage}
\etabu
\caption{Three different solutions of the hybrid sampler~\cite{Labat06b}, obtained by changing the random generator seed only. The actual values are marked as bullets.}
\label{fig:MCMCTimeshifted}
\end{figure}
%-----------------------------fin simu timeshift compensated MCMC---------------------------------
Figure~\ref{fig:MCMCTimeshifted} further illustrates that the hybrid sampler tends to produce unreliable estimated values. Here, Mendel's well-known BG sequence is adopted as a test signal~\cite{Kormylo83a}. The same IR is used 
%, especially for the BG sequence \xb in the case of a more difficult problem such as the well-known Mendel's test sequence~\cite{Kormylo83a}.
and the data are corrupted by Gaussian noise with $\sigma_{\epsilon}^2=4\times10^{-6}$, corresponding to a signal-to-noise ratio (SNR) of $12.80$\,dB. 
%Initialization of the Markov chain $\Theta^{(0)}$ is done according to:
%$$
%%(\sigma_{\epsilon}^{(0)},~\qb^{(0)}=\xb^{(0)} = \zerob,~\hb^{(0)}=\frac{\|\zb\|_1}M\delta(n-11),~\lambda^{(0)} = 0.1.
%\sigma_{\epsilon}^{(0)}=10^{-2}\!\!,~\qb^{(0)} = \zerob,~\hb^{(0)}=\frac{\|\zb\|_1}M\delta(n-11),~\lambda^{(0)} = 0.1.
%%\overline{|\zb|}
%$$
The three estimation results in Figure~\ref{fig:MCMCTimeshifted} are obtained from the same simulated data \zb, and the same initial states, the only difference being the seed value of the random generator. 
For each Markov chain, $1000$ samples are produced and the last 250 are averaged to compute the estimation.
Substantial variations exist from one estimated result to another, especially in the number and the positions of spikes in Figure~\ref{fig:MCMCTimeshifted}(a). Actually, it can be checked that each sequence $\stdacc{\qb^{(k)}}$ tends to remain constant for many iterations, instead of exploring the state space efficiently. It should be stressed that these results are typical, \ie they have not been selected on purpose.

This phenomenon is basically due to Step~1 of the sampler. % that explores the state space of $(\qb,\xb)$ with a low efficiency.
More precisely, the corresponding sampled chains hardly escape from local optima of the posterior distribution, because the consecutive samples are highly dependent. 
% Further in Section~\ref{sec:To a more efficient sampler} is discussed a more efficient method aiming to tackle the problem indirectly by marginalizing the amplitudes \xb out of the target posterior distribution, and by so doing releases the Markov chain from local maxima of the likelihood. Simulation comparison is available in Section~\ref{sec:Simulation results}.

%The next section introduces a Gibbs-type sampler that replaces the sampling step~\stepun by a more efficient move. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Generalized Gibbs on K-tuple variables}

\label{sec:Generalized Gibbs on K-tuple variables}
%Even after the improvements brought by Labat \etal~\cite{Labat06b} as shown in the previous section, the relative inefficiency of a plain Gibbs sampler~\cite{Cheng96} is evident.
This inefficiency in MCMC sampling schemes has been already noticed by Bourguignon and Carfantan in~\cite{Bourguignon05}. They propose a solution involving shifts of detected spikes to adjacent positions.
However, such a solution is a deterministic procedure aimed at increasing the posterior probability during the burn-in period of the Markov chain.
Therefore, it does not leave the posterior distribution invariant. In contrast, our goal here is to propose a valid sampling step that produces spike shifts.
%thus with a probability $\mu^{(k)}\rightarrow 0$ as $k\rightarrow\infty$ in a similar way as simulated annealing~\cite{Geman84}. 
%On the contrary, an exact Gibbs sampling is valid for each iteration of the Markov chain generation process and is mathematically ensured to converge to the invariant distribution. It is therefore of interest to introduce the spike-shifting in a Gibbs sampling framework.  
%However, such a solution is considered rather as a technique valid only for the burn-in period than an exact Metropolis-Hastings step and thus should be proposed with a probability $\mu^{(k)}\rightarrow 0$ as $k\rightarrow\infty$ to ensure the convergence to the invariant distribution~\cite{Bourguignon05}. $\mu^{(k)}$ should also be well tuned to ensure the non-existence of local maximum after the burn-in. 

Let us introduce the set $\Omega=\{0,\ldots,K-1\}$, and the notation $i+\omega= \{i+m \I m \in \omega\}$ for all $\omega\subset\Omega$. For any vector \vb, let us denote $\vb_\omega$ the subvector formed by the entries $v_m$, $m\in\omega$, and $\vb_{-\omega}$ the subvector formed by the remaining entries. Similarly, for any matrix \Vb, we introduce the submatrices $\Vb_\omega$ and $\Vb_{-\omega}$ formed by columns of \Vb, indexed by $\omega$ in the first case, and gathering the remaining columns in the second case.

The principle of the grouped Gibbs sampler to K-tuple variables consists in updating vector $\yb_{i+\Omega}$ for each $i=1,\ldots,M-K+1$ instead of the scalar updates of Step 1 in Table~\ref{Cheng0}.
%
%The subscript $i+\Omega$ (and $-(i+\Omega)$) denote the connected $K$-tuple index (and its complement respectively, \ie $\{1,\ldots, i-1, i+K, \ldots, M\}$). For any column vector \vb, $\vb_{i+\Omega}$ is its truncated vector with index belonging to $i+\Omega$ whereas for any matrix \Vv, $\Vv_{i+\Omega}$ is its sub-matrix whose column indexes belong to $i+\Omega$. For example, $\qb_{i+\Omega}$ (respectively, $\xb_{i+\Omega)}$) is the connected K-tuple Bernoulli variables (respectively, their associated Gaussian amplitudes). $\Hv_{i+\Omega}= [\Hv_i,\ldots,\Hv_{i+K-1}]$ denotes a $N\times K$ sub-matrix of \Hv.
%We also note $L_{i+\Omega}=\sum_{m\in i+\Omega} q_m$, the sum of the truncated Bernoulli vector $\qb_{i+\Omega}$.
%The algorithm in Table~\ref{Generalized Gibbs} inspires from the idea of~\cite{Bourguignon05} and generalizes it to a Gibbs sampling method. 
%For each $i=1,\ldots,M-K+1$, $\{\yb_{i+\Omega}\}$ is updated once at a time compared with the update of scalar quantities $\{\yb_i\}$ in the Step~1\, of Table~\ref{Cheng}. 
In particular, let us remark that such a joint update strategy allows a permutation of the Bernoulli vector $\qb_{i+\Omega}$ (from a configuration $(0,1)$ to $(1,0)$ in the case of $K=2$) within one iteration.

The resulting sampling scheme is given in Table~\ref{Generalized Gibbs}. 
The notation of $\qb_{i+\Omega}^{(k+1/2)}$ (and that of $\xb_{i+\Omega}^{(k+1/2)}$) denotes the fact that $\yb^{(k+1)}$ is latter re-sampled in step 2(a). 
Compared to the hybrid version of the original Gibbs sampler of Table~\ref{Cheng0}, only Step~1 has been modified. It is clear that in the case of $K=1$, we are driven back to the hybrid sampler.
%
\begin{table}[htb]
\caption{Grouped Gibbs sampler updating $K$-tuple variables}
\label{Generalized Gibbs}
\centering
\medskip
\bfmi[.5\linewidth]
\btabu{ll}
\textcircled{\tiny1}&For each  $i=1\ldots,M-K+1$,  \\
&(a)~draw $\qb_{i+\Omega}^{(k+1/2)}\I \yb_{-(i+\Omega)}^{(k)},\hb^{(k)},\sigma_{\epsilon}^{(k)},\lambda^{(k)}, \zb$\\
&(b)~draw $\xb_{i+\Omega}^{(k+1/2)}\I \yb_{-(i+\Omega)}^{(k)},\qb_{i+\Omega}^{(k+1/2)},\hb^{(k)},\sigma_{\epsilon}^{(k)}, \zb$\\
\textcircled{\tiny2}&(a) draw $\yb^{(k+1)}$ according to Table~\ref{algo:TimeShift}\\%\COMMENT{time-shift compensation}\\
&(b) draw $\hb^{(k+1)}\I\yb^{(k+1)},\sigma_{\hb}^{(k)},\sigma_{\epsilon}^{(k)}, \zb$\\
\textcircled{\tiny3}& draw $\sigma_{\epsilon}^{(k+1)}\I\xb^{(k+1)},\hb^{(k+1)}, \zb$\\
\textcircled{\tiny4}& draw $\lambda^{(k+1)}\I\qb^{(k+1)}$\\
\textcircled{\tiny5}& draw $\sigma_{\hb}^{(k+1)}\I\hb^{(k+1)}$
\etabu
\efmi
\end{table}

It should be noted that the complexity of the $K$-tuple sampler increases rapidly (indeed, exponentially) with $K$, while larger values of $K$ are expected to allow a more efficient exploration of the state space.
%better chances are expected to escape rapidly from local maxima, such as the case in Figure~\ref{fig:Configs}. 
The latter fact is illustrated by Table~\ref{tab:K-tupleIter}, which reports the number of iterations needed to visit the exact configuration of \qb for the first time, for different values of $K$ and of the random generator seed. For $K=1$, the required number of iterations varies in large amounts. Typically, several thousands of iterations are required. In contrast, the exact configuration of \qb is visited after a few tens of iterations for $K=2$. For $K=3$ or $K=4$, the first time visit happens immediatly in most cases.
\btabl[htb]
\caption{Number of iterations for $K$-tuple samplers to visit the true configuration for the first time in the example of Figure~\ref{fig:Configs}\label{tab:K-tupleIter}}
\centering
\medskip
\setlength{\tabcolsep}{3pt}
\btabu{|c|c|c|c|c|}
\hline
   \text{seed value}& $K=1$ & $K=2$ & $K=3$ & $K=4$ \\
\hline
$19$ &	2415 & 21  & 4 & 1 \\
\hline
$29$ &  1558 & 55 & 1 & 1   \\
\hline
$39$&  5137 & 12 & 1 & 1 \\
\hline
$49$&  7132 & 37 & 3 & 2 \\
\hline
\etabu
\etabl

In practice, the question remains to know what is the best trade-off between $K=1$ or larger values of $K$. In the first case, the convergence is slow in terms of iteration number, but the computing time for each iteration is relatively low. By increasing values of $K$, the computing time per iteration increases rapidly, while the iteration number decreases. Finding the best trade-off on a purely theoretical and general basis is probably hopeless. More modestly, we examine this question in Section~\ref{sec:Simulation results} on the sole basis of Mendel's example of Figure~\ref{fig:MCMCTimeshifted}. However, a prerequisite is to rely on an optimized implementation of the K-tuple algorithm. This is the goal of the next subsection.


%As a general observation, the convergence rate of a MCMC scheme is optimized by choosing a good trade-off between the two conflicting factors~\cite[Chapter 7]{Liu01}:
%\bit
%\item drawing each component conditional on others should be computationally simple (complexity per iteration); 
%\item the induced Markov chain should converge reasonably fast to its equilibrium law (number of iterations).
%\eit
%The motivation of the $K$-tuple sampling scheme for $K>1$ is to achieve higher convergence rate by considerably reducing iteration numbers while spending reasonably more time on each iteration.    
% include a table of iteration numbers for each K-tuple sampler
%In the next subsection, we will give the detailed implementation of the $K$-tuple $\{\yb_{i+\Omega}\}$ sampling Step~1 of Table~\ref{Generalized Gibbs}.

%Typically, for $K=M$, \qb is sampled jointly according to $\rho(\qb\I\hb,\sigma_{\epsilon},\lambda,\zb)$ while the complexity amounts to $\Oc(2^M)$. 
%------------------------------Implementation of the Generalised method---------------------------------------------------

\subsection{Numerical implementation}

Steps~2 to 5 being identical to those of the hybrid sampler, only the algorithmic implementation of the Step~1 is detailed in Table~\ref{algo:Ktuple}.
In the K-tuple sampler, Step~1 draws $(\qb_{i+\Omega},\xb_{i+\Omega})$ according to their joint conditional law:
\bal
P(\qb_{i+\Omega},\xb_{i+\Omega}\I \text{rest}) &\propto P(\zb\I \xb_{i+\Omega}, \xb_{-(i+\Omega)},\hb, \sigma_{\epsilon}^2) \, 
P(\xb_{i+\Omega}\I \qb_{i+\Omega},\sigma_x^2) \, P(\qb_{i+\Omega}\I \lambda)\notag\\
&\propto\exp\acc{-\frac1{2\sigma_{\epsilon}^2}\left\|\eb_{K,i}- \Hv_{i+\Omega}\xb_{i+\Omega}\right\|^2}
\, g(\xb_{i+\Omega};\sigma_x^2\diag(\qb_{i+\Omega})) \, 
%\lambda^{L_{i+\Omega}}(1-\lambda)^{K-L_{i+\Omega}}
\pth{\frac{\lambda}{1-\lambda}}^{\unb\T\qb_{i+\Omega}}
\label{jointlaw}
\eal
where %$L_{i+\Omega}=\unb\T\qb_{i+\Omega}$ and
%\Hv denotes the $N\times M$ Toeplitz matrix such that (\ref{eqn:direct}) also reads $\zb = \Hv\xb + \epsilonb$, $\hb_m$ the $m$th column of $\Hv$ and
%\balx
%\eb_{K,i}= \zb- \Hv\xb + \Hv_{K(i)}\xb_{K(i)},
$\eb_{K,i}= \zb- \Hv_{-(i+\Omega)}\xb_{-(i+\Omega)}$,
%\ealx
so that $\eb_{K,i}$ is a function of $(\xb_{-(i+\Omega)},\hb,\zb)$ but not of $\xb_{i+\Omega}$.

Let $\omega\subset\Omega$ be such that $\qb_{i+\omega}$ gathers the nonzero entries of $\qb_{i+\Omega}$, so that $\unb\T\qb_{i+\Omega}=\#\omega$.
From~\eqref{jointlaw}, it is straightforward to deduce that Step~1(b) amounts to draw a Gaussian law defined by $\xb_{i+\Omega\setminus\omega}=\zerob$ and, if $\omega\neq\emptyset,$ $\xb_{i+\omega}\sim\Nc(\mb_{i+\omega},\Rv_{i+\omega})$,
%\balx
%\xb_{(i+\Omega)\setminus(i+\omega)}&=\zerob,\\
%\xb_{i+\omega}&\sim\Nc(\mb_{i+\omega},\Rv_{i+\omega}),
%\ealx
where
\bal
\label{R}
\Rv_{i+\omega}\M&={\sigma_{\epsilon}^{-2}}\Hv_{i+\omega}\T \Hv_{i+\omega} + {\sigma_{\xb}^{-2}} \Iv_{\#\omega},
\\
\label{m}
\mb_{i+\omega}&={\sigma_{\epsilon}^{-2}}\Rv_{i+\omega}\Hv_{i+\omega}\T\eb_{K,i}.
\eal

By integration of~\eqref{jointlaw} with respect to $\xb_{i+\Omega}$, it is then possible to express the marginal conditional law of $\qb_{i+\Omega}$, along which the latter vector must be sampled according to Step~1(a):
$$
P(\qb_{i+\omega}=\unb,\qb_{i+\Omega\setminus\omega}=\zerob) = \frac{p_{i+\omega}}{\sum_{\omega'}p_{i+\omega'}}
$$
where
\beq
\label{p}
p_{i+\omega}=\sigma_{\xb}^{-(\#{\omega})}|\Rv_{i+\omega}|^{1/2}\exp\acc{\frac12\mb_{i+\omega}\T\Rv_{i+\omega}\M\mb_{i+\omega}}\pth{\frac{\lambda}{1-\lambda}}^{\#\omega}
\eeq

Compared to a direct implementation using~\eqref{R}-\eqref{p}, several sources of computational saving can be found. The main one exploits the fact that both \Hv and $\Hv\T\Hv$ are Toeplitz matrix in the ``zero boundary'' case of finite convolution. As a consequence, neither $\Hv_{i+\omega}$, nor $\Hv_{i+\omega}\T\Hv_{i+\omega}$, nor $\Rv_{i+\omega}$ depend on the position $i$. A shorter notation where $i+\omega$ is replaced by $\omega$ will thus be adopted. More importantly, it becomes possible to precompute and store $\Rv_{\omega}\M$ and $|\Rv_{\omega}|$ for all values of $\omega$. An even more efficient scheme manipulates Cholesky factors as detailed in Table~\ref{algo:Ktuple}, such that $\Uv_{\omega}\T \Uv_{\omega} = \Sv_{\omega}$, $\Uv_{\omega}$ being an upper-triangular matrix.

%--------------------------------Ktuple algorithm----------------------------------------------
%\begin{algorithm}
\begin{table}[h]
\caption{Implementation of $K$-tuple sampling, Step~1(a),(b) in Table~\ref{Generalized Gibbs}.}
\label{algo:Ktuple}
\centering
\medskip
\bfmi[.65\linewidth]
\begin{algorithmic}
\FORALL{nonempty $\omega \subset \Omega=\{0,\ldots,K-1\}$} %\COMMENT{$\Oc(2^K)$}
\STATE $L_{\omega} \gets \#\omega$
%\STATE $\Hv_{\omega} \gets \left[\hb_m\right]\I_{m\in 1+\omega}$
\STATE $\Sv_{\omega} \gets \sigma_{\epsilon}^{-2}\Hv_{\omega}\T\Hv_{\omega}+ \sigma_{\xb}^{-2} \texttt{eye}(L_{\omega})$ \COMMENT{$=\Rv_{\omega}\M$}
\STATE $\Uv_{\omega} \gets \texttt{chol}(\Sv_{\omega})$ \COMMENT{costs $\Oc(L_{\omega}^2)$ since $\Rv_{\omega}\M$ is Toeplitz}
%\STATE $\Wv_{\omega} \gets \text{inv}(\Uv_{\omega})$ \COMMENT{costs $\Oc(L_{\omega}^2)$ since $\Uv_{\omega}$ is triangular}
%\STATE $\Rv_{\omega} \gets \Uv_{\omega}\M\Uv_{\omega}^{-t}$ \COMMENT{costs $\Oc(L_{i+\Omega}^2)$}
\STATE $\alpha_{\omega} \gets |\Uv_{\omega}|\M$ \COMMENT{$=|\Rv_{\omega}|^{1/2}$}
%\STATE store $\left\{\Uv_{\omega}\M, |\Rv_{\omega}|^{1/2}\right\}$.
\ENDFOR
\FOR{$i=1$ to $M-K+1$}
\STATE $\eb_{K,i} \gets \zb- \Hv_{-(i+\Omega)}\xb_{-(i+\Omega)}$
\FORALL{nonempty $\omega \subset \Omega=\{0,\ldots,K-1\}$} %\COMMENT{$\Oc(2^K)$}
\STATE $c_{\omega} \gets \Uv_{\omega}\T\backslash \Hv_{\omega}\T\eb_{K,i}/\sigma_{\epsilon}^2$
\STATE $p_{i+\omega} \gets  \sigma_{\xb}^{-L_{\omega}} \alpha_{\omega} \exp(\norm{c_{\omega}}^2/2)(1/\lambda-1)^{-L_{\omega}}$ %\COMMENT{the normalization factor}
\ENDFOR
\STATE $p_{i+\emptyset} \gets 1$
\STATE \COMMENT{Step~1(a)}
\STATE sample $\qb_{i+\Omega}$ according to~\eqref{p}
%P(\qb_{i+\omega}=\unb,\qb_{i+\Omega)\setminus(i+\omega}=\zerob) \propto p_{i+\omega}$ 
\STATE \COMMENT{Step~1(b): $\xb_{i+\omega}\,\sim\,\mathcal{N}(\mb_{\omega},\Rv_{\omega})$}
%\STATE $L_{\omega} \gets \dim(\omega)$
\STATE $\xb_{i+\omega} \gets \Uv_{\omega}\backslash(c_{\omega} + \texttt{randn}(L_{\omega},1))$ 
\STATE $\xb_{i+\Omega\setminus\omega} \gets \zerob$ 
\ENDFOR
\end{algorithmic}
\efmi
\end{table}
%\end{algorithm}
%\end{algorithm}
%------------------------------------------------end algorithm--------------------------------------------------------------------------- 
%To evaluate the unnormalized marginal probability $p_{i+\omega}$ for each configuration $\omega$, $C_{\omega} = \mb_{\omega}\T\Rv_{\omega}\M\mb_{\omega}$ is obtained using the inverse of the Cholesky factors $\Uv_{\omega}\M$ (see comments in Table~\ref{algo:Ktuple}).   
%Once $\{p_{\omega}\}$ for all configurations are calculated, $\qb_{i+\Omega}$ and $\xb_{i+\Omega}$ can be re-sampled as vectors in the same manner as previously. 
%Since the Gaussian mean of $\xb_{i+\Omega}$  $\mb_{\omega} = 1/(\sigma_{\epsilon}^2)\Rv_{\omega}\Hv_{\omega}\T\eb_{i+\Omega}$
%$\Uv_{\omega}\M$ are also used to sample the Gaussian amplitudes $\xb_{i+\omega}$ as a multivariate Gaussian $\mathcal{N}(\mb_{\omega},\Rv_{\omega})$ in step~\stepun\,(b), for which
%\balx
%\Rv_{\omega}\M &= \frac1{\sigma_{\epsilon}^2}\Hv_{\omega}\T\Hv_{\omega}+ \frac1{\sigma_{\xb}^2} \Iv_{L_{\omega}} \\
%\mb_{\omega} &= \frac1{\sigma_{\epsilon}^2}\Rv_{\omega}\Hv_{\omega}\T\eb_{K,i}
%\ealx
%$$
%\xb_{i+\omega} \gets \mb_{\omega} + \Uv_{\omega}\M\texttt{randn}(L_{i+\Omega},1)=\Uv_{\omega}\M\left(\Uv_{\omega}^{-t}\Hv_{\omega}\T\eb_{K,i}/\sigma_{\epsilon}^2 + \texttt{randn}(L_{i+\Omega},1)\right)
%$$ 
%Note that $\left\{\Wv_{\omega}\, \texttt{randn}(\#\omega,1)\right\}$ generates a centered Gaussian of variance $\Rv_{\omega}$ since  $\Rv_{\omega} = \Uv_{\omega}\M\Uv_{\omega}\MT$.
%Finally it is again verified from the numerical implementation that the classical site-by-site Gibbs sampler is a special case ($K=1$) of the generalized $K$-tuple BG re-sampling algorithm in a sense that $\eb_{K,i}=\eb_i,\mb_{\omega} = \mu_i, \Rv_{\omega}=\sigma_1^2, p_{\omega}=\widetilde{\lambda}_i$ if $K=1$ and $\omega=\{0\}$ (See Algorithmin Table~\ref{algo:TimeShift}).

Using an implementation of Table~\ref{algo:Ktuple}, the iteration numbers of Table~\ref{tab:K-tupleIter} can be converted in computer time. This confirms that despite the larger time taken per iteration, the K-tuple sampler may be faster for $K>1$ than for $K=1$. However, because of the exponentially increasing cost per iteration as a function of $K$, the best trade-off cannot be reached for large values of $K$. Here and after, we have limited our study to values of $K$ not larger than 4. Some critical results of the deconvolution problem of the Mendel's sequence in Fig.~\ref{fig:MCMCTimeshifted} are reported in Tab.~\ref{tab:K1-4}: including the overall computational costs using $4$ samplers in the $K$-tuple family and the corresponding number of iterations until convergence.
More complete simulation results are provided in Section~\ref{sec:Simulation results}. 
\btabl[htb]
\caption{Convergence rate of the K-tuple ($K=1,\ldots,4$) Gibbs scheme for the example in Fig.~\ref{fig:MCMCTimeshifted} measured in iteration numbers and time in seconds. Iteration numbers are noted in each case, e.g. $4600$ for the simple hybrid sampler and $900$ for the $2$-tuple , etc. \label{tab:K1-4}}
\centering
\btabu{|c||c|c|c|c|}
\hline
		&		$K =1$ & $K =2$ & $K =3$ & $K =4$ \\
\hline
Number of iterations & 	4600 & 900 & 700& 600 \\
\hline
Computation time (s) & 915 & 202 & 285 & 452 \\
\hline
\etabu
\etabl

%Though this generalization seems direct to implement and efficient in terms of iteration numbers (see example in Table~\ref{tab:K-tupleIter}), the computational cost for each iteration grows exponentially with respect to $K$ and there exists no guideline as to how to choose a good compromise. 
%For example, some critical results of the deconvolution problem of the Mendel's sequence in Figure~\ref{fig:MCMCTimeshifted} are reported in Table~\ref{tab:K1-4}: including the overall computational costs using $4$ samplers in the $K$-tuple family and the corresponding number of iterations until convergence. %Figure~\ref{fig:K1-4}) illustrates the total computation time (measured in seconds) until convergence for the $4$ $K$-tuple samplers dealing with the same deconvolution problem in Figure~\ref{fig:MCMCTimeshifted}. 
%Considerable gain in iteration numbers (from $4600$ to $900$) is achieved by using the $2$-tuple sampler. The gain is however less significant if we further increment the $K$-tuple dimension. 
%Thus, the best compromise in the $K$-tuple sampler {\it series} is reached for $K=2$ in the given example. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\begin{figure}[htb]
%  \centering
%  \begin{sideways}\parbox{30mm}{time (s)}\end{sideways}
%  \figc[width=0.8\linewidth]{K1-4}\\
%\caption{Convergence rate of the K-tuple ($K=1,\ldots,4$) Gibbs scheme measured in seconds. Iteration numbers are noted in each case, e.g. $4600$ for the simple hybrid sampler and $900$ for the $2$-tuple , etc.}
%\label{fig:K1-4}
%\end{figure}
%\begin{table}[htb]
%\caption{Convergence rate of the K-tuple ($K=1,\ldots,4$) Gibbs scheme measured in seconds. Iteration numbers are noted in each case, \eg $4600$ for the simple hybrid sampler and $900$ for the $2$-tuple.}
%\label{tab:K1-4}
%\centering
%\medskip
%\setlength{\tabcolsep}{3pt}
%\btabu{|c||c|c|c|c|}
%\hline
%		&		$K =1$ & $K =2$ & $K =3$ & $K =4$ \\
%\hline
%Number of iterations & 	4600 & 900 & 700& 600 \\
%\hline
%Computation time (s) & 515 & 202 & 285 & 452 \\
%\hline
%\etabu
%\end{table}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% should add an example of convergence rate-versus K-tuple dimensions, maybe after the introduction of Gelman's convergence criteria in the simulation section.
%Comparing a simple hybrid sampler with a $2$-tuple sampler, 

%In summary, the $K$-tuple sampling scheme allows the spike-shifting in a strict MCMC framework and avoids the burn-in constraint. It is also proved that a $2$-tuple sampler achieves better compromise than the hybrid sampler (the case of $K=1$) for the given examples.
%More complete simulation results are provided for discussion in Section~\ref{sec:Simulation results}.
%it is still impossible to pick the best $K$ among the $K$-tuple family.  
%In Section~\ref{sec:Partially marginal Gibbs sampler} is discussed another sampling scheme to indirectly tackle the same inefficiency problem. The Markov Chain is shown to be released from local optima of the posterior distribution by marginalizing out \xb.
%Simulation comparison is available in Section~\ref{sec:Simulation results}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\begin{figure}[htb]
%  \centering
%\btabu{@{}c@{~~}c@{}}
%  \figc[width=.47\linewidth]{Iter_K=2tuple0} & \figc[width=.47\linewidth]{Iter_K=2tuple40} \\
%{\footnotesize  (a) Initial configuration.} & {\footnotesize (b) 40th iteration}\\[2mm]
% \figc[width=.47\linewidth]{Iter_K=2tuple74} & \figc[width=.47\linewidth]{Iter_K=2tuple79} \\
%{\footnotesize  (c)  74th iteration} & {\footnotesize (d) 79th iteration}\\
%\etabu
%\caption{BG sequences sampled by the K-tuple ($K=2$) Gibbs scheme in the example of Figure~\ref{fig:Configs}.}
%\label{fig:2tuple}
%\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\begin{figure}[htb]
%  \centering
%  \btabu{@{}c@{~~}c@{}}
%  \figc[width=.47\linewidth]{Iter_K=3tuple0} & \figc[width=.47\linewidth]{Iter_K=3tuple35} \\
%{\footnotesize  (a) Initial configuration.} & {\footnotesize (b) 35th iteration}\\[2mm]
%  \figc[width=.47\linewidth]{Iter_K=3tuple37} & \figc[width=.47\linewidth]{Iter_K=3tuple39} \\
%{\footnotesize  (c)  37th iteration} & {\footnotesize (d) 39th iteration}\\
%\etabu
%\caption{BG sequences sampled by the K-tuple ($K=3$) Gibbs scheme in the example of Figure~\ref{fig:Configs}.}
%\label{fig:3tuple}
%\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%----------On passe donc � un sampler plus efficace � marginaliser les ampli----------------------

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Partially marginal Gibbs sampler}

\label{sec:Partially marginal Gibbs sampler}
Here, another sampling scheme is proposed and studied, that indirectly tackles the same inefficiency by marginalizing out \xb in the step that samples \qb.

\subsection{Partial marginalization of \xb}
\label{subsec:Another Gibbs sampler}% we present here the modified Gibbs sampler, in order to produce more efficient samples
%The motivation to exclude \xb from the Gibbs sampler is twofold. Firstly, the estimation of $\xb$ conditional on $(\widetilde{\Theta},\zb)$ with $\widetilde{\Theta}=(\qb,\hb,\sigma_{\epsilon},\sigma_{h}, \lambda)$ is a linear estimation problem, easy to solve afterward. 
In principle, one alternative to the approach in Eq.~\eqref{eqn:MonteCarlo} could be to generate a Markov chain on $\widetilde{\Theta}=(\qb,\hb,\sigma_{\epsilon}^2,\sigma_{h}^2, \lambda)$ marginalizing \xb, and then build the estimates in two steps: (1) the marginal maximum \apost estimate on each binary parameter $q_i$ and the posterior mean estimate on the continuous parameters $\widetilde{\Theta}\setminus \qb$, (2) the posterior mean estimate of \xb conditional on $(\widetilde{\Theta},\zb)$, that is a linear estimation problem. 
Furthermore, a Gibbs sampler with target distribution $P(\widetilde{\Theta}\I\zb)$ is likely to be more efficient than the hybrid sampler, particularly with respect to the sampling of \qb. Theoretical foundations are available in~\cite{Liu94b} and~\cite[Chapter 6.7]{Liu01} regarding the convergence rate of the so called \emph{collapsed Gibbs sampler}, with the marginalized equilibrium distribution $P(\widetilde{\Theta}\I\zb) = \int P(\Theta\I\zb)d\xb$. According to~\cite[Chapter 6.7]{Liu01}, the collapsed Gibbs sampler produces a substantial gain in terms of convergence rate if the marginalized parameters form highly dependent pairs with other parameters of the Markov chain~(see also~\cite{Cappe99} and references therein). This is exactly the case of the pair $\{\xb,\qb\}$ in the spike train deconvolution problem. Park \etal~\cite{Park09} also illustrates through three examples that such a sampler must be implemented with care to be sure that the desired stationary distribution is preserved. 

%From a more practical point of view, the convergence rate should be accelerated by correcting the pathological behavior as illustrated in Figure~\ref{fig:MCMCTimeshifted}.

%To measure the Markov chain's convergence rate, the most widely used property is its geometric ergodicity. 
%Functional analysis shows that given a {\it Hilbert space} $L^2(\rho)=\{h(\cdot)\I \left\|h\right\|_{\rho}<\infty\}$ equiped with an inner-product $\left\langle h,g \right\rangle=E_{\rho}\{h(\cdot)g(\cdot)\}$, the Gibbs kernel $\kappa(\Theta,\Theta')$ induces a contraction operator on $L^2(\rho)$ such that
%\beq\label{eq:FO}
%\FD h(\Theta)=\int\kappa(\Theta,\Theta')h(\Theta')d\Theta'=E\{h(\Theta^{(1)})\I\Theta^{(0)}=\Theta\}.
%\eeq
%If $\kappa(\Theta,\Theta')$ is reversible with respect to the stationary distribution $\rho$, then the operator $\FD$ is {\it self-adjoint}, \ie $\left\langle \FD h,g \right\rangle=\left\langle h,\FD g \right\rangle$ and has thus real spectrum in $[-1,1]$. Its second largest eigenvalue (in absolute value) characterizes the geometric ergodicity (convergence rate) of the Markov chain~\cite{Liu94b}. 
%
%Then consider a subspace $L_0^2(\rho)=\{h(\cdot)\in L^2(\rho)\I E_{\rho}\{h\}=0\} \subset L^2(\rho)$. Clearly, $L_0^2(\rho)$ is also a {\it Hilbert space} and is invariant under the operator $\FD_0$ induced by $\FD$. Then the largest eigenvalue of $\FD_0$ (its spectral radius $r_0$) is the second largest of those of $\FD$ (both in absolute value). A general relationship between the norm and the spectral radius of an operator is :
%$$
%\lim_{n\rightarrow\infty}\left\|\FD_0^n\right\|^{1/n}=r_0,
%$$ 
%suggesting that by comparing the norms of the induced operator $\FD_0$, one can measure the convergence rate of different Markov chains.
%\begin{theorem}\label{theo:Norm}
%Let $\FD_s$ (respectively, $\FD_c$) denote the operator 
%\beq\label{eq:FO}
%\FD h(\Theta)=\int\kappa(\Theta,\Theta')h(\Theta')d\Theta'=E\{h(\Theta^{(1)})\I\Theta^{(0)}=\Theta\}.
%\eeq
%defined by the Gibbs kernel for the Markov chain $\{\Theta^{(t)}=(\theta_1^{(t)},\ldots,\theta_d^{(t)})\}$ (respectively, $\{\widetilde{\Theta}^{(t)}=(\theta_1^{(t)},\ldots,\theta_{d-1}^{(t)})\}$). The two stationary distributions satisfy $\rho(\widetilde{\Theta}) = \int\rho(\Theta)d\theta_d$, then the operator norms (geometric ergodicity) are ordered as
%$$
%\left\|\FD_c\right\|\leq\left\|\FD_s\right\|.
%$$
%\end{theorem}
%The proof is detailed by Liu \etal in~\cite{Liu94b}. 
%On the contrary, marginalizing \xb in step~\stepdeux(a) allows to escape more easily from a local maximum (configuration A) and in the mean time keeps the Metropolis-Hasting structure in lieu of a burn-in technique. A typical sampling sequence by marginalizing \xb is shown in Section~\ref{sect:NewGibbs} after introduction of the proposed Gibbs-type sampler. 
 
%Another issue relating the real computational convergence speed and the convergence rate

%\subsection{Proposed sampler}\label{sect:NewGibbs}
However, marginalizing \xb from the posterior distribution by $P(\widetilde{\Theta}\I\zb) = \int P(\Theta\I\zb)d\xb$ leads to practical difficulties.
Technically, the marginalization of \xb is only feasible on conditional laws of $(\qb, \lambda, \sigma_{\hb}^2)$. On the contrary, the conditional sampling of \hb becomes extremely difficult if \xb is integrated out, since $P(\hb\I\qb,\sigma_{\hb}^2,\sigma_{\epsilon}^2, \zb)$ is a multivariate, non Gaussian law with a complex structure. The same observation is made on the sampling of $\sigma_{\epsilon}^2$ if \xb is marginalized out.
Instead of a plain collapsed Gibbs sampler on $\widetilde{\Theta}$, the sampling scheme in Table~\ref{Di} is proposed to circumvent the direct conditional sampling of $\hb$ and $\sigma_{\epsilon}^2$.
%The six steps are iterated until convergence toward the posterior distribution $P(\Theta \I \zb)$. 
Compared with the hybrid sampler, the main difference appears in Step~1, where \xb has been analytically integrated out. Subsection~\ref{subsec:Marginal posterior distribution} provides an efficient way to implement Step~1. 
\begin{table}[hbt]
\caption{The proposed partially marginal sampler.}
\label{Di}
\centering
\medskip
\bfmi[.5\linewidth]
\btabu{ll}
\textcircled{\tiny1} &(a) for each $i=1\ldots,M$,  \notag\\
&draw $q_i^{(k+1/2)}\I\qb_{1:i-1}^{(k+1/2)},\qb_{i+1:M}^{(k)},\hb^{(k)},\sigma_{\epsilon}^{(k)},\lambda^{(k)}, \zb$\\ %\label{AGibbs:q}
                     &(b) draw $\xb^{(k+1/2)} \I\qb^{(k+1/2)},\hb^{(k)},\sigma_{\epsilon}^{(k)}, \zb$\\
\textcircled{\tiny2}~&(a) draw $\yb^{(k+1)}$ according to Table~\ref{algo:TimeShift}\\%\COMMENT{time-shift compensation}\\
                     &(b) draw $\hb^{(k+1)} \I\xb^{(k+1)},\sigma_{\hb}^{(k)},\sigma_{\epsilon}^{(k)}, \zb$ \\%\label{Di}\\
\textcircled{\tiny3} &draw $\sigma_{\epsilon}^{(k+1)}\I\xb^{(k+1)},\hb^{(k+1)}, \zb$\\ %\label{AGibbs:eps}
\textcircled{\tiny4} &draw $\lambda^{(k+1)}\I\qb^{(k+1)}$\\
\textcircled{\tiny5} &draw $\sigma_{\hb}^{(k+1)}\I\hb^{(k+1)}$ %\label{AGibbs:end}
\etabu
\efmi
\end{table}

%On the other hand, in step~\stepdeux, $\xb$ only plays the role of an auxiliary quantity, since it does not belong to the Markov chain anymore. Steps~\stepdeux(b),\,(c) are structurally identical to steps~\stepdeux, \steptrois\ in Table~\ref{Cheng}. In order to justify that the sampling scheme of Table~\ref{Di} is mathematically valid, we let $\widetilde{\Theta}'=(\qb,\hb',\sigma_{\epsilon}',\sigma_{h}, \lambda)$ and check the \emph{invariant condition}, \ie that
%\bal
%\label{eq:InvCond}
%\int_{\hb,\sigma_{\epsilon}}\kappa_2(\hb',\sigma_{\epsilon}'\I\widetilde{\Theta},\zb)\,\rho(\widetilde{\Theta}\I\zb)\,d\hb d\sigma_{\epsilon}
%=\rho(\widetilde{\Theta}'\I\zb)
%\eal
%holds for step~\stepdeux, where $\rho$ and $\kappa_2$ denote the invariant target distribution and the transition kernel corresponding to the step~\stepdeux, respectively. Furthermore, the analytical expression for the transition kernel writes:
%\bal
%\kappa_2&(\hb',\sigma_{\epsilon}'\I\widetilde{\Theta},\zb)=
%\int_{\xb}\rho(\xb\I \widetilde{\Theta},\zb)\,\rho(\hb' \I \xb,\widetilde{\Theta}\setminus\hb,\zb)\,\rho(\sigma_{\epsilon}'\I \xb,\widetilde{\Theta}' \setminus \sigma_{\epsilon}',\zb)d\xb.  \label{eq:noyau}
%\eal
%We emphasize that Eq.~\eqref{eq:noyau} is not obtained by Monte Carlo integration. It is rather the exact expression of the transition kernel of steps~\stepdeux\, to update the pair $\{\hb,\sigma_{\epsilon}\}$ while \xb is no longer retained in the chain. 
%The exact proof of Eq.~\eqref{eq:InvCond} is given by Ge \etal in~\cite{Ge08} and in this paper 
We give here a mathematical justification that the proposed partially marginal sampler identifies with a Gibbs sampler of the joint posterior law (Eq.~\eqref{eqn:JointPostdistribution}, that includes \xb), for a particular scanning scheme. The techniques are similar to those denoted as \emph{marginalization} and \emph{trimming} in~\cite{Park09}.

Let $(q_1,\xb), (q_2,\xb),\ldots,(q_M,\xb), \hb,\sigma_{\epsilon}^2,\lambda,\sigma_{\hb}^2$ be the scanning sequence of the Gibbs sampler. For the sampling of couples $(q_i,\xb)$, a two-stage procedure is considered by applying the Bayes rule $P(\xb, q_i \I \Theta \setminus \{\xb, q_i\} ,\zb) = P(q_i \I \Theta \setminus \{\xb, q_i\} ,\zb)P(\xb \I \Theta \setminus \xb, \zb)$: 
\bit
\item draw $q_i$ according to $P(q_i \I \qb_{1:i-1},\qb_{i+1:M},\hb,\sigma_{\epsilon}^2,\lambda,\zb)$;
\item draw \xb according to $P(\xb \I \qb,\hb,\sigma_{\epsilon}^2,\lambda,\zb)$.
\eit

The second stage is therefore repeated $M$ times within each iteration of the sampler, while all but the last one are useless since the correspondingly sampled values of \xb do not enter any subsequent operation. The $M-1$ corresponding sampling operations can thus be skipped and the resulting sequence of operations coincides exactly with the partially marginal sampler of Table~\ref{Di}. 
Another mathematical proof is given in~\cite{Ge08} by considering Table~\ref{Di} as a collapsed Gibbs sampler with  target distribution $P(\widetilde{\Theta}\I\zb)$ and then checking the \emph{invariant condition} of the related Markov chain. 
%Since each $q_i$ is sampled marginally to \xb,  
%Nevertheless, the demonstration using such a Gibbs sampler of a particular scanning scheme does not explain how should the partially marginal sampler outperform the hybrid sampler, especially on the issue of faster convergence of \qb.  

In the same conditions than that of Table~\ref{tab:K-tupleIter}, the partially marginal sampler takes about $20$ iterations.
Figure~\ref{fig:Escape} illustrates the way it escapes from a local optimal configuration within an acceptable number of iterations in the example of Figure~\ref{fig:Configs}. The true configuration is reached after 20 iterations. Moreover, it is observed from the configurations obtained at the 18th and 19th iterations that our scheme is able to radically modify \xb in one single step, a characteristic directly related to the marginalization of \xb in Step~1(a) of the sampler.

%The real spike is positioned at index $10$, while the chain starts from a local maximum, similar to the configuration A. From (b) to (c) in Figure~\ref{fig:Escape}, we note that the correct detection of $q_{10}$ drastically modifies the neighborhood amplitudes \xb, yielding final convergence in one more step.
\begin{figure}[htb]
  \centering
%\setlength{\tabcolsep}{0pt}
\btabu{@{}c@{~~}c@{}}
%  {\footnotesize initial configuration (1st iteration)}  &  {\footnotesize 18th iteration}\\
  \figc[width=.47\linewidth]{Iter3} & \figc[width=.47\linewidth]{Iter18} \\
{\footnotesize  (a) 1st iteration} & {\footnotesize (b) 18th iteration}\\[2mm]
%  {\footnotesize 19th iteration}  &  {\footnotesize 20th iteration}\\
  \figc[width=.47\linewidth]{Iter19} & \figc[width=.47\linewidth]{Iter20} \\
{\footnotesize  (c)  19th iteration} & {\footnotesize (d) 20th iteration}\\
\etabu
\caption{Sampled BG sequences obtained by our partially marginal scheme in the example of Figure~\ref{fig:Configs}.
The Markov chain escapes rapidly the initial configuration neighborhood.
%A local maximum configuration has been choosen as initial state to test its ability of visiting other probable states.
}
\label{fig:Escape}
\end{figure}

%Typically, theorem~\ref{theo:Norm} shows that in terms of number of iterations, the Markov chain $\{\widetilde{\Theta}^{(t)}\}$ converges faster than $\{\Theta^{(t)}\}$. However, the complexity per iteration should also be taken into account to evaluate the overall efficiency.
%The next section deals with the issue of the complexity per iteration, whereas the overall computational costs are compared in Section~\ref{sec:Simulation results} using simulation tests.


%---------------------------------------Marginal Post law------------------------------------------------------
\subsection{Marginal posterior distribution}
\label{subsec:Marginal posterior distribution}
%the manipulation of the Core algorithm as in SMLR
In this subsection, a numerical implementation of Step~1 of Table~\ref{Di} is proposed, inspiring from
the recursive method in~\cite{Champagnat96a} to evaluate $P(q_i=0,1|\widetilde{\Theta}\setminus q_i;\zb)$ sequentially. While the latter is based on the storage and update of an $L \times L$ ($L=\sum_i q_i$) matrix, we introduce an even less burdensome strategy by handling the Cholesky factor matrix. The issue of the complexity per iteration is dealt here, whereas the overall computational costs are compared in Section~\ref{sec:Simulation results} using simulation tests.
%Our contribution is to further adapt the recursive method to scheme~(\ref{Di}) by only updating and storing the Cholesky factor, from which one can directly sample \xb in step~\stepdeux(a).  
 
%Let \Hv denote the $N\times M$ Toeplitz matrix such that (\ref{eqn:direct}) also reads $\zb = \Hv\xb + \epsilonb$.   
First, the conditional posterior distribution of $q_i$ takes the following form, as detailed in~\cite{Champagnat96a}:
\balx%\label{eqn:postdq}
P(q_i|\widetilde{\Theta}\setminus q_i; \zb)%&\propto\rho(\zb|\qb_{-i},q_i,\hb,\sigma_{\epsilon})\rho(q_i|\lambda)\nonumber\\
\propto |\Bv|^{-1/2}\exp\left(-\frac12\zb\T\Bv\M\zb\right)\pth{\frac{\lambda}{1-\lambda}}^{q_i}\nonumber
\propto \exp\left(-f(q_i)/2\right)
\ealx
where
\bal
\Bv &= \Hv\diag(\qb)\Hv\T + \sigma_{\epsilon}^2\Iv_N, \label{eqn:B}\\
f(q_i) &= \zb\T\Bv\M\zb + \log\bars{\Bv} + 2q_i\log(1/\lambda-1). \label{eqn:f}
\eal
%Eq.~(\ref{eqn:postdq}) can be evaluated from $\frac{\rho(q_i=1|\widetilde{\Theta}\setminus q_i; \zb)}{\rho(q_i=0|\widetilde{\Theta}\setminus q_i; \zb)}$~\cite{Kormylo83a,Champagnat96a}. 
%We define
%\bal
%L(q_i|\widetilde{\Theta}\setminus q_i; \zb)&= 2\log\rho(q_i|\widetilde{\Theta}\setminus q_i; \zb)+N\log 2\pi \nonumber\\
%& = -\zb\T\Bv\M\zb - \log|\Bv| - 2q_i\log(1/\lambda-1), \nonumber
%\eal
%and $\Delta L=L(q_i=1|\widetilde{\Theta}\setminus q_i; \zb)-L(q_i=0|\widetilde{\Theta}\setminus q_i; \zb)$ is computed to obtain $\rho(q_i|\widetilde{\Theta}\setminus q_i; \zb)$,
After normalization, the marginal probability of $q_i$ reads:
\beqx
P(q_i\I\widetilde{\Theta}\setminus q_i, \zb)=\pth{1+\exp\left(-%\frac1{2}
(f(1-q_i)-f(q_i))/2\right)}\M,
%\rho(q_i=0|\widetilde{\Theta}\setminus q_i; \zb)=1-\rho(q_i=1|\widetilde{\Theta}\setminus q_i; \zb).
\eeqx
which reduces to the evaluation of $f(1-q_i)-f(q_i)$, using the following steps~\cite{Champagnat96a}:
\bal
\tau_i &= \delta_i + \sigma_{\epsilon}^{-2} \hb_i\T\widetilde{\Bv}_0\M \hb_i \\
\widetilde{\Bv}_i\M &= \widetilde{\Bv}_0\M - \sigma_{\epsilon}^{-2} \widetilde{\Bv}_0\M\hb_i \tau_i\M \hb_i\T \widetilde{\Bv}_0\M \label{eqn:deltaL}\\
\stdbars{\widetilde{\Bv}_i\M} &= \stdbars{\widetilde{\Bv}_0\M}\delta_i\tau_i \label{eq:DetBi}
\eal
with $\widetilde{\Bv} = \Bv/\sigma_{\epsilon}^2$. $\delta_i=\pm1$ depends on whether $1$ is added or removed at $q_i$ and $\widetilde{\Bv}_i$ and $\widetilde{\Bv}_0$ differ only at $q_i$. Let us note that $\widetilde{\Bv}_0, \widetilde{\Bv}_i > 0$, so that it can be deduced from~\eqref{eq:DetBi} that $\tau_i$ has the same sign as $\delta_i$.
%$\widetilde{\Bv}_i$ and $\widetilde{\Bv}_0$ differ only at component $q_i$.

Further simplifications are also introduced in~\cite{Champagnat96a} by exploiting the sparse nature of $\widetilde{\Bv}$ and applying the matrix inversion lemma. Noticing that the rank of $\Hv\diag(\qb)\Hv\T$ is only $L$, $\widetilde{\Bv}$ takes the alternate form $\widetilde{\Bv} = \sigma_{\epsilon}^{-2}\Gv\Gv\T+\Iv$, where $\Gv=\Hv\Dv$ is full rank, and $\Dv$ made of the nonzero columns of $\diag(\qb)$. By applying the matrix inversion lemma, we have
\bal
\widetilde{\Bv}\M = \Iv - \sigma_{\epsilon}^{-2} \Gv\Cv\M\Gv\T \label{eqn:Binv}
\eal
where $\Cv = \sigma_{\epsilon}^{-2}\Gv\T\Gv+\Iv$ is a $L\times L$ matrix.
Therefore, $\Cv\M$ ($L\times L$) can be stored and updated instead of $\widetilde{\Bv}\M$ ($N\times N$). %which amounts to $\Oc(L^2)$ operations instead of $\Oc(N^2)$.
In the case of adding a spike at $q_i$, the formula for the update of $\Cv_i\M$ up to a matrix permutation operation is the following: 
\bal
\bb &= -\sigma_{\epsilon}^{-2} \tau_i\M \Cv_0\M\Gv_0\T\hb_i \label{eqn:updateb} \\
\Cv_i\M &= \left[
\barr{cc}
\Cv_0\M+\bb\tau_i\bb\T & \bb \\[3pt]
\bb\T & \tau_i\M
\earr\label{eqn:updateC}
\right]
\eal
The case of removing a spike is straightforward since we have the following:
\bal \label{eq:UpdateCminus}
\Cv_0\M &= \left[
\barr{cc}
\Cv_i\M+\bb\tau_i\bb\T & \bb \\[3pt]
\bb\T & \tau_i\M
\earr
\right]
\eal
Notice that in the case of spike removal, \bb and $\tau_i$ are extracted from $\Cv_0\M$ rather than calculated as in Eq.~\eqref{eqn:updateb}.  
These operations are invariant up to a matrix permutation in cases of adding/removing a spike at an arbitrary location.  
%plus a matrix permutation operation and $\bb = -\sigma_{\epsilon}^{-2} \tau_i\M \Cv_0\M\Gv_0\T\hb_i$.
%On the other hand, it can be shown that $f(1-q_i)-f(q_i)$ takes the following form:

Finally, we propose to further reduce the computation and memory load, by updating and storing the Cholesky factor $\Fv$ of $\Cv\M$ ($\Fv\T\Fv = \Cv\M$ for an upper-triangular matrix \Fv). 
It is derived from Equations~\eqref{eqn:deltaL} and~\eqref{eqn:Binv} that:
\balx
%\\
\tau_i&=\delta_i+\sigma_{\epsilon}^{-2}\norm{\hb_i}^2-\sigma_{\epsilon}^{-4}\stdnorm{\Fv_0\Gv_0\T\hb_i}^2,\\
\zb\T\widetilde{\Bv}_0\M\hb_i&=\zb\T\hb_i-\sigma_{\epsilon}^{-2}(\Fv_0\Gv_0\T\zb)\T(\Fv_0\Gv_0\T\hb_i)
\ealx
and from Eq.~\eqref{eqn:f} and~\eqref{eqn:deltaL} that:
\bal\label{eq:deltaF}
f(1-q_i)-f(q_i)= \log(\delta_i\tau_i)-\sigma_{\epsilon}^{-4}\tau_i\M(\zb\T\widetilde{\Bv}_0\M\hb_i)^2+2\delta_i\log\pth{\UNsur{\lambda}-1}
\eal

Two extra advantages of the Cholesky decomposition are:
\bit  
\item since the conditional law in Step~1(b) reads 
\beq\label{eq:UpdateX}
(\xb\I\wt{\Thetab},\zb)\, \sim \,\Nc(\sigma_{\epsilon}^{-2}\Cv\M\Gv_0\T\zb,\Cv\M),
\eeq 
sampling \xb costs $\Oc(L^2)$ using \Fv (See Table~\ref{algo:Marginal}); 
\item the incremental update of $\Fv$ also costs $\Oc(L^2)$ using a {\it rank-$1$ Cholesky update} method. 
The examples of adding and removing a spike at the last location is given in the following for notational simplicity. It is however detailed in the appendix that unlike adding a spike at at an arbitrary location, the case of removing a spike at an arbitrary location requires two operations of Cholesky rank-$1$ update. 

For adding a spike at the last location, $\tau_i>0$ and equation~\eqref{eqn:updateC} also reads:
\beq\label{eq:CholUpdateC}
\Cv_i\M = \cro{
\barr{c@{\qquad}c}
\Fv_0\T & \zerob \\
\zerob\T & 0
\earr
}
\left[
\barr{c@{\qquad}c}
\Fv_0 & \zerob \\
\zerob\T & 0
\earr
\right]+\left[
\barr{c}
\bb\tau_i^{1/2} \\
\tau_i^{-1/2} 
\earr
\right]
\left[
\barr{c}
\bb\tau_i^{1/2} \\
\tau_i^{-1/2}
\earr
\right]\T
\eeq
\eit 
In the case of removing a spike at the last location, from Eq.~\eqref{eq:UpdateCminus} and the definition of Cholesky decomposition, we have:
\bal
\Fv_0\T \Fv_0 &= \left[
\barr{cc}
\Fv_i\T\Fv_i+\bb\tau_i\bb\T & \bb \\
\bb\T & \tau_i\M
\earr
\right] \\ 
\Fv_0 &= \left[
\barr{cc}
\Sv_0 &  \ab \\
\zerob\T  &  m
\earr
\right]
\eal
where $\Sv_0 $ is also an upper triangular matrix. It is then directly deduced that: 
\bal\label{eqn:UpdateFb}
\Sv_0\T \ab &= \bb \\
\ab\T \ab + m^2 &= \tau_i\M \\
\Sv_0\T \Sv_0 &= \Fv_i\T\Fv_i +\bb\tau_i\bb\T	 
\eal
Let $\texttt{cholupdate}(\Av, \db, '\pm')$ be the method that updates in $\Oc(L^2)$ the Cholesky factor of $\Av\T\Av\pm \db\db\T$, where \Av is already an upper triangular matrix. Then, $\Fv_i = \texttt{cholupdate}(\Sv_0, \bb\tau_i^{1/2}, '-')$. 
%\bit
%\item Evaluate $\rho(q_i=0,1\I\widetilde{\Theta}\setminus q_i; \zb)$ using \eqref{deltaf}, and sample $q_i^{(k+1)}$ accordingly;
%%\item Sample $u\sim\Uc(0,1)$ and $q_i^{(k+1)}=1$ if $u\leq\rho(q_i=1\I\widetilde{\Theta}\setminus q_i; \zb)$, $q_i^{(k+1)}=0$ otherwise;
%\item Update $\Fv$ if needed, \ie
%%(the Cholesky factor of $\Cv\M$ in Eq.~(\ref{eqn:updateC}))
%if $q_i^{(k+1)}\ne q_i^{(k)}$.%; otherwise $\Fv_i=\Fv_0$
%\eit
%--------------------------------Marginalized algorithm----------------------------------------------
\begin{table}[htb]
\caption{Step~1 of the partially marginal sampling scheme. The update of \Fv is given here up to a matrix permutation.}
\label{algo:Marginal}
\centering
\medskip
\bfmi[.9\linewidth]
\begin{algorithmic}
\STATE Initialize \Fv and \Gv
\STATE \COMMENT{Step~1(a): sequential sampling of $q_i$}
\FOR{$i=1$ to $M$}
\STATE $\delta_i \gets (-1)^{q_i}$
\STATE $\tau_i \gets \delta_i+\sigma_{\epsilon}^{-2}\norm{\hb}^2-\sigma_{\epsilon}^{-4}\stdnorm{\Fv\Gv\T\hb_i}^2$
\STATE $\phi \gets \zb\T\hb_i-\sigma_{\epsilon}^{-2}(\Fv\Gv\T\zb)\T(\Fv\Gv\T\hb_i)$ \quad \quad \COMMENT{$\phi=\zb\T\Bv\hb$}
\STATE $\Delta f \gets \texttt{log}(\delta_i\tau_i)-\sigma_{\epsilon}^{-4}\tau_i\M\phi^2+2\delta_i\log\pth{1/{\lambda}-1}$ \quad \COMMENT{see Eq.\eqref{eq:deltaF}}
\STATE sample $u \sim \Uc([0,1])$ \COMMENT{Uniform distribution in $[0,1]$}
\IF{$u > (1+\exp(-\Delta f/2))\M$} 
\STATE $q_i \gets q_i + \delta_i$ \quad \COMMENT{$q_i + \delta_i$ is accepted}
\IF{$\delta_i=1$}
\STATE $\bb \gets -\sigma_{\epsilon}^{-2} \tau_i\M \Fv\T\Fv\Gv\T\hb_i$ \quad \COMMENT{see Eq.~\eqref{eqn:updateb}}
\STATE $\Fv \gets \texttt{cholupdate}([\Fv\ \zerob; \zerob\T\ 0],[\bb \tau_i^{1/2};\tau_i^{-1/2}], '+')$ \quad \COMMENT{see Eq.\eqref{eq:CholUpdateC}}
\ELSE[the case of $\delta_i= -1$ removing a spike, see appendix for detail]
%%\STATE 
%\STATE 
%$
%\left[
%\barr{ll}
%\Sv &  \ab \\
%\mathbf{0}\T  &  m
%\earr
%\right] \leftarrow \Fv
%$
%\STATE $\bb \gets \Sv\T \ab,$ \quad $\tau_i\M \gets \ab\T \ab + m^2$ \quad \COMMENT{see Eq.\eqref{eqn:UpdateFb}}
%%\STATE $[\Fv\ \zerob; \zerob\T\ 0] \gets \texttt{cholupdate}(\Fv, \sqrt{\tau_i}[\bb;1/\tau_i], '-')$ \quad \COMMENT{$\Oc(L^2)$}
%\STATE $\Fv \gets \texttt{cholupdate}(\Sv, \tau_i^{1/2}\bb, '-')$ \quad \COMMENT{$\Oc(L^2)$}
\STATE $L \gets \text{sum}(\qb)$
\STATE $\bb \gets \Fv\T\Fv(:,i), \vb \gets \bb_{-i}, \tau\M \gets \bb_i$
\IF[other than the last location, an extra update needed]{$i< L$}
\STATE $\eb \gets \Fv(i, 1+i:L)\T$
\STATE $\Fv(1+i:L, 1+i:L) \gets \texttt{cholupdate}(\Fv(1+i:L, 1+i:L), \eb, '+')$ \quad \COMMENT{$\Oc((L-i)^2)$}
\ENDIF
\STATE $\Fv(i, :) \gets [], \Fv(:, i) \gets []$ \COMMENT{Eliminate the $i$th row and column}
\STATE $\Fv \gets \texttt{cholupdate}(\Fv, \tau^{1/2}\vb, '-')$ \quad \COMMENT{$\Oc(L^2)$}
\ENDIF
\STATE $\Gv \gets \Hv(:\,,\texttt{find}(\qb))$
\ENDIF \quad \COMMENT{No change of $(q_i,\Fv)$ otherwise}
\ENDFOR
\STATE \COMMENT{Step~1(b): sampling of \xb}
\STATE $\xb \gets \zerob$ 
\STATE $L \gets \text{sum}(\qb)$
\STATE $\xb(\texttt{find}(\qb)) \gets \Fv\T\left(\sigma_{\epsilon}^{-2}\Fv\Gv\T\zb + \texttt{randn}(L,1)\right)$ \quad \COMMENT{see Eq.~\eqref{eq:UpdateX}}
\end{algorithmic}
\efmi
\end{table}
%------------------------------------------------end algorithm--------------------------------------------------------------------------- 

\section{Simulation results}

\label{sec:Simulation results}

%We dedicate the first subsection to the introduction of a convergence criterion, based on which, simulation tests are conducted in the following to compare the convergence evolution of the above-mentioned MCMC schemes for the problem of BG deconvolution.
%maybe a Bernoulli-laplace model should be tested also
\subsection{Convergence diagnostic}
In order to compare empirically the convergence speed of the different samplers, we have resorted to Brooks and Gelman's iterated graphical method to assess convergence ~\cite{Brooks98}. This diagnostic method is based upon the covariance estimation of $m$ independent Markov chains $\{\Phi_{jt}, j=1,\ldots,m;t=1,\ldots,n\}$ of equal length $n$. Let $\overline{\Phi}_{j.}$ (respectively, $\overline{\Phi}_{..}$) denote the local (respectively, global) mean of the chains. The intra-chain and inter-chain variances are defined as covariance matrix averages:
\balx%\label{eqn:variance}
\Vv_{\text{intra}} &= \frac1{m(n-1)}\sum_{j=1}^m\sum_{t=1}^n(\Phi_{jt}-\overline{\Phi}_{j.})(\Phi_{jt}-\overline{\Phi}_{j.})\T \nonumber\\
\Vv_{\text{inter}} &= \frac1{m-1}\sum_{j=1}^m(\overline{\Phi}_{j.}-\overline{\Phi}_{..})(\overline{\Phi}_{j.}-\overline{\Phi}_{..})\T
\ealx
that characterize the convergence behavior.
%In particularl, the example in Figure~\ref{fig:MCMCTimeshifted} shows signs of local convergence after burn-in (decrease in the spectral norm of $\|\Vv_{\text{intra}}\|_2$), but not of global convergence (since $\|\Vv_{\text{inter}}\|_2$ does not decreases towards zero), thus causing inconsistent solutions.
Brooks and Gelman~\cite{Brooks98} proposed to evaluate the \emph{multivariate potential scale reduction factor} (MPSRF):
$$
\text{MPSRF}=\frac{n-1}{n}+\frac{m+1}{m}\lambda(\Vv_{\text{intra}}\M\Vv_{\text{inter}}) 
$$
where $\lambda(\cdot)$ returns the largest eigenvalue of the covariance matrix. Convergence is diagnosed when MPSRF is close to one (\eg $\text{MPSRF}<1.2$ as proposed in~\cite{Brooks98}). In order to get a graphical evolution of the convergence, each chain is divided into batches of $b$ samples, and the MPSRF is calculated upon the second halves of the Markov chains $\{\Phi_{jt}\}, t=1,\ldots, kb$ of increasing lengths $kb$, while the first halves are regarded as a burn-in period. From an empirical basis, it is also advised to select $b\approx n/20$ in~\cite{Brooks98}.  

\subsection{Simulation tests}
A test scenario is designed here to compare the generalized K-tuple Gibbs sampler (including the classical hybrid sampler in the case $K=1$) of Section~\ref{sec:Generalized Gibbs on K-tuple variables} and the partially marginal sampler of Section~\ref{sec:Partially marginal Gibbs sampler}, in terms of robustness with respect to different random initial conditions. Time-shift and scale ambiguities are taken into account in all the methods, as described in Section~\ref{subsec:Classical Gibbs Sampling}.
Let us reexamine the example of the benchmark sequence (Mendel's sequence, for which $M=300$ and $\lambda=0.1$.) in Figure~\ref{fig:MCMCTimeshifted}.
To concentrate on the convergence quality of the Bernoulli sequence, we evaluate the MPSRF evolutions of $\{\qb^{(k)}\}$ for ten independent Markov chains, \ie $m=10$. 
%we adopt $\Phi_{jt}=\qb_{jt}$. Ten independent Markov chains are considered in parallel (\ie $m=10$).
Figure~\ref{fig:12dBCompare} is compares the convergence rate of $K$-tuple samplers for different values of $K$: while MPSRF falls under $1.2$ after $515$ seconds of simulation for the classical hybrid sampler ($K=1$), this threshold is reached after about $285$ seconds for the 3-tuple sampler and in less than 202 seconds for the 2-tuple sampler.
 %(as also reported in Table~\ref{tab:K1-4}).
It is interesting to observe that both the $3$ and $4$-tuple sampler have a more stable convergence behavior in the convergence zone while the hybrid sampler presents the most oscillating MPSRF values. These results confirm the discussion in Section~\ref{sec:Generalized Gibbs on K-tuple variables}: augmenting $K$ improves the convergence quality in terms of iteration numbers %(also shown in Table~\ref{tab:K1-4})
while in the meantime increases the computational load per iteration. The overall performance relies on the trade-off between the two criteria. Thus, the best compromise in the $K$-tuple sampler {\it series} is reached for $K=2$ in the given example.  
%-----------------------------------------Test of convergence--------------------------------------------------
\begin{figure}[htb]
\mongraphe{\includegraphics[width=160mm]{Variance}}%
{\hspace{60mm}
	{\footnotesize time in seconds (in log-scale)}
}%
{\hspace{22mm}
	{\footnotesize MPSRF (in log-scale)}
}%
{\hspace{90mm}\vspace{50mm}%
	\renewcommand{\arraystretch}{0.8}
	\begin{tabular}{@{}l@{~}l@{}}
		{\scriptsize $\diamond$}      &{\scriptsize K-tuple sampler, $K=1$ (\ie hybrid sampler)}\\
		{\scriptsize $+$}             &{\scriptsize K-tuple sampler, $K=2$}\\
		{\scriptsize $\circ$}         &{\scriptsize K-tuple sampler, $K=3$} \\
		{\scriptsize $\ast$}          &{\scriptsize K-tuple sampler, $K=4$}
	\end{tabular}
}
\caption{Evolution of MPSRF for the $K$-tuple sampler family ($K=1,2,3$) in the case of Mendel's sequence ($\mathrm{SNR}=12.80$\,dB).\label{fig:12dBCompare}}
\end{figure}

Next we compare the convergence diagnostic results on the $2$-tuple sampler and the partially marginal sampler. Each MPSRF in Figure~\ref{fig:GibbsTypeVS2Tuple} is evaluated for every $b=100$ iterations. Altogether 300 iterations (116 seconds of simulation) are required for the partially marginalized sampler to reach convergence, not to mention the overall lower level of MPSRF in the convergence zone compared with the $2$-tuple sampler.   
%Equally shown in Figure~\ref{fig:GibbsTypeVS2Tuple} is the computation load per iteration for each scheme since MPSRF is evaluated for every 100 samples. 
Thus the partially marginalized sampler takes fewer iterations and less time to converge than the best selected $K$-tuple sampler in Figure~\ref{fig:12dBCompare}. 
%the evolution of the computation load per iteration for both methods are easily deduced. 
We also point out that the asymptotic complexity per iteration of the partially marginal sampler is lower than that of the $2$-tuple sampler in the given simulation example, though the contrary is observed for the first 100 samples (corresponding to the heating period) in which the number of detected spikes $L$ is relatively important. 
%----------------------------------------GibbsTypeVS2Tuple----------------------------
\begin{figure}[htb]
%  \centering
%  \renewcommand{\arraystretch}{0.8}
%\btabu{@{}c@{}}
%  \begin{sideways}\parbox{30mm}{\footnotesize MPSRF (in $\log$ scale)}\end{sideways}
%\figc[width=.8\linewidth]{GibbsTypeVS2Tuple} \\
%{\footnotesize time in seconds}	
%\etabu
%\\[-5.5cm]
%%\footnotesize
%\hfill
%\raisebox{-1mm}{\mbox{
%\renewcommand{\arraystretch}{0.8}
%\btabu{@{}c@{~}l@{~~~~~~~~~~~~~~~~~}}
%{\scriptsize $\square$}	   &{\scriptsize partially marginal sampler}\\
%{\scriptsize $+$} &{\scriptsize K-tuple ($K=2$) sampler}
%\etabu}}
%\\[4.5cm]
\mongraphe{\includegraphics[width=160mm]{GibbsTypeVS2Tuple}}%
{\hspace{70mm}
	{\footnotesize time in seconds}
}%
{\hspace{22mm}
	{\footnotesize MPSRF (in log-scale)}
}%
{\hspace{110mm}\vspace{55mm}%
	\renewcommand{\arraystretch}{0.8}
	\begin{tabular}{@{}l@{~}l@{}}
		{\scriptsize $\square$}      &{\scriptsize partially marginal sampler}\\
		{\scriptsize $+$}             &{\scriptsize K-tuple ($K=2$) sampler}
	\end{tabular}
}
\caption{Evolution of MPSRF for the partially marginal sampler and the $2$-tuple sampler in the case of Mendel's sequence. MPSRF has been computed every 100 samples for both methods.  \label{fig:GibbsTypeVS2Tuple}}
\end{figure}

%{\setlength{\baselineskip}{.8\baselineskip}
To further illustrate its robustness, deconvolution results obtained by the partially marginal sampler are shown in Figure~\ref{fig:MCMCSLMR} under identical initial conditions as in Figure~\ref{fig:MCMCTimeshifted}. Only one of the results is reported, as they are undistinguishable from each other (perfectly consistent estimations). And this is true for all the samplers once their MPSRF values drop below the $1.2$ threshold. 

%We noticed also that instead of calculating the posterior Gaussian mean based on $\{\widehat{\qb},\widehat{\hb},\zb\}$, the estimator $\widehat{\xb}$ can also be obtained by averaging the auxiliarily sampled \xb in step~\stepdeux(a). In fact, as $\{\widetilde{\Theta}^{(k)}\}$ is shown to converge to the distribution $\rho(\widetilde{\Theta}\I\zb)$ in Section~\ref{sect:NewGibbs}, step~\stepdeux(a) samples \xb conditional on $\widetilde{\Theta}$ and \zb such that taken marginally $\xb^{(k)}$ follows the distribution $\rho(\xb\I\zb)$.  

%---------------------------------- Gibbs type results----------------------------------
\begin{figure}[htb]
  \centering
\setlength{\tabcolsep}{0pt}
\btabu{cc}
   \figc[height=5cm]{SMLRX20}%,height=3cm}
&  \figc[height=5cm]{SMLRIR20} \\
		{\footnotesize estimation of \xb}
&  	{\footnotesize estimation of \hb}
\etabu
\caption{Estimation result using the proposed sampler, under the same conditions as in Figure~\ref{fig:MCMCTimeshifted}. The 10 independent Markov chains yield undistinguishable estimates.}
\label{fig:MCMCSLMR}
\end{figure}
%-----------------------------------------------------------------------------------------
In the given example of the Mendel's sequence (for which $M=300$), the partially marginal sampler achieves a better compromise between the two criteria~\cite{Liu01} than the most performant sampler in the generalized $K$-tuple family: the time required to reach convergence is reduced by half.
%This however dose not mean that the $K$-tuple sampling algorithm is completely of no interest compared with the partially marginal sampler. 
%The issue of complexity related to $L$ is further discussed in the following example, in which we show that the partially marginal sampler do not necessarily outperform the $K$-tuple sampler in dealing with longer datasets.  
However, the partially marginal sampler do not necessarily outperform the $K$-tuple sampler in dealing with longer datasets.
Notice that the computational complexity for a Cholesky update operation in the latter case is $\Oc(L^2)$ ($L=\sum q_i\approx \lambda M$, the number of detected spikes), thus proportional to the data length of the BG sequence, whereas the complexity per iteration of the $K$-tuple sampler is non-increasing with respect to the number of detected spikes $L$. To illustrate the disadvantage of the partially marginal sampler, it is necessary to test on simulation examples of longer observation data \zb while fixing the Bernoulli parameter $\lambda$. Figure~\ref{fig:CostIter} compares the complexity per iteration as a function of BG sequence length $M$ for the two algorithms. The $2$-tuple sampler potentially outperforms the partially marginal sampler due to its linear growth of complexity per iteration with respect to $M$. 
%Actually, measured in total convergence time in the example of $M=3000$, the $2$-tuple sampler is about $6$ times faster than the partially marginal method. We have reason to believe that this factor will further increase by augmenting the dataset length. 
Figure~\ref{fig:ConvergTime} resumes the convergence time for all the $5$ samplers ($K=1,\ldots, 4$ in solide lines and the partially marginalised sampler in dashed line) on simulation data whose lengths vary from $100$ to $1600$ and $\lambda = 10\%$. While the partially marginal sampler outperforms all the rest algorithms in areas of shorter data length ($M<800$), its quadratic computation cost per iteration as shown in Fig.~\ref{fig:CostIter} yields the contrary in applications with longer data lengths. We also note that for the given bernoulli parameter $\lambda = 10\%$, the hybrid sampler is never a good choice inside the $K$-tuple family regardless of the data length. 
  
\begin{figure}[htb]
\mongraphe{\hspace{5mm} \includegraphics[height=90mm]{CostEvolution}}%
{\hspace{65mm} {\footnotesize spike train length $M$}}
{\hspace{30mm} {\footnotesize cost per iteration in second}}
{\hspace{25mm}\vspace{80mm}
	\renewcommand{\arraystretch}{0.8}
	\begin{tabular}{@{}l@{~}l@{}}
		{\scriptsize $+$ }              &{\scriptsize K-tuple sampler, $K=2$}\\
		{\scriptsize $\square$}         &{\scriptsize Partially marginal sampler}
	\end{tabular}
}
\caption{Computation cost per iteration comparison between the $2$-tuple Gibbs sampler and the partially marginal sampler. The two methods are tested on the same simulated observation data \zb, for which the length of the spike train ranges from $100$ to $3000$ and the Bernoulli parameter $\lambda$ is fixed to $10\%$. Polynomial interpolations (of degree $1$ and $2$ respectively) are traced to show the quasi-linear and -quadratic evolution of complexity per iteration with respect to the spike train length $M$. }
\label{fig:CostIter}
\end{figure}

\begin{figure}[htb]
\mongraphe{\hspace{5mm} \includegraphics[height=100mm]{ConvergTime}}%
{\hspace{60mm} {\footnotesize spike train length $M$ (in log-scale)}}
{\hspace{40mm}{\footnotesize time in second}}
{\hspace{20mm}\vspace{80mm}
	\renewcommand{\arraystretch}{0.8}
	\begin{tabular}{@{}l@{~}l@{}}
		{\scriptsize $\diamond$}        &{\scriptsize K-tuple sampler, $K=1$ (\ie hybrid sampler)}\\
		{\scriptsize $+$ }              &{\scriptsize K-tuple sampler, $K=2$}\\
		{\scriptsize $\circ$  }   			&{\scriptsize K-tuple sampler, $K=3$} \\
		{\scriptsize $\ast$}            &{\scriptsize K-tuple sampler, $K=4$} \\
		{\scriptsize $\square$}         &{\scriptsize Partially marginal sampler}
	\end{tabular}
}
\caption{Convergence time comparison for the $5$ samplers applied on spike trains of increasing length. \label{fig:ConvergTime}}
\end{figure}

%\section{Conclusions}\label{sec:Conclusions}
%In this paper, we have examined and analysed the Gibbs samplers in dealing with the spike train deconvolution problem. To avoid being stuck in local optimal solutions, a new Gibbs-type sampler is proposed and proved to converge to the marginalized posterior distribution. Special techniques are introduced to keep the invariant condition and to evaluate the marginal posterior law for each {\it Bernoulli} variable $q_i$ at reasonable computation costs. Simulation results confirm the assumed robustness and efficiency of the proposed sampler.

%Our main perspective in the future is to generalize this approach. Reducing the sample dimension can achieve better convergence rate by trading off the simplicity in sampling. The overall gain is considerable by marginalizing variables that form highly-dependant pairs with others in the Markov chain. Another condition of this approach is that marginalized variable should be easily estimated ultimately. 
%In fact, the proposed sampler is equivalent to the classical Gibbs sampler with $\Theta=\{\qb,\xb,\hb,\lambda,\sigma_{\epsilon},\sigma_{\hb}\}$, the only difference being that each $q_i$ is sampled according to $\int\rho(q_i|\Theta\setminus q_i;\zb)\dd\xb$, independent of \xb.

% -------------------------------------------------------------------------
\section{Conclusion}\label{sec:Conclusion}
This paper proposes two distinct methods in dealing with the identified inefficiency concerning the Bernoulli label \qb sampling in the BG deconvolution problem. Detailed algorithms are given for both methods as well as their performance on simulation tests. 
It is shown that both methods are mathematically valid MCMC sampling schemes and both achieve better convergence properties in comparison to the hybrid sampler, the latter is based on Cheng \etal{}'s Gibbs scheme with time-shift compensations and scale re-sampling. Both proposed methods demonstrate better trade-off for the convergence rate by reducing drastically the iteration steps needed to converge. On the simulation test of Mendel's sequence, the partially marginal sampler is shown to out-perform the whole class of K-tuple samplers, among which the optimized parameter is found at $K=2$. This is however due to the relatively small size of the simulation problem and the contrary is observed by simply augmenting the simulated data length. 
 
We conclude with some extensions of the proposed techniques. By combining the idea of a K-tuple joint sampling and that of partial marginalization of \xb, \ie updating $\qb_{i+\Omega} \I \qb_{-(i+\Omega)},\hb,\sigma_{\epsilon},\lambda$ in Step~1 of Table~\ref{Di}, one might further reduce the number of iterations necessary to escape from local optima at the cost of higher computation load per iteration. Such a solution should be adopted when the deconvolution problem presents more difficulties, characterized by either low SNR, band-limited IR and/or less sparse spike trains. 
%From the results in Figure~\ref{fig:12dBCompare} and Figure~\ref{fig:GibbsTypeVS2Tuple}, it is clear that the combination is unnecessary in this particular example. 
%Regarding the burn-in computation burden due to over-detected spike numbers $L$, a trivial solution would be to initialize the Markov chain using a \emph{heated} sample produced by the reference Gibbs scheme and thus restrain $L$ to smaller values. 

%Another interesting issue is concerned with the relaxation of the Gaussian assumption of \xb. Actually, the generalization of the K-tuple sampling is straightforward as long as Cheng \etal{}'s method is applicable; on the contrary, integration of \xb as a non-Gaussian random vector is less evident if not impossible. These might be the limits of the second approach.  
\renewcommand{\theequation}{A-\arabic{equation}}    
\setcounter{equation}{0}  % reset counter
\section*{APPENDIX}
To update the Cholesky factor \Rv on removing a spike at the location $i$ such that $\Kv\T\Kv = \Rv\T\Rv - \vb\vb\T$, we specify : 
\bal
\Rv = \left[
\barr{ccc}
\Av & \bb & \Cv \\
%\hline
0   & \delta & \eb\T \\
%\hline
0   &  0     & \Fv  
\earr
\right], 
\Kv = \left[
\barr{ccc}
\Gv & 0 & \Hv \\
0   & 0 & 0 \\
0   &  0     & \Jv  
\earr
\right],\label{eq:RK}\\
\left\{
\barr{ll}
\vb &= \sqrt{\tau}\Rv\T \cdot \Rv(:,j) \\
\tau\M &= \bb\T \bb + \delta^2
\earr
\right. \label{eq:vtau}
\eal
where $\{\Av, \Fv, \Gv, \Jv\}$ are upper-triangular matrices and $\delta = \Rv_{i,i}$ the $i$th diagonal value of \Rv. $(\Kv_{-i,-i})$ is therefor the updated Cholesky factor. 

A most evident solution consists an operation $\Kv = \text{Cholupdate}(\Rv, \vb, '-')$ followed by an extraction to get the updated Cholesky factor (of lowered dimension) :
\balx
\left[
\barr{cc}
\Gv & \Hv \\
0   & \Jv  
\earr
\right], 
\ealx
which however systematically fails due to the fact that $\Kv\T \Kv$ is no longer definite positive. In fact, one necessary condition~\cite{Seeger08} for a successful rank-$1$ \emph{downdate} operation (update by subtracting $\vb\vb\T$) is that $\rho^2 = 1- \pb\T\pb > 0$ for which \pb is defined by 
\bal
\Rv\T \pb = \vb. 
\eal
Without inverting $\Rv\T$, it can be deduced from ~\eqref{eq:vtau} that $\pb = \sqrt{\tau} \Rv(:,j)$ and thus $\rho^2 = 0$. The algorithm of \emph{Cholupdate} by a series of \emph{Givens} rotation should break down according to~\cite{Seeger08}. 

The problem of non-definite positiveness is avoided by a strategy that first lower the dimension and then perform the \emph{Cholupdate} algorithm. 
From ~\eqref{eq:RK}, we obtain: 
\balx
\Rv\T \Rv = \left[
\barr{c|c|c}
\Av\T\Av & \Av\T\bb & \Av\T\Cv \\ \hline
\bb\T\Av & \tau\M   & \bb\T\Cv + \delta\eb\T \\ \hline
\Cv\T\Av & \Cv\T\bb + \delta\eb   & \Cv\T\Cv + \eb\eb\T + \Fv\T\Fv  
\earr
\right]
= \left[
\barr{c|c|c}
\Gv\T\Gv & 0 & \Gv\T\Hv \\ \hline
0   & 0 & 0 \\ \hline
\Hv\T\Gv   &  0     & \Hv\T\Hv +\Jv\T\Jv   
\earr
\right]+ \vb\vb\T
\ealx
The submatrices on both sides by eliminating the $i$th row and column yields :
\balx
\left[
\barr{c|c}
\Av\T\Av &  \Av\T\Cv \\ \hline
\Cv\T\Av & \Cv\T\Cv + \eb\eb\T + \Fv\T\Fv  
\earr
\right]
= \left[
\barr{c|c}
\Gv\T\Gv  & \Gv\T\Hv \\  \hline
\Hv\T\Gv  & \Hv\T\Hv +\Jv\T\Jv   
\earr
\right]+ \vb_{-i}\vb_{-i}\T
\ealx
for which $\vb_{-i}$ denotes \vb without the $i$th component. 
The following relation is key to the algorithm that we propose : 
\bal
\left[
\barr{cc}
\Av &  \Cv \\
0   &  \Fv
\earr
\right]\T 
\left[
\barr{cc}
\Av &  \Cv \\
0   &  \Fv
\earr
\right] + \left[
\barr{c}
0 \\ \eb
\earr
\right] \left[
\barr{c}
0 \\ \eb
\earr
\right]\T
= \left[
\barr{cc}
\Gv  & \Hv \\
0    & \Jv   
\earr
\right]\T \left[
\barr{cc}
\Gv  & \Hv \\
0    & \Jv   
\earr
\right] + \vb_{-i}\vb_{-i}\T
\eal
The two steps that leads to a spike removing update are : 
\ben[align=left, leftmargin=*, noitemsep]
\item an update on \Fv (of reduced dimension): 
\balx
\Mv = \text{Cholupdate}(\Fv, \eb, '+');  
\ealx
\item a downdate that keeps the positiveness of the product matrix and will not fail: 
\balx
\left[
\barr{cc}
\Gv  & \Hv \\
0    & \Jv   
\earr
\right] = \text{Cholupdate}(\left[
\barr{cc}
\Av &  \Cv \\
0   &  \Mv
\earr
\right], \vb_{-i}, '-'). 
\ealx
it can be proven that $\rho^2 = 1-\pb\T\pb$ is strictly positive in this case. 
\een
 
%\bibliographystyle{IEEEbib}
\bibliographystyle{ieeeji}
%\bibliography{bibenabr,revuedef,revueabr,baseAJ,baseKZ,gpipubli}
\begin{thebibliography}{10}

\bibitem{Cheng96}
Q.~Cheng, R.~Chen, and T.-H. Li,
\newblock ``Simultaneous wavelet estimation and deconvolution of reflection
  seismic signals'',
\newblock {\em {{IEEE} {T}rans. {G}eosci. {R}emote {S}ensing}}, vol. 34, pp.
  377--384, {M}ar. 1996.

\bibitem{Kormylo83a}
J.~J. Kormylo and J.~M. Mendel,
\newblock ``Maximum-likelihood seismic deconvolution'',
\newblock {\em {{IEEE} {T}rans. {G}eosci. {R}emote {S}ensing}}, vol. GE-21, no.
  1, pp. 72--82, {J}an. 1983.

\bibitem{Bourguignon05}
S.~Bourguignon and H.~Carfantan,
\newblock ``Bernoulli-{G}aussian spectral analysis of unevenly spaced
  astrophysical data'',
\newblock in {\em IEEE Workshop Stat. Sig. Proc.}, Bordeaux, France, {J}uly
  2005, pp. 811--816.

\bibitem{Champagnat96a}
F.~Champagnat, Y.~Goussard, and J.~Idier,
\newblock ``Unsupervised deconvolution of sparse spike trains using stochastic
  approximation'',
\newblock {\em {{IEEE} {T}rans. {S}ignal {P}rocessing}}, vol. 44, no. 12, pp.
  2988--2998, {D}ec. 1996.

\bibitem{Liu01}
J.~S. Liu,
\newblock {\em Monte {C}arlo Strategies in Scientific Computing},
\newblock Springer Series in Statistics. Springer Verlag, New York, NY, 2001.

\bibitem{Robert04}
C.~P. Robert and G.~Casella,
\newblock {\em Monte {C}arlo Statistical Methods},
\newblock Springer Texts in Statistics. Springer Verlag, New York, NY, 2nd
  edition, 2004.

\bibitem{Labat06b}
C.~Labat and J.~Idier,
\newblock ``Sparse blind deconvolution accounting for time-shift ambiguity'',
\newblock in {\em {{P}roc. {IEEE} {ICASSP}}}, Toulouse, France, {M}ay 2006, pp.
  616--619.

\bibitem{Veit08}
T.~Veit, J.~Idier, and S.~Moussaoui,
\newblock ``{R}��chantillonnage de l'�chelle dans les algorithmes {MCMC} pour
  les probl�mes inverses bilin�aires'',
\newblock {\em {{T}raitement du {S}ignal}}, vol. 25, no. 4, pp. 329--343, 2008.

\bibitem{Chi84}
C.~Y. Chi and J.~M. Mendel,
\newblock ``Improved maximum-likelihood detection and estimation of
  {B}ernoulli-{G}aussian processes'',
\newblock {\em {{IEEE} {T}rans. {I}nf. {T}heory}}, vol. 30, pp. 429--435, 1984.

\bibitem{Ge08}
D.~Ge, J.~Idier, and E.~{Le Carpentier},
\newblock ``A new {MCMC} algorithm for blind {B}ernoulli-{G}aussian
  deconvolution'',
\newblock in {\em EUSIPCO}, Lausanne, Switzerland, {S}ep. 2008.

\bibitem{Dobigeon07}
N.~Dobigeon, J.-Y. Tourneret, and J.~D. Scargle,
\newblock ``Joint segmentation of multivariate astronomical times series:
  {B}ayesian sampling with a hierarchical model'',
\newblock {\em {{IEEE} {T}rans. {S}ignal {P}rocessing}}, vol. 55, no. 2, pp.
  414--423, {F}eb. 2007.

\bibitem{Brooks98}
S.~P. Brooks and A.~Gelman,
\newblock ``General methods for monitoring convergence of iterative
  simulations'',
\newblock {\em {{J}. {C}omput. {G}raph. {S}tatist.}}, vol. 7, no. 4, pp.
  434--455, 1998.

\bibitem{Liu94b}
J.~S. Liu, W.~H. Wong, and A.~Kong,
\newblock ``Covariance structure of the {G}ibbs sampler with applications to
  the comparisons of estimators and augmentation schemes'',
\newblock {\em {{B}iometrika}}, vol. 81, pp. 27--40, 1994.

\bibitem{Cappe99}
O.~Capp�, A.~Doucet, M.~Lavielle, and E.~Moulines,
\newblock ``Simulation-based methods for blind maximum-likelihood filter
  identification'',
\newblock {\em {{S}ignal {P}rocessing}}, vol. 73, no. 1-2, pp. 3--25, {F}eb.
  1999.

\bibitem{Park09}
T.~Park and D.~A. van Dyk,
\newblock ``{Partially Collapsed Gibbs Samplers: Illustrations and
  Applications}'',
\newblock {\em Journal of Computational and Graphical Statistics}, vol. 18, no.
  2, pp. 283--305, 2009.

\bibitem{Seeger08}
M.~Seeger,
\newblock ``{Low Rank Updates for the Cholesky Decomposition}'',
\newblock Tech. {R}ep., Department of EECS, University of California at
  Berkeley, 485 Soda Hall, Berkeley CA, 2008.

\end{thebibliography}
\end{document}

sigproc09.bbl
\begin{thebibliography}{}

\end{thebibliography}
