blind source separation ( bss ) aims to recover unknown source signals from observed mixtures with or without very limited information about their mixing process . bss problems have been addressed in many previous studies , for example , [ 112 ] , which were motivated by several real - world applications . in a cocktail - party problem , microphones receive noisy mixtures of acoustic signals that propagate along multiple paths from their sources . in a real scenario , the number of audio sources may be greater than the number of microphones , audio sources may have different timbres and similar pitches , and audio signals may be only locally stationary . a convolutive and under - determined mixing approach needs to be adopted to model this problem . some of these   operate in the time - domain by solving the alternative finite impulse response ( fir ) inverse model using independent component analysis ( ica ) methods . another method is to extract meaningful features from the time - frequency ( tf ) representations of mixtures . this approach seems to be more efficient than the ica - based techniques especially when the number of microphones is lower than the number of sources . acoustic signals are usually sparse in the tf domain , so the source signals can be separated efficiently even if they are partially overlapped and the problem is under - determined . these features can be extracted using several techniques , including tf masking [ 15 , 16 ] , frequency bin - wise clustering with permutation alignment ( fbwc - pa ) [ 17 , 18 ] , subspace projection , hidden markov models ( hmm ) , interaural phase difference ( ipd ) , nonnegative matrix factorization ( nmf ) [ 22 , 23 ] , and nonnegative tensor factorization ( ntf ) . nonnegative matrix factorization   is a feature extraction method with many real - world applications . ozerov and fevotte   developed the em - nmf algorithm , which is suitable for unsupervised convolutive and possibly under - determined unmixing of audio sources using only stereo observations . their model of the sources was based on the generalized wiener filtering model [ 2729 ] , which assumes that each source is locally stationary and that it can be expressed in terms of superimposed amplitude - modulated gaussian components . thus , a power spectrogram of each source can be factorized into lower - rank nonnegative matrices , which facilitates the use of nmf for estimating the frequency and temporal profiles of each latent source component . in the tf representation , , this assumption is very weak for any adjacent bins because real audio signals have locally smooth frequency and temporal structures . motivated by several papers on smoothness [ 26 , 28 , 30 , 31 ] in bss models , we attempt to further improve the em - nmf algorithm by enforcing local smoothness both in the frequency and temporal profiles of the nmf factors . similar to [ 28 , 30 , 32 ] , we introduce a priori knowledge to the nmf - based model using a bayesian framework , although our approach is based on a gibbs prior with a markov random field ( mrf ) model to describe pairwise interactions among adjacent bins in spectrograms . as demonstrated in , the mrf model with green s function , which is well known in many tomographic image reconstruction applications , can improve the em - nmf algorithm . in this paper , we extend the results presented in   using other smoothing functions , particularly a more flexible simultaneous autoregressive ( sar ) model that is more appropriate in term of hyperparameter estimation and computational complexity . let i microphones receive signals that can be modeled as a noisy convolutive mixture of j audio signals . the signal received by the i - th microphone ( \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$i=1 , \ldots , i$$\end{document } ) can be expressed as1\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \tilde x_i(t)=\sum_{j=1}^j \sum_{l=0}^{l-1 } \tilde a_{ijl } \tilde s_j(t - l)+\tilde n_i(t ) , $ $ \end{document}where \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\tilde a_{ijl}$$\end{document } represents the corresponding mixing filter coefficient , \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\tilde s_j(t)$$\end{document } is the j - th source signal ( \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$j=1 , \ldots , j$$\end{document } ) , \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\tilde n_i(t)$$\end{document } is the additive noise , and l is the length of the mixing filter . in the tf domain , the model ( 1 ) can be expressed as2\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ x_{ift}=\sum_{j=1}^j a_{ijf } s_{jft}+n_{ift } , \quad \hbox{or equivalently } , \quad \varvec{x}_f=\varvec{a}_f \varvec{s}_f + \varvec{n}_f , $ $ \end{document}where \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$${\varvec{x}_f=[x_{ift}]_f \in { \mathbb c}^{i \times t } , \varvec{a}_f=[a_{ijf}]_f \in { \mathbb c}^{i \times j } , \varvec{s}_f=[s_{jft}]_f \in { \mathbb c}^{j \times t } , \varvec{n}_f=[n_{ift}]_f \in { \mathbb c}^{i \times t}}$$\end{document } , and \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$f=1 , \ldots , f$$\end{document } is the index of a frequency bin . the noise nift is assumed to be stationary and spatially uncorrelated , i.e,3\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ n_{ift } \sim { \mathcal{n}}_c(0 , \varvec{\upsigma}_n ) , $ $ \end{document}where \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\varvec{\upsigma}_n=\hbox{diag}\left ( [ \sigma^2_{i } ] \right)$$\end{document } and \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\mathcal{n}_c(0 , \varvec{\upsigma}_n)$$\end{document } is a proper complex gaussian distribution with a zero - mean and the covariance matrix \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\varvec{\upsigma}_n$$\end{document}. benaroya et  al . described an audio source \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\tilde{s}(t)$$\end{document } as a superimposed amplitude - modulated gaussian process:4\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \tilde{s}(t)=\sum_{r=1}^r h_r(t ) \tilde{w}_r(t ) , $ $ \end{document}where hr(t ) is a slowly varying amplitude parameter in the r - th component ( \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$r=1 , \ldots , r$$\end{document } ) , and \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\tilde{w}_r(t)$$\end{document } is a stationary zero - mean gaussian process with the power spectral density r2 ( f ) . the tf representation of ( 4 ) leads to5\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ s(f , t ) \sim { \mathcal{n}}_c \left ( 0 , \sum_{r=1}^r h_r(t ) \sigma^2_{r}(f ) \right ) . $ $ \end{document}the power spectrogram of ( 5 ) is given by |sft|  =  r=1rwfrhrt , where wfr thus , the spectrogram of the j - th source \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\tilde{s}_j(t)$$\end{document } can be factorized as follows:6\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ |\varvec{s}_j|^2=\varvec{w}_j \varvec{h}_j , $ $ \end{document}where \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$${\varvec{s}_j \in { \mathbb c}^{f \times t } , \varvec{w}_j \in { \mathbb r}_+^{f \times r_j } , \varvec{h}_j \in { \mathbb r}_+^{r_j \times t } , r_j}$$\end{document } is the number of latent components in the j - th source , and \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$${\mathbb r}_+$$\end{document } is the nonnegative orthant of the euclidean space . the column vectors of \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\varvec{w}_j$$\end{document } represent the frequency profiles of the j - th source , while the row vectors of \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\varvec{h}_j$$\end{document } are the temporal profiles . .   transformed the model ( 5 ) to the following form:7\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ s(f , t)=\sum_{r=1}^r c_r ( f , t ) $ $ \end{document}where \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$c_r(f , t ) \sim \mathcal{n}_c \left ( 0 , h_r(t ) \sigma^2_{r}(f ) \right)=\mathcal{n}_c \left ( 0 , \varvec{\upsigma}_{ft}^{(c ) } \right)$$\end{document}. thus,8\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \varvec{\upsigma}_{ft}^{(c)}=\hbox{diag}\left ( [ w_{fr } h_{rt}]_r \right)=\hbox{diag}\left ( \left [ |c_{rft}|^2 \right]_r \right ) , $ $ \end{document}and \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\left [ \sum\nolimits_{r=1}^{r } |c_{rft}|^2 \right]=\varvec{w}\varvec{h}$$\end{document } , where \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$${\varvec{w } \in { \mathbb r}^{f \times r}_+ , \varvec{h } \in { \mathbb r}^{r \times t}_+}$$\end{document}. consequently , the model ( 2 ) can be expressed as9\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ x_{ift}=\sum_{r=1}^{r } \bar{a}_{irf } c_{rft}+n_{ift } , \quad\hbox { where } \quad s_{jft}=\sum_{r \in { \mathfrak{r}}_j } c_{rft } , $ $ \end{document}\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$${r=|\mathfrak{r}|}$$\end{document } is the number of entries in the set \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$${{\mathfrak{r}}=\bigcup\nolimits_{j=1}^j \mathfrak{r}_j}$$\end{document } , and \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$${\bar{\varvec{a}}_f=[\bar{a}_{irf } ] \in { \mathbb c}^{i \times r}}$$\end{document } is created from the columns of the matrix \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\varvec{a}_f$$\end{document}. for example , assuming \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$${\forall j : \mathfrak{r}_j=\{1 , \ldots , \bar{r}\}}$$\end{document } , we have \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$r = j \bar{r}$$\end{document } , and \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$${\bar{\varvec{a}}_f=[\varvec{a}_f , \ldots , \varvec{a}_f ] \in { \mathbb c}^{i \times r}}$$\end{document } is the augmented mixing matrix   created from \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\bar{r}$$\end{document } matrices \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\varvec{a}_f$$\end{document}. from ( 8) and ( 9 ) , we have \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$s_{jft } \sim \mathcal{n}_c(0 , \varvec{\upsigma}_{ft}^{(s)})$$\end{document } where \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$${\varvec{\upsigma}_{ft}^{(s)}=\hbox{diag}\left ( \left [ \sum\nolimits_{r \in \mathfrak{r}_j } w_{fr } h_{rt } \right]_j \right)}$$\end{document}. to estimate the parameters \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$${\mathcal{a}=[a_{ijf } ] \in { \mathbb c}^{i \times j \times f } , \mathcal{c}=[c_{rft } ] \in { \mathbb c}^{r \times f \times t } , \varvec{w } \in { \mathbb r}_+^{f \times r } , \varvec{h } \in { \mathbb r}_+^{r \times t}}$$\end{document } , and \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$${\varvec{\upsigma}_n \in { \mathbb r}_+^{i \times i}}$$\end{document } , we formulate the following posterior:10\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ p({\mathcal{c}},\varvec{w},\varvec{h}|{\mathcal{x}},{\mathcal{a}},\varvec{\upsigma}_n)= \frac{p({\mathcal{x}}|{\mathcal{c } } , { \mathcal{a}},\varvec{\upsigma}_n ) p({\mathcal{c}}|\varvec{w},\varvec{h } ) p(\varvec{w } ) p(\varvec{h})}{p({\mathcal{x}}|{\mathcal{a } } , \varvec{\upsigma}_n ) } , $ $ \end{document}from which we obtain11\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \begin{aligned } \ln p({\mathcal{x}},{\mathcal{c}},\varvec{w},\varvec{h}|{\mathcal{a}},\varvec{\upsigma}_n)&= . \end{aligned } $ $ \end{document } from ( 3 ) and ( 9 ) , we have the joint conditional pdf for \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\mathcal{x}$$\end{document}:12\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ p({\mathcal{x}}|{\mathcal{c}},{\mathcal{a}},\varvec{\upsigma}_n)=\prod_{i , f , t } { \mathcal{n}}_c \left ( \sum_{r=1}^{r } \bar{a}_{irf } c_{rft } , \sigma^2_{i } \right)=\prod_{f , t } { \mathcal{n}}_c(\varvec{a}_f \varvec{s}_{ft } , \varvec{\upsigma}_n ) . $ $ \end{document } based on ( 12 ) , the log - likelihood term in ( 11 ) can be expressed as13\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \begin{aligned } \ln p({\mathcal{x}}|{\mathcal{c}},{\mathcal{a}},\varvec{\upsigma}_n)&= - \sum_{f , t } ( \varvec{x}_{ft } - \bar{\varvec{a}}_f \varvec{c}_{ft})^h \varvec{\upsigma}_n^{-1}(\varvec{x}_{ft } - \bar{\varvec{a}}_f \varvec{c}_{ft } ) - \sum_{f , t } \ln \det \varvec{\upsigma}_n \\ & = - \sum_{f , t } ( \varvec{x}_{ft } - \varvec{a}_f \varvec{s}_{ft})^h \varvec{\upsigma}_n^{-1}(\varvec{x}_{ft } - \varvec{a}_f \varvec{s}_{ft } ) \\ & \quad-\sum_{f , t } \ln \det \varvec{\upsigma}_n+\hbox{const},\\ \end{aligned } $ $ \end{document}where \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$${\varvec{c}_{ft}=[c_{1 ft } , \ldots , c_{rft}]^t \in { \mathbb c}^r , \varvec{s}_{ft}=[s_{1 ft } , \ldots , s_{jft}]^t \in { \mathbb c}^j}$$\end{document } , and \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$${\varvec{x}_{ft}=[x_{1 ft } , \ldots , x_{ift}]^t \in { \mathbb c}^i}$$\end{document}. the joint conditional pdf for \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\mathcal{c}$$\end{document } comes from the model ( 5):14\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \begin{aligned } p({\mathcal{c}}|\varvec{w},\varvec{h})&= \prod_{r=1}^r \prod_{f=1}^f \prod_{t=1}^t { \mathcal{n}}_c \left ( 0 , w_{fr } h_{rt } \right ) \\ & = \prod_{r=1}^r \prod_{f=1}^f \prod_{t=1}^t |\pi w_{fr } h_{rt}|^{-1 } \exp \left \{- \frac{|c_{rft}|^2 } { w_{fr } h_{rt } } \right \}. \end{aligned } $ $ \end{document}from ( 14 ) , we have the log - likelihood functional for \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\mathcal{c } : $ $ \end{document}15\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \ln p({\mathcal{c}}|\varvec{w},\varvec{h})=- \sum_{r , f , t } \left ( \ln(w_{fr}h_{rt})+\frac{|c_{rft}|^2}{w_{fr}h_{rt } } \right ) + \hbox{const}. $ $ \end{document}the negative log - likelihood in ( 15 ) is the itakura - saito ( is ) divergence , which is particularly useful for measuring the goodness of fit between spectrograms . the is divergence is the special case of the -divergence when \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\beta \rightarrow -1$$\end{document } . the priors \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$p(\varvec{w})$$\end{document } and \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$p(\varvec{h})$$\end{document } in ( 10 ) can be determined in many ways . proposed the determination of priors using markov chains and the inverse gamma distribution . in our approach , we propose to model the priors with the gibbs distribution , which is particularly useful for enforcing local smoothness in images . let us assume that prior information on the total smoothness of the estimated components w and h is modeled using the following gibbs distributions:16\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ p(\varvec{w})=\frac{1}{z_w } \exp \left \{-\alpha_w u(\varvec{w } ) \right \ } , \quad p(\varvec{h})=\frac{1}{z_h } \exp \left \{-\alpha_h u(\varvec{h } ) \right \ } $ $ \end{document}where zw and zh are partition functions , w and h are regularization parameters , and u(p ) is a total energy function , which measures the total roughness in p. the function u(p ) is often formulated with respect to the mrf model , which is commonly used in image reconstruction for modeling local smoothness . the functions u(w ) and u(h ) can be determined for the matrices w and h in the following way:17\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ u(\varvec{w})= \sum_{f , r } \sum_{l \in s_f } \nu_{fl } \psi \left(w_{fr } - w_{lr } , \delta_w \right ) , $ $ \end{document}18\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ u(\varvec{h})= \sum_{t , r } \sum_{l \in s_t } \nu_{tl } \psi \left(h_{rt } - h_{rl } , \delta_h \right ) . $ $ \end{document}in the first - order interactions ( nearest neighborhood ) , we have sf  =  { f    1 , f  +  1 } and the weighting factor fl  =  1 , and st  =  { t    1 , t  +  1 } with tl  =  1 . in the second - order interactions , sf  =  { f    2 ,  f    1 ,  f  +  1 ,  f  +  2 } and st  =  { t    2 ,  t    1 ,  t  +  1 ,  t  +  2}. the parameters w and h are scaling factors , while \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\psi \left(\xi , \delta \right)$$\end{document } is a potential function of  that can take different forms . the potential functions that can be applied to the em - nmf algorithm are listed in table  1.table  1potential functionsauthor(s ) ( name)functions : ( ,  )reference(gaussian)(/ ) 
 besag ( laplacian ) 
 \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\left|\xi/\delta \right|$$\end{document } 
 bouman and sauer ( ggmrf)|/| 
 geman and mcclure 
 \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\frac{16}{3 \sqrt{3}}\frac{(\xi/\delta)^2}{(1+(\xi/\delta)^2 ) } $ $ \end{document } 
 geman and reynolds 
 \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\frac{|\xi/\delta|}{1+|\xi/\delta|}$$\end{document } 
 green 
 \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\delta \ln [ \cosh(\xi/\delta)]$$\end{document } 
 hebert and leahy 
 \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\delta \ln [ 1+(\xi/\delta)^2]$$\end{document } 
  according to lange , a robust potential function in the gibbs prior should have the following properties : nonnegative , even , 0  at   =  0 , strictly increasing for   >  0 , unbounded , and convex with bounded first - derivative . of the functions listed in table  1 , green s function satisfies all of these properties , and consequently , it was selected for use in the tests in . unfortunately , the application of green s function to both matrices w and h demands the determination of two hyperparameters w and h , and two penalty parameters w and h . moreover , data - driven hyperparameter estimation usually involves an approximation of the partition functions \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\varvec{z}_w$$\end{document } and \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\varvec{z}_h$$\end{document } , which is not easy in this task . the gaussian function \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\psi \left(\xi , \delta \right)=\left(\xi/\delta \right)^2$$\end{document } , as shown in table  1 , does not have a bounded first - derivative , but its scaling parameter  may be merged with a penalty parameter . consequently , only two parameters need to be determined . the mrf model with a gaussian potential function is actually the sar model [ 4244 ] , which is used widely in many scientific fields [ 44 , 45 ] to represent the interactions among spatial data with gaussian noise . let \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$${\varvec{w}_r \in { \mathbb r}^f_+}$$\end{document } be the r - th column of the matrix w , and \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$${\underline { { \bf h}}_r \in { \mathbb r}^{1 \times t}_+}$$\end{document } be the r - th row of the matrix h. random variables in the vectors \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\varvec{w}_r$$\end{document } and \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\underline{{\bf h}}_r$$\end{document } can be modeled using the following stochastic equations:19\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \varvec{w}_r=\varvec{s}^{(w ) } \varvec{w}_r+\epsilon , \quad \underline { { \bf h}}_r=\underline { { \bf h}}_r \varvec{s}^{(h)}+\epsilon , $ $ \end{document}where \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$${\varvec{s}^{(w ) } \in { \mathbb r}^{f \times f}}$$\end{document } and \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$${\varvec{s}^{(h ) } \in { \mathbb r}^{t \times t}}$$\end{document } are symmetric matrices of spatial dependencies between the random variables , \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\epsilon \sim \mathcal{n}(0 , gaussian noise , and i is an identity matrix with a corresponding size . according to [ 45 , 46 ] , the spatial dependence matrices can be expressed as \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\varvec{s}^{(w)}=\gamma \varvec{z}^{(w)}$$\end{document } and \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\varvec{s}^{(h)}=\gamma \varvec{z}^{(h)}$$\end{document } , where  is a constant that ensures that the matrices \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\varvec{c}^{(w)}={\bf i}_f - \varvec{s}^{(w)}$$\end{document } and \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\varvec{c}^{(h)}={\bf i}_t - \varvec{s}^{(h)}$$\end{document } are positive - definite , while \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\varvec{z}^{(w)}=[z_{mf}^{(w)}]$$\end{document } and \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\varvec{z}^{(h)}=[z_{tn}^{(h)}]$$\end{document } are binary symmetric band matrices indicating the neighboring entries in \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\varvec{w}_r$$\end{document } and \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\underline { { \bf h}}_r$$\end{document } , respectively . in the first - order interactions , we have z1,2(w )  =  zf , f-1(w )  =  zm , m-1(w )  =  zm , m+1(w )  =  1 for \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$m \in \{2 , \ldots , f-1 \ } ; z_{2,1}^{(h)}=z_{t-1,t}^{(h)}=z_{n-1,n}^{(h)}=z_{n+1,n}^{(h ) } = 1$$\end{document } for \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$n \in \{2 , \ldots , t-1 \}$$\end{document } , and zmf  =  ztn  =  0 otherwise . in the p - order interactions , each entry wfr and hrt has the corresponding sets of neighbors : { wf-,r } ,  { wf+,r } ,  { hr , t- } ,  { hr , t+ } with \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\nu=1 , \ldots , p$$\end{document}. as a consequence , \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\varvec{z}^{(w)}$$\end{document } and \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\varvec{z}^{(h)}$$\end{document } are symmetric band matrices with p sub - diagonals and p super - diagonals , the entries of which are equal to ones , but zeros otherwise . the matrices \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\varvec{c}^{(w)}$$\end{document } and \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\varvec{c}^{(h)}$$\end{document } are positive - definite , if   <  ( 2p ) for p - order interactions [ 45 , 46 ] . we selected \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\gamma=(2p)^{-1 } - \tilde \epsilon$$\end{document } , where \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\tilde \epsilon$$\end{document } is a small constant , for example , \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\tilde \epsilon=10^{-16}$$\end{document}. in the sar model , gibbs priors ( 16 ) may be expressed as their joint multivariate gaussian priors:20\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ p(\varvec{w})= \prod_{r=1}^r z_w^{-1 } \exp \left \{- \frac{\alpha_w}{2 } ||\varvec{c}^{(w ) } \varvec{w}_r||_2 ^ 2 \right \ } , $ $ \end{document}21\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ p(\varvec{h})= \prod_{r=1}^r z_h^{-1 } \exp \left \{- \frac{\alpha_h}{2 } ||\underline { { \bf h}}_r \varvec{c}^{(h)}||_2 ^ 2 \right \ } $ $ \end{document}where \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$z_w=\left ( \frac{2\pi}{\alpha_w } \right)^{f/2 } \left ( \prod\nolimits_{f=1}^f \lambda_f^2(\varvec{c}^{(w ) } ) \right)^{-1/2}$$\end{document } , and \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\lambda_f(\varvec{c}^{(w)})$$\end{document } is the f - th eigenvalue of the matrix \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\varvec{c}^{(w)}$$\end{document}. similarly , \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$z_h=\left ( \frac{2\pi}{\alpha_h } \right)^{t/2 } \left ( \prod_{t=1}^t \lambda_t^2(\varvec{c}^{(h ) } ) \right)^{-1/2}$$\end{document}. if \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$p=1 : \lambda_f(\varvec{c}^{(w ) } ) \cong 1 - \cos \left ( \frac{\pi f}{f } \right)$$\end{document } and \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\lambda_t(\varvec{c}^{(h ) } ) \cong 1 - \cos \left ( \frac{\pi t}{t } \right)$$\end{document } , which simplifies the hyperparameter estimation . the em algorithm   is applied to maximize \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\ln p(\mathcal{x},\mathcal{c},\varvec{w},\varvec{h}|\mathcal{a},\varvec{\upsigma}_n)$$\end{document } in ( 11 ) . to calculate the e - step , the log - likelihood functional ( 13 ) is transformed to the following form22\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \begin{aligned } \ln p({\mathcal{x}}|{\mathcal{c}},{\mathcal{a}},\varvec{\upsigma}_n)&= - t \sum_f \hbox{tr}\left \{\varvec{\upsigma}_n^{-1 } { \bf r}_f^{(xx ) } \right \}+ t \sum_f \hbox{tr}\left \{\varvec{a}_f^h \varvec{\upsigma}_n^{-1 } { \bf r}_f^{(xs ) } \right \ } \\ & \quad + t \sum_f \hbox{tr}\left \{\varvec{\upsigma}_n^{-1 } \varvec{a}_f ( { \bf r}_f^{(xs)})^h \right \ } - t \sum_f \hbox{tr}\left \{\varvec{a}_f^h \varvec{\upsigma}_n^{-1 } \varvec{a}_f { \bf r}_f^{(ss ) } \right \ } \\ & \quad -\sum_{f , t } \ln \det \varvec{\upsigma}_n , \end{aligned } $ $ \end{document}where the correlation matrices are given by \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$${\bf r}_f^{(xx)}=\frac{1}{t } \sum_t \varvec{x}_{ft } \varvec{x}_{ft}^h , { \bf r}_f^{(ss)}=\frac{1}{t } \sum_t \varvec{s}_{ft } \varvec{s}_{ft}^h$$\end{document } , and the cross - correlation \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$${\bf r}_f^{(xs)}=\frac{1}{t } \sum_t \varvec{x}_{ft } \varvec{s}_{ft}^h$$\end{document}. ozerov et  al .   observed that the set { rf(xx ) ,  rf(xs ) ,  rf(ss ) ,  |crft| } provides sufficient statistics for the exponential family , so the sources \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\varvec{s}_{ft}$$\end{document } and the latent components \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\varvec{c}_{ft}$$\end{document } can be estimated by computing the conditional expectations of the natural statistics . according to , we have the following posterior estimates:23\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \hat{\varvec{s}}_{ft}=\varvec{\upsigma}_{ft}^{(s ) } \varvec{a}_f^h ( \varvec{a}_f \varvec{\upsigma}_{ft}^{(s ) } \varvec{a}_f^h+\varvec{\upsigma}_n)^{-1 } \varvec{x}_{ft } , $ $ \end{document}24\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \hat{\varvec{\upsigma}}^{(s)}_{ft}=\varvec{\upsigma}^{(s)}_{ft } - \varvec{\upsigma}^{(s)}_{ft } \varvec{a}_f^h ( \varvec{a}_f \varvec{\upsigma}^{(s)}_{ft } \varvec{a}_f^h+\varvec{\upsigma}_n)^{-1 } \varvec{a}_f \varvec{\upsigma}^{(s)}_{ft}. $ $ \end{document}similarly , for the latent components , we have25\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \hat{\varvec{c}}_{ft}=\varvec{\upsigma}^{(c)}_{ft } \bar{\varvec{a}}_f^h ( \bar{\varvec{a}}_f \varvec{\upsigma}^{(c)}_{ft } \bar{\varvec{a}}_f^h+\varvec{\upsigma}_n)^{-1 } \varvec{x}_{ft } , $ $ \end{document}26\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \hat{\varvec{\upsigma}}^{(c)}_{ft}=\varvec{\upsigma}^{(c)}_{ft } - \varvec{\upsigma}^{(c)}_{ft } \bar{\varvec{a}}_f^h ( \bar{\varvec{a}}_f \varvec{\upsigma}^{(c)}_{ft } \bar{\varvec{a}}_f^h+\varvec{\upsigma}_n)^{-1 } \bar{\varvec{a}}_f \varvec{\upsigma}^{(c)}_{ft}. $ $ \end{document}the conditional expectations for the sufficient statistics are as follows:27\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \hat{{\bf r}}_f^{(xx)}={\bf r}_f^{(xx ) } , \quad \hat{{\bf r}}_f^{(xs)}=\frac{1}{t } \sum_t \varvec{x}_{ft } { \mathcal{e}}({\varvec{s}}_{ft}^h)=\frac{1}{t } \sum_t \varvec{x}_{ft } \hat{\varvec{s}}_{ft}^h , $ $ \end{document}28\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \hat{{\bf r}}_f^{(ss)}=\frac{1}{t } \sum_t { \mathcal{e}}(\varvec{s}_{ft } ) { \mathcal{e}}({\varvec{s}}_{ft}^h)+\hat{\varvec{\upsigma}}^{(s)}_{ft}=\frac{1}{t } \sum_t \hat{\varvec{s}}_{ft } \hat{\varvec{s}}_{ft}^h+\hat{\varvec{\upsigma}}^{(s)}_{ft } , $ $ \end{document}29\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ |c_{rft}|^2 \leftarrow { \mathcal{e}}(c_{rft } ) { \mathcal{e}}(c_{rft}^h)+(\hat{\varvec{\upsigma}}^{(c)}_{ft})_{rr}= |\hat{c}_{rft}|^2+(\hat{\varvec{\upsigma}}^{(c)}_{ft})_{rr}. $ $ \end{document}detailed derivations of the formulae ( 23)(26 ) are presented in the `` appendix '' . from the m - step , we have \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\frac{\partial}{\partial \varvec{a}_f } \ln p(\mathcal{x},\mathcal{c},\varvec{w},\varvec{h}|\mathcal{a},\varvec{\upsigma}_n)=2t(- \varvec{\upsigma}_n^{-1 } { \bf r}^{(xs)}_f+\varvec{\upsigma}_n^{-1 } \varvec{a}_f { \bf r}^{(ss)}_f)=0$$\end{document } , which gives \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\varvec{a}_f=\hat{{\bf r}}^{(xs)}_f ( \hat{{\bf r}}^{(ss)}_f)^{-1}$$\end{document}. from\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \frac{\partial}{\partial \varvec{\upsigma}_n^{-1 } } \ln p({\mathcal{x}},{\mathcal{c}},\varvec{w},\varvec{h}|{\mathcal{a}},\varvec{\upsigma}_n)=0 , $ $ \end{document}we have\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \varvec{\upsigma}_n=\hbox{diag}\left \{{\bf r}^{(xx)}_f - \varvec{a}_f ( \hat{{\bf r}}^{(xs)}_f)^h - \hat{{\bf r}}^{(xs)}_f \varvec{a}_f^h+\varvec{a}_f \hat{{\bf r}}^{(ss)}_f \varvec{a}_f^h \right \}. $ $ \end{document}from \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\frac{\partial}{\partial w_{fr } } \ln p(\mathcal{x},\mathcal{c},\varvec{w},\varvec{h}|\mathcal{a},\varvec{\upsigma}_n)=0$$\end{document } , we have30\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ w_{fr}=\frac{1}{t } \sum_{t=1}^t \frac{|c_{rft}|^2}{h_{rt } } - \alpha_w \nabla_{w_{fr } } u(\varvec{w } ) . $ $ \end{document}similarly , from \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\frac{\partial}{\partial h_{rt } } \ln p(\mathcal{x},\mathcal{c},\varvec{w},\varvec{h}|\mathcal{a},\varvec{\upsigma}_n)=0$$\end{document } , we get31\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ h_{rt}=\frac{1}{f } \sum_{f=1}^f \frac{|c_{rft}|^2}{w_{fr } } - \alpha_h \nabla_{h_{rt } } u(\varvec{w } ) . $ $ \end{document}the terms \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\nabla_{w_{fr } } u(\varvec{w})$$\end{document } and \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\nabla_{h_{rt } } u(\varvec{w})$$\end{document } in ( 30 ) and ( 31 ) take the following forms with respect to the potential functions : gaussian ( sar model):32\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \nabla_{w_{fr } } u(\varvec{w})= \left [ ( \varvec{c}^{(w)})^t \varvec{c}^{(w ) } \varvec{w } \right]_{fr } , $ $ \end{document}33\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \nabla_{h_{rt } } u(\varvec{h})= \left [ \varvec{h } \varvec{c}^{(h ) } ( \varvec{c}^{(h)})^t \right]_{rt } , $ $ \end{document}gr function ( proposed by green ) : 34\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \nabla_{w_{fr } } u(\varvec{w})= \sum_{l \in s_f } \nu_{fl } \tanh \left ( \frac{w_{fr } - w_{lr}}{\delta_w } \right ) , $ $ \end{document}35\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \nabla_{h_{rt } } u(\varvec{h})= \sum_{l \in s_t } \nu_{tl } \tanh \left ( \frac{h_{rt } - h_{rl}}{\delta_h } \right ) . $ $ \end{document}hl function ( proposed by hebert and leahy ) : 36\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \nabla_{w_{fr } } u(\varvec{w})= \sum_{l \in s_f } \nu_{fl } \frac{2 \delta_w ( w_{fr } - w_{lr})}{\delta^2_w+(w_{fr } - w_{lr})^2 } , $ $ \end{document}37\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \nabla_{h_{rt } } u(\varvec{h})= \sum_{l \in s_t } \nu_{tl } \frac{2 \delta_h ( h_{rt } - h_{rl})}{\delta^2_h+(h_{rt } - h_{rl})^2}. $ $ \end{document } gaussian ( sar model):32\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \nabla_{w_{fr } } u(\varvec{w})= \left [ ( \varvec{c}^{(w)})^t \varvec{c}^{(w ) } \varvec{w } \right]_{fr } , $ $ \end{document}33\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \nabla_{h_{rt } } u(\varvec{h})= \left [ \varvec{h } \varvec{c}^{(h ) } ( \varvec{c}^{(h)})^t \right]_{rt } , $ $ \end{document } gr function ( proposed by green ) : 34\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \nabla_{w_{fr } } u(\varvec{w})= \sum_{l \in s_f } \nu_{fl } \tanh \left ( \frac{w_{fr } - w_{lr}}{\delta_w } \right ) , $ $ \end{document}35\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \nabla_{h_{rt } } u(\varvec{h})= \sum_{l \in s_t } \nu_{tl } \tanh \left ( \frac{h_{rt } - h_{rl}}{\delta_h } \right ) . $ $ \end{document } hl function ( proposed by hebert and leahy ) : 36\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \nabla_{w_{fr } } u(\varvec{w})= \sum_{l \in s_f } \nu_{fl } \frac{2 \delta_w ( w_{fr } - w_{lr})}{\delta^2_w+(w_{fr } - w_{lr})^2 } , $ $ \end{document}37\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \nabla_{h_{rt } } u(\varvec{h})= \sum_{l \in s_t } \nu_{tl } \frac{2 \delta_h ( h_{rt } - h_{rl})}{\delta^2_h+(h_{rt } - h_{rl})^2}. $ $ \end{document } experiments were conducted using selected sound recordings taken from the stereo audio source separation evaluation campaign ( sisec)1 in 2007 . this campaign aimed to evaluate the performance of source separation algorithms using stereo under - determined mixtures . we selected the benchmarks given in table  2 , which included speech recordings ( three male voices female3 ) , three nonpercussive music sources  nodrums , and three music sources that included drums  wdrums . the mixed signals were recordings that lasted 10  s , which were sampled at 16  khz ( the standard settings of recordings from the  under - determined speech and music mixtures " datasets in the sisec2008 ) . for each benchmark , the number of true sources was three ( j  =  3 ) but it only had two microphones ( i  =  2 ) , that is , stereo recordings . synthetic convolutive mixtures were obtained for a meeting room with a 250  ms reverberation time using omnidirectional microphones with 1  m spacing.table  2benchmarksinstantaneousconvolutive 
 male3_inst_mix 

 male3_synthconv_250ms_1m_mix 

 female3_inst_mix 

 female3_synthconv_250ms_1m_mix 

 nodrums_inst_mix 

 nodrums_synthconv_250ms_1m_mix 

 wdrums_inst_mix wdrums_synthconv_250ms_1m_mix 
  the spectrograms were obtained by a short - time fourier transform ( stft ) using half - overlapping sine windows . to create the spectrograms and recover the time - domain signals from stft coefficients , we used the corresponding stft_multi and istft_multi matlab functions from the sisec2008 webpage2 . for instantaneous and convolutive mixtures , the em - nmf algorithm was taken from ozerov s homepage3 , while the mrf - em - nmf algorithm was coded and extensively tested by ochal . the proposed algorithm is based on an alternating optimization scheme , which is intrinsically non - convex , and hence , its initialization plays an important role . an incorrect initialization may result in slow convergence and early stagnation at an unfavorable local minimum of the objective function . as done in many nmf algorithms , the factors w and h are initialized with uniformly distributed random numbers , whereas the entries in the matrix a are drawn from a zero - mean complex gaussian distribution . after w and h have been initialized , the covariance matrices \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\varvec{\upsigma}_{ft}^{(s)}$$\end{document } and \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\varvec{\upsigma}_{ft}^{(c)}$$\end{document } given by ( 8) can be computed . a noise covariance matrix \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\varvec{\upsigma}_n$$\end{document } is needed to update the e - step . the e - step in mrf - em - nmf is identical to that in em - nmf , and hence , all of these techniques can be used in this experiment . the initial matrix \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\varvec{\upsigma}_n$$\end{document } was determined based on the empirical variance of the observed power spectrograms . the mrf - em - nmf and em - nmf algorithms were initialized using the same random values ( given as \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\bar{r}$$\end{document } ) and run for 1,500 iterations . the choice of the parameters { w ,  h ,  w ,  h } used in the gibbs distributions also affected the performance . motivated by iterative thresholding strategies , we used the following strategies : linear thresholding:\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \alpha(k)=\alpha \frac{k}{k_{\rm max } } , $ $ \end{document}nonlinear thresholding:\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \alpha(k)=\frac{\alpha}{2 } \left ( 1+\tanh \left ( \frac{k - \nu k_{\rm max}}{\tau k_{\rm max } } \right ) \right ) , $ $ \end{document}fixed thresholding:\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \alpha(k)=\left \{\begin{array}{cc}\alpha & \hbox{if } \quad k > k_1 , \\ 0 & \hbox{otherwise } \end{array } \right . $ $ \end{document}where k is the current iteration , kmax is the maximum number of iterations , \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\tau \in ( 0,1)$$\end{document } is the shape parameter , \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\nu \in ( 0,1)$$\end{document } is the shift parameter , k1 is the threshold , and  can be equal to w or h . all of the above thresholding strategies aim to relax smoothing during the early iterations when the descent directions in the updates are sufficiently steep and to emphasize smoothing if noisy perturbations become significantly detrimental to the overall smoothness . these strategies are motivated by standard regularization rules that apply to ill - posed problems . we tested all of the thresholding strategies using instantaneous and convolutive mixtures , and we obtained the best performance with fixed thresholding using k1  =  kmax/2 . linear thresholding:\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \alpha(k)=\alpha \frac{k}{k_{\rm max } } , $ $ \end{document } nonlinear thresholding:\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \alpha(k)=\frac{\alpha}{2 } \left ( 1+\tanh \left ( \frac{k - \nu k_{\rm max}}{\tau k_{\rm max } } \right ) \right ) , $ $ \end{document } fixed thresholding:\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \alpha(k)=\left \{\begin{array}{cc}\alpha & \hbox{if } \quad k > k_1 , \\ 0 & \hbox{otherwise } \end{array } \right . $ $ \end{document}where k is the current iteration , kmax is the maximum number of iterations , \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\tau \in ( 0,1)$$\end{document } is the shape parameter , \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\nu \in ( 0,1)$$\end{document } is the shift parameter , k1 is the threshold , and  can be equal to w or h . all of the above thresholding strategies aim to relax smoothing during the early iterations when the descent directions in the updates are sufficiently steep and to emphasize smoothing if noisy perturbations become significantly detrimental to the overall smoothness . these strategies are motivated by standard regularization rules that apply to ill - posed problems . we tested all of the thresholding strategies using instantaneous and convolutive mixtures , and we obtained the best performance with fixed thresholding using k1  =  kmax/2 . the parameters w and h in the mrf models can be estimated using standard marginalization procedures or by maximizing the type ii ml estimate for ( 10 ) . however , these techniques have a huge computational cost for the nonlinear potential functions in the mrf models . for practical reasons , they are not very useful for the gr or hr functions . in this study , we tested all of the benchmarks in table  2 and the following potential functions : the first- and second - order gaussian , gr , and hr . for the gaussian functions , we tested all combinations of the regularization parameters w and h from the discrete set { 0.001 , 0.005 , 0.01 , 0.05 , 0.1}. for gr and hl , the regularization parameters could take only two values , { 0.001 , 0.01 } , although the parameters w and h were tested with the following values : { 0.1 , 1 , 10}. the optimal values of the smoothing parameters are summarized in table  3.table  3parameters of the mrf - em - nmf algorithm for each test case shown in fig .  1 
 benchmarksmoothinginstantaneous mixtureconvolutive mixture 
 \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\bar{r}$$\end{document } 
 w 
 h w 
 h 

 \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\bar{r}$$\end{document } 
 w 
 h 
 w 
 h 
 malegr120.010.011140.010.010.110malehl120.0010.00111040.0010.0111male1-gaussian120.0010.0140.050.05male2-gaussian120.0010.0140.050.01femalegr120.010.01101040.010.0111femalehl120.0010.00111040.0010.0010.110female1-gaussian120.0010.00140.10.001female2-gaussian120.0010.00140.050.005nodrumsgr40.010.0110140.010.01100.1nodrumshl40.010.00111040.010.010.10.1nodrums1-gaussian40.0010.0140.0010.05nodrums2-gaussian40.010.00140.0050.01wdrumsgr40.010.0111040.010.0111wdrumshl40.010.00111040.0010.0110.1wdrums1-gaussian40.0010.00140.0010.1wdrums2-gaussian40.0010.00140.0050.1the notations " 1-gaussian " and " 2-gaussian " represent the first- and second - order gaussian functions , respectively parameters of the mrf - em - nmf algorithm for each test case shown in fig .  1 
  the notations " 1-gaussian " and " 2-gaussian " represent the first- and second - order gaussian functions , respectively the separation results were evaluated in terms of the signal - to - distortion ratio ( sdr ) and the signal - to - interference ratio ( sir ) . figure  1 shows the sdrs and sirs averaged for the sources , which were estimated using the em - nmf and mrf - em - nmf with various smoothing functions based on instantaneous and convolutive mixing models . for each sample in table  2 and each smoothing function , the smoothing parameters were tuned optimally for a given fixed initializer . this unsupervised learning approach evaluated the efficiency of the smoothing functions with respect to a given recording scenario . however , the smoothing parameters need to be determined with a supervised learning framework in practice . to test this option , each recording in table  2 was divided into two 5  s excerpts during the training and testing stages . for each training excerpt , the smoothing parameters and initializer were selected to maximize the sdr performance . the results obtained during the testing stage with the instantaneous mixtures are shown in fig . 1source separation results obtained with the mrf - em - nmf ( first- and second - order gaussian , gr , and hl functions ) and em - nmf ( no smoothing ) algorithms after 1,500 iterations : a mean sdr ( db ) for instantaneous mixture , b mean sdr ( db ) for convolutive mixture , c mean sir ( db ) for instantaneous mixture , d mean sir ( db ) for convolutive mixture . the smoothing parameters were tuned separately for each mixture in table  2 
 fig . 2source separation results obtained in the testing stage with the mrf - em - nmf ( first- and second - order gaussian , gr , and hl functions ) and em - nmf ( no smoothing ) algorithm after 1,500 iterations : a mean sdr ( db ) , b mean sir ( db ) . 5  s excerpts were used in the training and testing stages source separation results obtained with the mrf - em - nmf ( first- and second - order gaussian , gr , and hl functions ) and em - nmf ( no smoothing ) algorithms after 1,500 iterations : a mean sdr ( db ) for instantaneous mixture , b mean sdr ( db ) for convolutive mixture , c mean sir ( db ) for instantaneous mixture , d mean sir ( db ) for convolutive mixture . the smoothing parameters were tuned separately for each mixture in table  2 
  source separation results obtained in the testing stage with the mrf - em - nmf ( first- and second - order gaussian , gr , and hl functions ) and em - nmf ( no smoothing ) algorithm after 1,500 iterations : a mean sdr ( db ) , b mean sir ( db ) . 5  s excerpts were used in the training and testing stages for comparison , table  4 shows the average sdr results produced and the running time taken when using several state - of - the - art algorithms , which were applied to the mixtures in table  2 . the generalized gaussian prior ( ggp ) algorithm   and the statistically sparse decomposition principle ( ssdp ) algorithms   were applied to the instantaneous mixtures . the convolutive mixtures were unmixed with the ipd , two versions of the fbwc - pa [ 17 , 18 ] algorithm , and the convolutive nmf . note that the last method in this list was based on supervised learning , whereas the others were unsupervised learning algorithms . in this case , the first 8  s excerpts of the 10  s source recordings were used for learning , while the remainder was used for testing.table  4mean sdr ( db ) and running time ( s ) for sources estimated from the mixtures shown in table  2 
 benchmarkmixturemalefemalenodrumswdrumstimemrf - em - nmf ( hr)inst8.069.9524.0721.722487mrf - em - nmf ( gr ) inst7.698.8626.6521.282498em - nmf inst2.626.511.719.872456ggp inst8.48.5713.910.35sabm+ssdp inst4.253.825.839.432mrf - em - nmf ( hr)conv1.062.21.171.72760mrf - em - nmf ( gr ) conv1.42.11.21.562762em - nmf conv0.951.60.20.442720ipd conv1.531.432.22.71200fbwc - pa conv0.14.430.772.5340generalized fbwc - pa conv5.957.451.20.698convnmf conv0.70.473.858.13347 mean sdr ( db ) and running time ( s ) for sources estimated from the mixtures shown in table  2 
  the averaged elapsed time measured using matlab 2008a for 1,500 iterations with \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\bar{r}=12$$\end{document } , executed on a 64-bit intel quad core cpu 3  ghz with 8  gb ram was almost the same for the mrf - em - nmf and em - nmf algorithms ( see table  4 ) . the simulations demonstrate that mrf smoothing improved the source separation results in almost all test cases . the mrf - em - nmf algorithm delivered the best mean sdr performance of all the algorithms tested with instantaneous mixtures . this was justified by the smooth frequency and temporal structures of non - percussive music spectrograms . if the source spectrograms were not very smooth ( as with the percussive audio recordings ) , mrf smoothing gave only a slight improvement ( see figs .  1 , 2 ) in the first - order mrf interactions , and even a slight deterioration in the higher - order mrf interactions . according to fig .  1 , the hl function delivered the most promising sdr results , which were stable with a wide range of parameters . in each case with the instantaneous mixtures , the best results were produced with the same hyperparameter values , w  =  1 and h  =  10 , and almost the same penalty parameter values , w and h . the sar model also improved the results compared with the standard em - nmf algorithm . moreover , the sar model was tuned using only two penalty parameters , and the partition function of the associated gibbs prior could be derived using a closed - form expression , which might be very useful for data - driven hyperparameter estimation . the source separation results produced with the mrf - em - nmf algorithm for convolutive and under - determined mixtures were better than those obtained with the em - nmf algorithm . unfortunately , the sdr values showed that these results were still a long way from being perfect , even after 1,500 iterations , and thus , further research is needed in this field . it is likely that some additional prior information could be imposed , especially on a mixing operator , which might increase the efficiency considerably . it should be noted that the sdr performance with both mixtures could still be improved by refining the associated parameters , especially in the mrf models , and by using more efficient initializers . this study demonstrated that imposing mrf smoothing on the power spectrograms of audio sources estimated from under - determined unmixing problems may improve the quality of estimated audio sounds considerably . this was justified because any type of meaningful prior information improves the performance , especially with under - determined problems . this study addressed the application of mrf smoothing in the em - nmf algorithm , but this type of smoothing could be applied to many other related bss algorithms based on feature extraction from power spectrograms . clearly , further studies are needed to improve this technique for convolutive mixtures and to integrate regularization parameter estimation techniques in the main algorithm . the conditional expectations of the natural statistics can be derived from the a posteriori distributions \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$p(\varvec{s}_{ft}|\varvec{x}_{ft})$$\end{document } and \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$p(\varvec{c}_{ft}|\varvec{x}_{ft})$$\end{document}. thus,38\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \begin{aligned } p(\varvec{s}_{ft}|\varvec{x}_{ft})&= \frac{p(\varvec{x}_{ft } , \varvec{s}_{ft})}{p(\varvec{x}_{ft } ) } \\ & = \frac{(\pi^{i+j } \det \varvec{\upsigma}_{ft})^{-1 } \exp \left \{- \left [ \begin{array}{c } \bar{\varvec{x}}_{ft } \\ \bar{\varvec{s}}_{ft } \end{array } \right]^h ( \varvec{\upsigma}_{ft})^{-1 } \left [ \begin{array}{c } \bar{\varvec{x}}_{ft } \\ \bar{\varvec{s}}_{ft } \end{array } \right ] \right \}}{(\pi^{i } \det \varvec{\upsigma}_{ft}^{(x)})^{-1 } \exp \left \{- ( \bar{\varvec{x}}_{ft})^h ( \varvec{\upsigma}_{ft}^{(x)})^{-1 } \bar{\varvec{x}}_{ft } \right \ } } \\ & = ( \pi^{j } \det { \bf \upgamma}_{ft})^{-1 } \exp \left \{- { \bf \uppsi}_{ft } \right \ } , \end{aligned } $ $ \end{document}where \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\bar{\varvec{x}}_{ft}=\varvec{x}_{ft } - \mathcal{e}(\varvec{x}_{ft } ) , \bar{\varvec{s}}_{ft}=\varvec{s}_{ft } - \mathcal{e}(\varvec{s}_{ft } ) , { \bf \upgamma}_{ft}=\varvec{\upsigma}_{ft}^{(s ) } - \varvec{\upsigma}_{ft}^{(sx ) } ( \varvec{\upsigma}_{ft}^{(x)})^{-1 } \varvec{\upsigma}_{ft}^{(xs)}$$\end{document},\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \varvec{\upsigma}_{ft}=\left [ \begin{array}{cc } \varvec{\upsigma}_{ft}^{(x ) } & \varvec{\upsigma}_{ft}^{(xs ) } \\ \varvec{\upsigma}_{ft}^{(sx ) } & \varvec{\upsigma}_{ft}^{(s ) } \end{array } \right ] , \quad { \bf \uppsi}_{ft}=\left [ \begin{array}{c } \bar{\varvec{x}}_{ft } \\ \bar{\varvec{s}}_{ft } \end{array } \right]^h ( \varvec{\upsigma}_{ft})^{-1 } \left [ \begin{array}{c } \bar{\varvec{x}}_{ft } \\ \bar{\varvec{s}}_{ft } \end{array } \right ] - ( \bar{\varvec{x}}_{ft})^h ( \varvec{\upsigma}_{ft}^{(x)})^{-1 } \bar{\varvec{x}}_{ft}. $ $ \end{document}we can transform \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\varvec{\upsigma}_{ft}^{-1}$$\end{document } in ( 38 ) into the following form:\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \varvec{\upsigma}_{ft}^{-1}=\left [ \begin{array}{cc } \left ( \varvec{\upsigma}_{ft}^{(x ) } - \varvec{\upsigma}_{ft}^{(xs ) } ( \varvec{\upsigma}_{ft}^{(s)})^{-1 } \varvec{\upsigma}_{ft}^{(sx ) } \right)^{-1 } & - ( \varvec{\upsigma}_{ft}^{(x)})^{-1 } \varvec{\upsigma}_{ft}^{(xs ) } { \bf \upgamma}_{ft}^{-1 } \\ - { \bf \upgamma}_{ft}^{-1 } \varvec{\upsigma}_{ft}^{(sx)}(\varvec{\upsigma}_{ft}^{(x)})^{-1 } & { \bf \upgamma}_{ft}^{-1 } \end{array } \right ] . we have\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ ( \varvec{\upsigma}_{ft}^{(x ) } - \varvec{\upsigma}_{ft}^{(xs ) } ( \varvec{\upsigma}_{ft}^{(s)})^{-1 } \varvec{\upsigma}_{ft}^{(sx)})^{-1}=(\varvec{\upsigma}_{ft}^{(x)})^{-1}+(\varvec{\upsigma}_{ft}^{(x)})^{-1 } \varvec{\upsigma}_{ft}^{(xs ) } { \bf \upgamma}_{ft}^{-1 } \varvec{\upsigma}_{ft}^{(sx ) } ( \varvec{\upsigma}_{ft}^{(x)})^{-1 } , $ $ \end{document}and finally,39\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \begin{aligned } { \bf \uppsi}_{ft}&= \left ( \bar{\varvec{s}}_{ft } - \varvec{\upsigma}_{ft}^{(sx ) } ( \varvec{\upsigma}_{ft}^{(x)})^{-1 } \bar{\varvec{x}}_{ft } \right)^h { \bf \upgamma}_{ft}^{-1 } \left ( \bar{\varvec{s}}_{ft } - \varvec{\upsigma}_{ft}^{(sx ) } ( \varvec{\upsigma}_{ft}^{(x)})^{-1 } \bar{\varvec{x}}_{ft } \right ) \\ & = \left ( { \varvec{s}}_{ft } - \hat{\varvec{s}}_{ft } \right)^h ( \hat{\varvec{\upsigma}}^{(s)}_{ft})^{-1 } \left ( { \varvec{s}}_{ft } - \hat{\varvec{s}}_{ft } \right ) , \end{aligned } $ $ \end{document}where40\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \hat{\varvec{s}}_{ft}={\mathcal{e}}(\varvec{s}_{ft})+\varvec{\upsigma}_{ft}^{(sx ) } ( \varvec{\upsigma}_{ft}^{(x)})^{-1 } \left ( \varvec{x}_{ft } - { \mathcal{e}}(\varvec{x}_{ft})\right ) , $ $ \end{document}41\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \hat{\varvec{\upsigma}}^{(s)}_{ft}={\bf \upgamma}_{ft}=\varvec{\upsigma}_{ft}^{(s ) } - \varvec{\upsigma}_{ft}^{(sx ) } ( \varvec{\upsigma}_{ft}^{(x)})^{-1 } \varvec{\upsigma}_{ft}^{(xs)}. $ $ \end{document}thus , \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$p(\varvec{s}_{ft}|\varvec{x}_{ft})=\mathcal{n}_c(\varvec{s}_{ft } ; \hat{\varvec{s}}_{ft } , \hat{\varvec{\upsigma}}^{(s)}_{ft})$$\end{document}. from ( 5 ) , it follows that \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\mathcal{e}(\varvec{s}_{ft})=0$$\end{document } , so \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\mathcal{e}(\varvec{x}_{ft})=0$$\end{document}. since the zero - mean noise nft from ( 3 ) is not correlated with \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\varvec{s}_{ft}$$\end{document } , we have42\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \begin{aligned } \varvec{\upsigma}_{ft}^{(xs)}&= { \mathcal{e } } \left ( ( \varvec{x}_{ft } - { \mathcal{e}}(\varvec{x}_{ft } ) ) ( \varvec{s}_{ft } - { \mathcal{e}}(\varvec{s}_{ft}))^h \right ) \\ & = { \mathcal{e } } \left ( ( \varvec{a}_f \varvec{s}_{ft}+{\bf n}_{ft})\varvec{s}_{ft}^h \right)=\varvec{a}_f { \mathcal{e}}(\varvec{s}_{ft}\varvec{s}_{ft}^h)+{\mathcal{e}}({\bf n}_{ft}\varvec{s}_{ft}^h)=\varvec{a}_f \varvec{\upsigma}_{ft}^{(s)}. \end{aligned } $ $ \end{document}inserting ( 42 ) and \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\varvec{\upsigma}_{ft}^{(sx)}=(\varvec{\upsigma}_{ft}^{(xs)})^h= \varvec{\upsigma}_{ft}^{(s ) } \varvec{a}_f^h$$\end{document } into ( 40 ) and ( 41 ) , we obtain the update rules ( 23 ) and ( 24 ) , respectively . analyzing \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$p(\varvec{c}_{ft}|\varvec{x}_{ft})$$\end{document } , one can obtain \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$p(\varvec{c}_{ft}|\varvec{x}_{ft})=\mathcal{n}_c(\varvec{c}_{ft } ; \hat{\varvec{c}}_{ft } , \hat{\varvec{\upsigma}}^{(c)}_{ft})$$\end{document } , which yields the update rules ( 25 ) and ( 26 ) .
convolutive and under - determined blind audio source separation from noisy recordings is a challenging problem . 
 several computational strategies have been proposed to address this problem . 
 this study is concerned with several modifications to the expectation - minimization - based algorithm , which iteratively estimates the mixing and source parameters . 
 this strategy assumes that any entry in each source spectrogram is modeled using superimposed gaussian components , which are mutually and individually independent across frequency and time bins . in our approach , we resolve this issue by considering a locally smooth temporal and frequency structure in the power source spectrograms . 
 local smoothness is enforced by incorporating a gibbs prior in the complete data likelihood function , which models the interactions between neighboring spectrogram bins using a markov random field . 
 simulations using audio files derived from stereo audio source separation evaluation campaign 2008 demonstrate high efficiency with the proposed improvement .