Multichannel Audio Source Separation with Probabilistic Reverberation Priors

Simon Leglaive, Roland Badeau, Senior Member, IEEE, Gaël Richard Senior Member, IEEE Abstract-Incorporating prior knowledge about the sources and/or the mixture is a way to improve under-determined audio source separation performance. A great number of informed source separation techniques concentrate on taking priors on the sources into account, but fewer works have focused on constraining the mixing model. In this paper we address the problem of under-determined multichannel audio source separation in reverberant conditions. We target a semi-informed scenario where some room parameters are known. Two probabilistic priors on the frequency response of the mixing filters are proposed. Early reverberation is characterized by an autoregressive model while according to statistical room acoustics results, late reverberation is represented by an autoregressive moving average model. Both reverberation models are defined in the frequency domain. They aim to transcribe the temporal characteristics of the mixing filters into frequency-domain correlations. Our approach leads to a maximum a posteriori estimation of the mixing filters which is achieved thanks to an expectation-maximization algorithm. We experimentally show the superiority of this approach compared with a maximum likelihood estimation of the mixing filters.

Index Terms-Multichannel audio source separation, probabilistic priors, mixing model, MAP estimation, EM algorithm.

I. INTRODUCTION

A UDIO source separation is the task that aims to recover a set of source signals from the observation of one or several mixtures. Audio source separation can be used as a pre-processing step for classification or recognition tasks. Indeed, more relevant features can be computed on the separated sources or higher-level information can be extracted from parametric models used for source separation. In music information retrieval, audio source separation can for example help improving the performance of tasks such as instrument recognition [START_REF] Heittola | Musical instrument recognition in polyphonic audio using source-filter model for sound separation[END_REF], main melody extraction [START_REF] Durrieu | A musically motivated midlevel representation for pitch estimation and musical audio source separation[END_REF], [START_REF] Tachibana | Melody line estimation in homophonic music audio signals based on temporal-variability of melodic source[END_REF] or singing voice detection [START_REF] Leglaive | Singing voice detection with deep recurrent neural networks[END_REF]. Audio source separation is more challenging when the separated sources are meant to be listened to by humans, such as for interference removal in close-microphone live recordings, remixing or upmixing applications.

When a punctual source is emitting in an enclosed space, the signal recorded by a microphone will include multiple reflections of the source signal on the surfaces and objects in the room. The acoustic channel between a source and a microphone is characterized by a Room Impulse Response (RIR) and represents how the sound propagates in the room. In S. Leglaive, R. Badeau and G. Richard are with LTCI, CNRS, Télécom ParisTech, Université Paris-Saclay, 75013, Paris, France (e-mail: firstname.lastname@telecom-paristech.fr).

The research leading to this paper was partly supported by the French National Research Agency (ANR) as a part of the EDISON 3D project (ANR-13-CORD-0008-02).

the context of audio source separation the mixture is referred to as convolutive due to the filters involved in the mixing process. Each RIR between a microphone and a source is called a mixing filter. Non-punctual source models have also been proposed using spatial covariance matrices in [START_REF] Vincent | Probabilistic modeling paradigms for audio source separation[END_REF], [START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF].

Source separation is commonly achieved in a Time-Frequency (TF) domain. Firstly because using a short-term representation can be useful to describe the spectro-temporal evolution of the sources. Secondly because it is possible to exploit the sparsity of the sources in this domain. And thirdly because expressing the mixture in this domain is easier: it can be defined according to a frequency-dependent mixing matrix built from the frequency response of the mixing filters.

Recent and now widely used approaches for tackling underdetermined audio source separation are based on variance modeling frameworks [START_REF] Vincent | Probabilistic modeling paradigms for audio source separation[END_REF]. In these probabilistic approaches, the Short-Term Fourier Transform (STFT) coefficients of each source are modeled as latent random variables following a complex circularly symmetric distribution with a time and frequency-dependent scale parameter. Within this framework, Non-negative Matrix Factorization (NMF) techniques are popular to represent the spectro-temporal characteristics of the sources [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence: With application to music analysis[END_REF]- [START_REF] Simsekli | Alpha-stable matrix factorization[END_REF]. The interpretability of a non-negative decomposition makes it possible to incorporate deterministic constraints or probabilistic priors on the sources [START_REF] Ozerov | A general flexible framework for the handling of prior information in audio source separation[END_REF], [START_REF] Vincent | From blind to guided audio source separation: How models and side information can improve the separation of sound[END_REF]. One can for example consider specific spectral structures as a source/filter model [START_REF] Durrieu | A musically motivated midlevel representation for pitch estimation and musical audio source separation[END_REF]. User input can be taken into account in the form of humming [START_REF] Smaragdis | Separation by humming: User-guided sound extraction from monophonic mixtures[END_REF]. Even the musical score can be incorporated within the decomposition [START_REF] Ganseman | Evaluation of a score-informed source separation system[END_REF], [START_REF] Hennequin | Score informed audio source separation using a parametric model of non-negative spectrogram[END_REF]. Deep neural networks have also been used recently in the context of these variance modeling frameworks [START_REF] Nugraha | Multichannel audio source separation with deep neural networks[END_REF].

Compared with the large number of constraints that have been proposed about the sources, only a few methods consider priors on the spatial mixing parameters (i.e. the mixing filters or the spatial covariance matrices). In [START_REF] Otsuka | Multichannel sound source dereverberation and separation for arbitrary number of sources based on Bayesian nonparametrics[END_REF] the authors consider a complex Wishart prior on the inverse of the spatial covariance matrices. Geometrically calculated or pre-measured steering vectors representing the direct path between the source and the microphones are used to parametrize the prior. In [START_REF] Duong | Spatial location priors for gaussian model based reverberant audio source separation[END_REF] the authors consider Inverse-Wishart and Gaussian priors over the spatial covariance matrices. The prior aims to represent that, on average over all possible source and microphone positions in the room, the spatial covariance matrix is equal to the sum of two terms: the outer-product of the steering vector and a term modeling the spatial correlations between the microphones due to late reverberation. The method requires the knowledge of the source positions and certain room characteristics.

In this article we propose two probabilistic priors over the frequency response of the mixing filters. The first one corresponds to an autoregressive model of early reverberation. It is based on our previous works [START_REF] Leglaive | A priori probabiliste anéchoïque pour la séparation sous-déterminée de sources sonores en milieu réverbérant[END_REF], [START_REF]Multichannel audio source separation with probabilistic reverberation modeling[END_REF] in which we neglected the influence of late reverberation. This is not the case in this paper as we consider a second prior for late reverberation. This prior aims to transcribe the time-domain exponential decay of the late reverberation power into frequency-domain correlations. Based on our late reverberation model in the frequency domain [START_REF]Autoregressive moving average modeling of late reverberation in the frequency domain[END_REF], we use an Autoregressive Moving-Average (ARMA) representation to characterize the late part of the mixing filters. This prior requires the knowledge of the reverberation time, the volume of the room and the total wall area. The priors are taken into account within an Expectation-Maximization (EM) algorithm so as to perform a Maximum A Posteriori (MAP) estimation of the mixing filters. We show that the source separation performance is improved compared to the standard Maximum Likelihood (ML) estimation.

Section II introduces the general source separation framework on which our work is built. It is mainly based on the previous methods [START_REF] Ozerov | Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation[END_REF], [START_REF] Ozerov | Multichannel nonnegative tensor factorization with structured constraints for userguided audio source separation[END_REF] where the mixing filters are estimated in the ML sense. In section III we introduce the reverberation models that will be used in section IV to define the priors on the mixing filters. The models are shown to be consistent with real data. Source separation experiments are detailed in section V where we compare the performance for ML and MAP estimations of the mixing filters. We finally draw conclusions in section VI.

II. SOURCE SEPARATION FRAMEWORK

We consider a mixture of J source signals s j (t) ∈ R on I channels denoted by x(t) = [x 1 (t), ..., x I (t)] T ∈ R I where (•) T is the transposition operator. The mixture can be decomposed according to the following additive model:

x(t) = J j=1 y j (t) + b(t), (1)
where y j (t) = [y 1j (t), ..., y Ij (t)] T ∈ R I is the j-th source image, i.e. the vector of size I containing the image of the source signal s j (t) at each microphone i, and b

(t) = [b 1 (t), ..., b I (t)] T ∈ R I is

A. Mixing Model

Using the STFT, we can write the mixture for all (f, n) ∈ {0, ..., F -1} × {0, ..., N -1} as:

x f n = J j=1 y j,f n + b f n , (2)
where

x f n = [x 1,f n , ..., x I,f n ] T , y j,f n = [y 1j,f n , ..., y Ij,f n ] T and b f n = [b 1,f n , ..., b I,f n ]
T are the complex-valued vectors containing the STFT coefficients 1 of x i (t), y ij (t) and b i (t) respectively, at TF point (f, n). Each source image can be characterized by the associated source signal and a set of mixing filters. A widely used approximation in source separation consists in considering that the mixing filters are short compared with the length of the STFT analysis window. Under this hypothesis the convolutive mixing can be approximated by an instantaneous mixing in each frequency band [START_REF] Parra | Convolutive blind separation of non-stationary sources[END_REF]. A source image thus writes:

y j,f n = a j,f s j,f n , (3)
where a j,f = [a 1j,f , ..., a Ij,f ] T ∈ C I contains the frequency response of the mixing filters between source j and the I microphones and s j,f n ∈ C represents the STFT of source j. This short mixing filters assumption is experimentally discussed in section V-C. From ( 2) and (3) the mixture can be written in matrix form:

x f n = A f s f n + b f n , (4)
where

A f = [a ij,f ] ij ∈ C
I×J is referred to as the mixing matrix and s f n = [s 1,f n , ..., s J,f n ] T ∈ C J contains the coefficients of the J source STFTs. Compared with the mixing equation ( 2) which only involves the source images, the convolutive model (4) highlights the way the source signals are mixed, through the introduction of the mixing matrix. This formalism is thus well suited for representing the physics of the mixing process such as the propagation of sound in a room. However it reaches its limits when the sources cannot be considered as punctual or when the mixing process involves long reverberations for example. In these cases, other models less oriented on the representation of the mixing physical phenomena can be used. One possible approach consists in directly modeling the source images using spatial covariance matrices [START_REF] Vincent | Probabilistic modeling paradigms for audio source separation[END_REF], [START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF]. As we focus in this paper on incorporating physically motivated priors on the mixing filters, we will only consider the convolutive model (4).

B. Local Gaussian Source Model 1) Source Distribution:

The local Gaussian model is widely used in audio source separation [START_REF] Févotte | Maximum likelihood approach for blind audio source separation using time-frequency gaussian source models[END_REF], [START_REF] Ozerov | Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation[END_REF], [START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF], [START_REF] Ozerov | A general flexible framework for the handling of prior information in audio source separation[END_REF], [START_REF] Kounades-Bastian | A variational EM algorithm for the separation of timevarying convolutive audio mixtures[END_REF]. Each coefficient s j,f n is modeled as a random variable following a centered proper complex Gaussian distribution:

s j,f n ∼ N c (0, v j,f n ). (5)
N c (µ, Σ) is the multivariate proper complex Gaussian distribution with Probability Density Function (PDF):

N c (x; µ, Σ) = 1 det(πΣ) exp[-(x -µ) H Σ -1 (x -µ)], ( 6
)
where Σ is the positive definite covariance matrix, (•) H is the Hermitian transposition operator and det(•) the matrix determinant. The term proper means that the pseudocovariance matrix E[(x -µ)(xµ) T ] is zero where E[•] is 1 Due to the Hermitian symmetry property that holds for real-valued signals, the STFT is a redundant TF transform. Therefore, the set of frequency bins {0, ..., F -1} only includes the positive frequencies.

the mathematical expectation. Moreover if the mean vector µ is also zero, the distribution is circularly symmetric. The sources are assumed mutually independent, we can thus write the distribution of the source vector introduced in (4) as:

s f n ∼ N c (0, Σ s,f n ), Σ s,f n = diag [v j,f n ] j , (7)
where diag([c m ] m ) is the diagonal matrix constructed from the coefficients c m for m = 1, ..., M .

2) NMF Source Power Parametrization: The variance v j,f n in the source model [START_REF] Vincent | Probabilistic modeling paradigms for audio source separation[END_REF] represents the short-term Power Spectral Density (PSD) of source j; it is the PSD of each frame of the short-term analysis of the source signal. The popularity of the local Gaussian model comes in particular from the possibility of introducing a parametric modeling of the source PSD. The importance of such a parametrization lies in the underdetermined nature of the source separation problem. Indeed we have to estimate the F N parameters {v j,f n } f n for each source. The idea is thus to represent the PSD of each source by a set of parameters whose cardinality is lower than F N . In this work we chose the widely used rank-reduction technique based on NMF. The idea of representing the short-term PSD of audio sources by a non-negative decomposition can be traced back to [START_REF] Benaroya | Non negative sparse representation for Wiener based source separation with a single sensor[END_REF]. The Gaussian generative model based on NMF was then introduced in [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence: With application to music analysis[END_REF] and extended to multichannel audio source separation in [START_REF] Ozerov | Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation[END_REF]. It is based on the following factorization of the matrix

V j = [v j,f n ] f n ∈ R F ×N + : V j = W j H j , (8)
where

W j ∈ R F ×Kj + and H j ∈ R Kj ×N +
. K j is the factorization rank and is generally chosen such that K j (F +N ) F N . W j is a matrix containing spectral templates while H j contains the activation of these templates over the time frames. The non-negativity constraint of this decomposition generally leads to physically meaningful spectral templates. For example they can correspond to the spectrum of the music notes in the source signal. The rows of H j then represent the activations of the notes over time frames.

3) Local Stationarity: One important hypothesis that will be used in the following consists in assuming the independence of the source STFT frames and the local stationarity of the source signals (i.e. over the support of each frame) [START_REF] Liutkus | Gaussian processes for underdetermined source separation[END_REF]. For some category of signals, particularly Gaussian and generally harmonisable α-stable [START_REF] Liutkus | Generalized Wiener filtering with fractional power spectrograms[END_REF], this local stationarity assumption implies the independence of the TF points of the STFT such that:

p({s f n } f,n ) = f,n p(s f n ). (9)
This widely used assumption is introduced for convenience because it strongly simplifies the models. However it is clearly not consistent with an STFT representation of the signals, in particular due to the framing and overlapping procedure.

C. Noise Model

As in [START_REF] Ozerov | Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation[END_REF], we assume a time-stationary and isotropic noise following a centered proper complex Gaussian distribution:

b f n ∼ N c (0, Σ b,f ), Σ b,f = σ 2 b,f I I , (10)
where σ 2 b,f > 0 and I m is the identity matrix of size m.

D. Parameter Estimation

Let X = {x f n } f,n be the set of observed data and η = {W j } j , {H j } j , {A f } f , {σ 2 b,f } f the set of parameters. In the most general case the parameters are estimated in a MAP sense:

η = arg max η ln p(η|X) = arg max η ln p(X|η) + ln p(η), (11)
where ln p(X|η) is the log-likelihood and ln p(η) the logprior over the parameters. This estimation can be done by maximizing a lower bound of the criterion [START_REF] Vincent | From blind to guided audio source separation: How models and side information can improve the separation of sound[END_REF], typically with an EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]. Let S = {s f n } f,n denote the set of latent or hidden variables. At the E-step of the EM algorithm we compute the following lower bound, from the current estimation η of the parameters:

Q MAP (η|η ) = Q ML (η|η ) + ln p(η), (12)
where Q ML (η|η ) is defined as the conditional expectation of the complete-data log-likelihood. From ( 4), ( 5) and ( 10) it writes:

Q ML (η|η ) =E S|X,η [ln p(X, S|η)] c = -N F -1 f =0 ln det(Σ b,f ) + Trace Σ -1 b,f Rxx,f -Σ -1 b,f A f RH xs,f -Σ -1 b,f Rxs,f A H f + Σ -1 b,f A f Rss,f A H f - J j=1 F -1 f =0 N -1 n=0 ln(v j,f n ) + pj,fn v j,f n , (13)
where

E S|X,η [•] is the conditional expectation of S given X, c
= denotes equality up to a constant (independent of η in ( 13)) and the statistics denoted by letters with a hat are defined at the E-step in Algorithm 1. If no prior is considered on the parameters the second term in the right-hand side of ( 11) and ( 12) disappears. The estimation is thus done in an ML sense. We see that for both MAP and ML estimations the E-step can be reduced to the computation of the conditional expectation in [START_REF] Ganseman | Evaluation of a score-informed source separation system[END_REF]. In the M-step we then maximize (13) (for ML estimation) or ( 12) (for MAP estimation) in order to obtain the new estimation of the parameters. These two steps are iterated until convergence. An alternative consists in only increasing and not maximizing the Q-function at the Mstep. In that case the algorithm is referred to as a Generalized EM (GEM). In this work we will only consider priors on the mixing filters, i.e. on the set {A f } f , so that only the update of the mixing matrix will be modified in the M-step compared with the ML case. We summarize in Algorithm 1 one iteration of the GEM algorithm. For more details about the derivation see [START_REF] Ozerov | Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation[END_REF], [START_REF] Ozerov | Multichannel nonnegative tensor factorization with structured constraints for userguided audio source separation[END_REF]. For ML estimation the update of the mixing filters at line 11 of Algorithm 1 is given by:

A f = Rxs,f R-1 ss,f . (14)
Finally, up to an additive constant independent of v j,f n , we can recognize in the last line of ( 13) the Itakura-Saito Algorithm 1: One iteration of the GEM algorithm E-step:

1: Σ s,f n = diag([v j,f n ] j ) with v j,f n = [W j H j ] f n 2: Σ x,f n = A f Σ s,f n A H f + Σ b,f 3: G s,f n = Σ s,f n A H f Σ -1 x,f n 4: ŝfn = G s,f n x f n 5: Σ post s,f n = (I J -G s,f n A f )Σ s,f n 6: Rss,fn = ŝfn ŝH f n + Σ post s,f n 7: pj,fn = [ Rss,fn ] j,j 8: Rss,f = 1 N n Rss,fn 9: Rxx,f = 1 N n x f n x H f n 10: Rxs,f = 1 N n x f n ŝH f n
M-step:

11: Update A f according to (14) (ML estimation) or Algo- rithm 2 (MAP estimation) 12: Σ b,f = Trace( Rxx,f -A f RH xs,f -Rxs,f A H f +A f Rss,f A H f )I I /I 13: W j , H j = IS-NMF( Pj ) with Pj = [p j,f n ] f n (IS) divergence [7] between v j,f n = [W j H j ] f n and the posterior mean of the source power spectrograms pj,fn = E S|X,η [|s j,f n | 2 ]
. Therefore the update of the source parameters at line 13 of Algorithm 1 is done by computing an NMF on Pj = [p j,f n ] f n ∈ R F ×N + using the IS divergence. It can be done with the standard multiplicative update rules (see [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence: With application to music analysis[END_REF]).

E. Source Reconstruction

The sources are estimated in the Minimum Mean Square Error (MMSE) sense:

ŝfn = E s f n |x f n [s f n ]. (15)
We know that s f n and b f n are two independent proper complex Gaussian random vectors. From the linearity of the normal law and (4) we can write the posterior distribution of the sources as

s f n |x f n ∼ N c (ŝ f n , Σ post s,f n )
where ŝfn and Σ post s,f n are defined at lines 4 and 5 of Algorithm 1 respectively. The sources and mixing filters are however estimated up to a frequency-dependent scale factor. This ambiguity can be fixed with a normalization strategy as in [START_REF] Ozerov | Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation[END_REF]. Nevertheless we prefer to provide as the output of the source separation system the reconstructed source images:

ŷj,fn = âj,f ŝj,fn , (16)
where ŝj,fn = [ŝ f n ] j and âj,f is the j-th column of the estimated mixing matrix Âf . This strategy implicitly solves the scale ambiguity. The time-domain signals are then reconstructed by inverse STFT.

III. REVERBERATION MODELING

The aim of this work is to propose a new way of estimating the mixing matrix A f . Based on physically motivated reverberation models, we define probabilistic priors over the frequency response of the mixing filters a ij,f . We then derive a MAP estimation of the mixing matrix. Compared with the ML Fig. 1: Simulated room impulse response using the Roomsimove toolbox [START_REF] Vincent | Roomsimove[END_REF]. The room is a 10 × 6.6 × 3 m shoebox. The reverberation time is 250 ms. The source to microphone distance is 1.38 m. The mixing time is defined in [START_REF] Duong | Spatial location priors for gaussian model based reverberant audio source separation[END_REF].

estimation given by ( 14), we take the dependencies between frequency points into account for estimating the mixing matrix. As explained in this section, these frequency correlations are induced by the specific temporal dynamics of the mixing filters, which actually are room responses. For the sake of clarity we first introduce the reverberation models by considering a single room response and by using specific notations not related to the previously introduced source separation problem. We then use these models in section IV to define priors on the frequency responses a ij,f of the mixing filters in the context of source separation. Let h(t) = h e (t) + h l (t), t = 0, ..., T -1, denote the whole RIR and h e (t), h l (t) the early and late parts respectively, having disjoint temporal supports (see Fig. 1). The room frequency response (RFR) is similarly defined for f = 0, ..., T -1 by

H(f ) = H e (f ) + H l (f ) where H (•) (f ) is the T -point Discrete Fourier Transform (DFT) of h (•) (t).
As the RIR is real-valued, the RFR satisfies the Hermitian symmetry, i.e.

H (•) (T -f ) = H (•) (f ) *
where (•) * denotes the complex conjugate. We assume that T is even and we will refer to as positive frequencies the set of indices {0, ..., T /2}.

The time instant when late reverberation starts is referred to as the mixing time. It is usually defined according to the volume V of the room in m 3 [START_REF] Lindau | Perceptual evaluation of physical predictors of the mixing time in binaural room impulse responses[END_REF]:

t 0 = C 0 √ V f s samples, ( 17
)
where C 0 = 2 × 10 -3 is a normalization constant, f s is the sampling rate in Hz and • the floor function. The time instant separating early and late reverberations in Fig. 1 is set according to this definition of the mixing time.

A. Early Reverberation 1) Model of Early Contributions: As represented in Fig. 1 the early part of the room response can be represented by a sum of early contributions. Let us consider R early contributions, each one is associated with an attenuation term ρ k and a delay τ k , k = 0, ..., R -1, such that H e (f ) ≈ G(f ) with:

G(f ) = R-1 k=0 ρ k δ f k where δ k = e -j2πτ k /T . ( 18
)
From ( 18) it follows that {G(f )} f =R,...,T -1 satisfies a recursive equation of the form (see, e.g., [START_REF] Kumaresan | On the zeros of the linear prediction-error filter for deterministic signals[END_REF]):

R r=0 ϕ e r G(f -r) = 0, (19)
such that {ϕ e r } R r=0 and {δ k } R-1 k=0 are respectively the coefficients and roots of the same polynomial of order R. Without loss of generality we assume that ϕ e 0 is equal to one. We consider that H e (f ) follows ( 19) up to a deviation κ(f ):

R r=0 ϕ e r H e (f -r) = κ(f ). (20)
κ(f ) is modeled as a centered proper complex white Gaussian noise with variance σ2 κ . From [START_REF]Autoregressive moving average modeling of late reverberation in the frequency domain[END_REF], {H e (f )} f can be seen as following an autoregressive (AR) model of order R, denoted AR(R). We can finally write the joint PDF of the sequence of H e (f ) for f = 0, ..., T /2:

p({H e (f )} f ) = p(H e (0), ..., H e (R -1)) T /2 f =R N c - R r=1 ϕ e r H e (f -r), σ 2 κ . (21)
2) Anechoic Case (R = 1): Theoretically the order R of the AR model [START_REF]Autoregressive moving average modeling of late reverberation in the frequency domain[END_REF] should be equal to the number of early contributions we consider. For example, for the RIR represented in Fig. 1 we should choose R = 13. However by considering that the direct path is dominating the early echoes, it is sufficient to choose R = 1. In this case, according to ( 18)- [START_REF]Autoregressive moving average modeling of late reverberation in the frequency domain[END_REF] we can write:

H e (f ) = δ 0 H e (f -1) + κ(f ), (22)
where δ 0 = e -j2πτ0/T and τ 0 = (r 0 /c)f s with r 0 the distance between the source and the microphone in meters and c the speed of sound in m.s -1 . This AR(1) model also expresses that, as the temporal support of the early part of the RIR is concentrated on the time instants close to zero, the associated frequency response tends to be smooth. From equation [START_REF] Ozerov | Multichannel nonnegative tensor factorization with structured constraints for userguided audio source separation[END_REF] we can write:

|H e (f )| ≈ |H e (f -1)|, (23)
arg(H e (f )) ≈ arg(H e (f -1)) -2πτ 0 /T. (24)
We represent in Fig. 2 the phase and magnitude relations between H e (f ) and H e (f -1) according to [START_REF] Parra | Convolutive blind separation of non-stationary sources[END_REF] and [START_REF] Févotte | Maximum likelihood approach for blind audio source separation using time-frequency gaussian source models[END_REF]. H e (f ) corresponds to the early part of the RIR represented in Fig. 1 (red part), it contains 13 early contributions. Nonetheless we observe that considering an order R = 1 is a reasonable approximation.

B. Late Reverberation

Late reverberation corresponds to a stage of propagation where we usually consider the sound as diffuse. It means that the sound energy is uniformly distributed in the room and over all directions [START_REF] Schultz | Diffusion in reverberation rooms[END_REF]. It corresponds to the part of the RIR where many reflections occur. Contrary to what we did for early reverberation, there are so many reflections that we cannot characterize each one individually. We conversely have to use statistical methods to describe late reverberation. 23) and [START_REF] Févotte | Maximum likelihood approach for blind audio source separation using time-frequency gaussian source models[END_REF]. Blue dots: data from the early part of the RIR represented in Fig. 1, red lines: line x = y.

1) Statistical Characterization: It is well known that for a diffuse sound field the reverberation power decays exponentially [START_REF] Jot | Analysis and synthesis of room reverberation based on a statistical time-frequency model[END_REF]. This temporal dynamics induces specific frequency correlations between the coefficients of an RFR. Schroeder in [START_REF] Schroeder | Frequency-correlation functions of frequency responses in rooms[END_REF] defined the frequency correlation functions of frequency responses in rooms. However he neglected the influence of the direct path and early echoes, resulting in a mismatch between theoretical and experimental results, as observed in [START_REF] Gustafsson | Source localization in reverberant environments: Modeling and statistical analysis[END_REF]. We showed in [START_REF]Autoregressive moving average modeling of late reverberation in the frequency domain[END_REF] that considering the validity of the exponential decay only for late reverberation yields a more accurate match between theory and experiments.

We first characterize the temporal dynamics of the late RIR by defining the Power Temporal Profile (PTP) for t = 0, ..., T -1:

hl (t) = E[|h l (t)| 2 ] = P 2 0 e -2t/τ 1 t∈{t0,...,T -1} (t), (25)
where 1 t∈T (t) is the indicator function which equals 1 if t ∈ T , 0 otherwise, P 2 0 is a constant related to the total power of late reverberation and τ is linked to the reverberation time 2 T 60 through:

τ = T 60 f s 3 ln(10) . (26)
It is important to mention that different realizations of the room responses can be interpreted as different observations at several source and microphone positions in the room. According to the theory of statistical room acoustics, {H l (f )} f is a proper centered and Wide Sense Stationary (WSS) complex Gaussian random process. We define the autocovariance function (ACVF) γ(m) and the PSD φ(t) of this process by:

γ(m) = E[H l (f )H l (f -m) * ]; (27)
φ(t) = 1 T E[|F T {H l (f )}| 2 ], (28)
where F T {•} is the T -point DFT. We have to mention that as we work in discrete time and frequency, all signals are Tperiodic. Strictly speaking the RFR is thus a T -periodic WSS random process. Moreover φ(t) has to be understood as a discretized PSD function. One can refer to [START_REF] Gu | Parametric modeling of wideband and ultrawideband channels in frequency domain[END_REF] for a review of some properties of periodic random processes. It can be shown that the PSD is related to the PTP by [START_REF]Autoregressive moving average modeling of late reverberation in the frequency domain[END_REF]:

φ(t) = T hl (T -t). (29)
From the Wiener-Khinchin theorem and ( 29) we obtain the theoretical ACVF:

γ(m) = σ 2 rev 1 -e 2/τ 1 -e 2(T -t0+1)/τ 1 -e (j2πm/T +2/τ )(T -t0+1)
1 -e j2πm/T +2/τ , (30) where σ 2 rev is defined by:

σ 2 rev = E[|H l (f )| 2 ] = C rev 1 -α παS , (31)
with α the average absorption coefficient (without dimension) and S the total wall area in m 2 . C rev = 75 is an empirical constant that has been shown to be consistent with real data in [START_REF]Autoregressive moving average modeling of late reverberation in the frequency domain[END_REF]. The average absorption coefficient can be computed using Norris-Eyring's formula [38, p. 24]:

α = 1 -e -24 ln(10)V /(cST60) . ( 32
)
2) ARMA Model: As proposed in [START_REF]Autoregressive moving average modeling of late reverberation in the frequency domain[END_REF] we consider an invertible and causal ARMA(P, Q) model for representing the late part of the RFR {H l (f )} f :

Φ(L)H l (f ) = Θ(L) (f ), (33)
where

Φ(L) = P p=0 ϕ l p L p , Θ(L) = Q q=0 θ q L q with ϕ l 0 = θ 0 = 1 and L is the lag operator, i.e. LH l (f ) = H l (f -1). (f ) is a centered proper complex white Gaussian noise of variance σ 2 for f ∈ [0, ..., T -1]
and is extended by T -periodicity elsewhere. The important point here is that the ARMA parameters can be estimated from the sole knowledge of the theoretical ACVF given by [START_REF] Lindau | Perceptual evaluation of physical predictors of the mixing time in binaural room impulse responses[END_REF], without the need of any data. Finally from this ARMA model we can write the following relationship involving the joint PDF of the sequence of H l (f ) for f = 0, ..., T /2:

p({H l (f )} f ) ∝ T /2 f =0 N c Φ(L) Θ(L) H l (f ); 0, σ 2 . ( 34
)
It is not an equality because we have omitted a term related to the Jacobian of the inverse ARMA filter transformation which does not depend on {H l (f )} f .

3) Validation: To validate the late reverberation model we verify on data the consistency of the theoretical ACVF and its ARMA parametrization. We consider a 10×6.6×3 m shoebox with reverberation time T 60 = 250 ms. We simulate 196 RIRs from the image source method using the Roomsimove toolbox [START_REF] Vincent | Roomsimove[END_REF]. According to a uniform grid of points in the room, the relative source/microphone position varies between the 196 RIRs and each RIR corresponds to one realization of the same random process. From this set of realizations we conduct a Monte-Carlo simulation to obtain the empirical ACVF for the late part of the RFR. We also compute the theoretical ACVF [START_REF] Lindau | Perceptual evaluation of physical predictors of the mixing time in binaural room impulse responses[END_REF] from the same room parameters as the ones used to simulate the RIRs. Finally we compute an ARMA(7, 2) parametrization from the theoretical ACVF (see [START_REF]Autoregressive moving average modeling of late reverberation in the frequency domain[END_REF] and [39, ch. 2] for the estimation procedure). We represent in Fig. 3 the three ACVFs: computed from the data, from expression (30), and the ARMA parametrization. We observe a good match between theory and experiments. The validity of the model is also verified on measured RIRs in [START_REF]Autoregressive moving average modeling of late reverberation in the frequency domain[END_REF].

IV. PRIORS ON THE MIXING FILTERS

We now incorporate the reverberation models introduced in section III into the source separation framework presented in section II. We consider that the mixing matrix A f can be decomposed as:

A f = A e,f + A l,f , (35)
where

A e,f = [a e ij,f ] ij ∈ C I×J and A l,f = [a l ij,f ] ij ∈ C I×J correspond
respectively to the early and late parts of the mixing filters.

A. Early Reverberation Prior 1) Prior: We consider that, independently for all i, j, the filter a e ij,f follows the AR(1) model [START_REF] Ozerov | Multichannel nonnegative tensor factorization with structured constraints for userguided audio source separation[END_REF]. According to (21) for R = 1 and by considering no prior on a e ij,0 for all i, j we can write:

ln p({A e,f } f ) = -IJ(F -1) ln(πσ 2 κ ) - 1 σ 2 κ F -1 f =1 A e,f -∆ • A e,f -1 2 F , (36)
where

∆ = [δ ij ] ij ∈ C I×J , || • || 2
F is the Frobenius norm and • is the element-wise matrix product. δ ij = e -j2πτij /La with L a the length of the mixing filters and τ ij = (r ij /c)f s where r ij is the distance between the j-th source and the i-th microphone in meters.

2) Hyperparameters: This prior is parametrized by the set of AR coefficients {δ ij } i,j and the variance σ 2 κ . The AR coefficients can be estimated within the GEM algorithm or fixed if we know the distance r ij between each source j and microphone i. In a Bayesian context the variance σ 2 κ is fixed and expresses how confident we are about the fact that a e ij,f is close to δ ij a e ij,f -1 .

B. Late Reverberation Prior 1) Prior: We consider that each filter a l ij,f independently follows for all i, j the ARMA model [START_REF] Schultz | Diffusion in reverberation rooms[END_REF]. We can also write that the filter vector a l j,f = [a l 1j,f , ..., a l Ij,f ] T ∈ C I follows:

Φ(L)a l j,f = Θ(L) f , (37)
with f ∼ N c (0, Σ ,f = σ 2 I I ). As Σ ,f is proportional to the identity matrix we do not consider any spatial correlation between the mixing filters. Finally, according to [START_REF] Jot | Analysis and synthesis of room reverberation based on a statistical time-frequency model[END_REF] we can write the following late reverberation prior:

ln p({A l,f } f ) c = F -1 f =0 -J ln det(πΣ ,f ) -Trace Φ(L) Θ(L) A l,f H Σ -1 ,f Φ(L) Θ(L) A l,f . (38)
2) Hyperparameters: This prior is parametrized by the ARMA coefficients {ϕ l p } p , {θ q } q and the variance σ 2 . We will use an ARMA(7, 2) model. The ARMA coefficients are learned from the theoretical ACVF in [START_REF] Lindau | Perceptual evaluation of physical predictors of the mixing time in binaural room impulse responses[END_REF]. This function only depends on the reverberation time, the room volume and the total wall area through the exponential decay factor τ (26), the mixing time t 0 (17), the variance σ 2 rev (31) and the absorption coefficient α [START_REF] Kumaresan | On the zeros of the linear prediction-error filter for deterministic signals[END_REF]. These parameters are assumed to be known. As for the early reverberation prior, the variance σ 2 is fixed and expresses how confident we are about the prior.

C. MAP Estimation

As explained in section II-D, only the M-step of the GEM algorithm (summarized in Algorithm 1) is modified in order to take the priors on the mixing filters into account. The MAP estimation of the set of mixing filters A e = {A e,f } f and A l = {A l,f } f at the M-step is done by minimizing:

L(A e , A l ) = -Q ML (η|η ) -ln p({A e,f } f ) -ln p({A l,f } f ), (39)
where Q ML (η|η ) is defined in [START_REF] Ganseman | Evaluation of a score-informed source separation system[END_REF] with [START_REF] Gustafsson | Source localization in reverberant environments: Modeling and statistical analysis[END_REF] and ln p({A l f } f ) in [START_REF] Naylor | Speech dereverberation[END_REF]. The whole procedure is summarized in Algorithm 2 and is explained below. The final GEM algorithm for MAP estimation corresponds to Algorithm 1 where line 11 is replaced by Algorithm 2.

A f = A e,f + A l,f , ln p({A e f } f ) is defined in
1) M-step for A e : Canceling the gradient of L(A e , A l ) with respect to A e,f leads to the updates summarized in Algorithm 2 Part 1. ⊗ denotes the Kroenecker product and vec(•) concatenates the columns of a matrix into a single column vector. We see that at each frequency the update of A e,f depends on A e,f -1 and/or A e,f +1 . This coordinate descent should thus be repeated over the whole set of frequencies until a stopping criterion is reached (e.g. convergence or maximum number of iterations).

2) Update of the AR coefficients: The AR coefficients {δ ij } i,j that parametrize the early reverberation prior (36) can be either fixed if we know the relative distance between each source and microphone, or estimated at the M-step of the GEM algorithm, from the current estimation of the early mixing filters {a e ij,f }. We choose the second option because firstly, knowing the source to microphone distance can be difficult in practice, and secondly we did not observe a significant improvement over the source separation results by setting the AR coefficients to their true values. We thus estimate in the M-step the AR coefficients that minimize the cost function [START_REF] Kay | Spectral estimation[END_REF] under the constraint |δ ij | = 1 imposed by the model. This is equivalent to maximizing the prior [START_REF] Gustafsson | Source localization in reverberant environments: Modeling and statistical analysis[END_REF] under the same constraint. It is important to note that this prior is here seen as a log-likelihood because we assume that the filters {a e ij,f } are observed from their current estimation. The update of the AR coefficients is obtained using the method of Lagrange multipliers. It is given in Part 2 of Algorithm 2.

3) M-step for A l : Unfortunately we cannot easily cancel out the gradient of L(A e , A l ) with respect to A l,f due to the MA part of the ARMA model. We thus have to consider another descent method to minimize L(A e , A l ) with respect to A l . We studied various approaches including a gradient descent with optimal step size, the Barzilai and Borwein approach [START_REF] Barzilai | Two-point step size gradient methods[END_REF] and a conjugate gradient method without or with preconditioning [START_REF] Golub | Matrix computations[END_REF]. We obtained the best results in terms of convergence speed with the third approach (with preconditioning). We detail the method here.

Our objective is to minimize L(A e , A l ) with respect to the vector of parameters of size IJF :

a l = vec A l,0 T , vec A l,1 T , ..., vec A l,F -1 T T . ( 40
)
We decompose the gradient vector g (column vector of size IJF ) as the concatenation of the F gradients:

g = [g T 0 , g T 1 , ..., g T F -1 ] T , (41)
with g f = vec 1 2 ∇ A l,f L(A e , A l ) . We can show from (39) that:

g f =N vec Σ -1 b,f (A e,f Rss,f -Rxs,f ) + N RT ss,f ⊗ Σ -1 b,f vec A l,f + Φ * (L -1 ) Θ * (L -1 ) (I J ⊗ Σ -1 ,f ) Φ(L) Θ(L) vec A l,f , (42)
where Φ * (L -1 ) = P p=0 (ϕ l p ) * L -p and Θ * (L -1 ) = Q q=0 θ * q L -q . We want to solve g = 0 which is here equivalent to solving a positive definite linear system. For that purpose we use the Preconditioned Conjugate Gradient (PCG) method [START_REF] Golub | Matrix computations[END_REF]. From [START_REF] Vinyes | MTG MASS dataset[END_REF] and the fact that Σ ,f is diagonal we can define the following preconditioning matrix:

D =      D 0 0 • • • 0 0 D 1 • • • 0 . . . . . . . . . . . . 0 0 • • • D F -1      , ( 43
)
Algorithm 2: MAP update of the mixing filters Part 1: Update {A e,f } f by coordinate descent 1: while stopping criterion not reached do

2: vec(A e,0 ) = N RT ss,0 ⊗ I I + 1 σ 2 κ (I J ⊗ Σ b,0 ) -1 × vec N Rxs,0 -N A l,0 Rss,0 + 1 σ 2 κ Σ b,0 ∆ * • A e,1
3: for all f ∈ {2, ..., F -2} do

4: vec(A e,f ) = N RT ss,f ⊗ I I + 2 σ 2 κ (I J ⊗ Σ b,f ) -1 × vec N Rxs,f -N A l,f Rss,f + 1 σ 2 κ Σ b,f ∆ • A e,f -1 + ∆ * • A e,f +1 5: end for 6: vec(A e,F -1 ) = N RT ss,F -1 ⊗ I I + 1 σ 2 κ (I J ⊗ Σ b,F -1 ) -1 × vec N Rxs,F -1 -N A l,F -1 Rss,F -1 + 1 σ 2 κ Σ b,F -1 ∆ • A e,F -2 7: end while Part 2: Update {δ ij } i,j 8: δ ij = F -1 f =1 a e ij,f (a e ij,f -1 ) * F -1 f =1 a e ij,f (a e ij,f -1 ) *
Part 3: Update {A l,f } f with the PCG method 9: Initialize g from ( 41) and ( 42) [START_REF] Vincent | BSS Eval Toolbox Version 3.0 for Matlab[END_REF] 11: while stopping criterion not reached do 41) and (42) 16: g p = D -1 g where D is given by ( 43) and ( 44)

10: Initialize ω = [ω T 0 , ω T 1 , ..., ω T F -1 ] T with ω f = D -1 f g f and D f given by
12: γ = [γ T 0 , γ T 1 , ..., γ T F -1 ] T with γ f = N ( RT ss,f ⊗Σ -1 b,f )ω f + Φ * (L -1 ) Θ * (L -1 ) (I J ⊗Σ -1 ,f ) Φ(L) Θ(L) ω f 13: µ = (ω H g)/(ω H γ) 14: a l ← a l -µω 15: Compute g from (
17: α = -(γ H g p )/(ω H γ) 18: ω ← g p + αω 19: end while Part 4: Update {A f } f 20: A f = A e,f + A l,f
where the diagonal matrix D f is given by:

D f = N diag( RT ss,f ⊗ Σ -1 b,f ) + 1 σ 2 I IJ N ψ s=0 |ψ s | 2 . ( 44
)
The parameters {ψ s } N ψ s=0 are the coefficients of the polynomial Ψ(L) which approximates the transfer function of the inverse ARMA model Φ(L)/Θ(L). For the chosen ARMA(7, 2) model, we obtain a very accurate approximation by setting N ψ = 2048. We summarize the PCG method for the update of {A l,f } f in Part 3 of Algorithm 2.

V. EXPERIMENTS

We present in this section the dataset and the evaluation criteria used for the experiments. We discuss the short mixing filters assumption used for writing the source image model [START_REF] Tachibana | Melody line estimation in homophonic music audio signals based on temporal-variability of melodic source[END_REF]. We also describe our initialization strategy before running the GEM algorithm. We explain how the algorithm parameters are chosen and we finally detail the source separation results. Audio examples are available from our demo web page 3 .

A. Dataset

Our experiments were conducted from the audio tracks without effects provided by the Musical Audio Signal Separation (MASS) dataset [START_REF] Vinyes | MTG MASS dataset[END_REF]. We created 8 stereo mixtures by simulating mixing filters with the Roomsimove toolbox [START_REF] Vincent | Roomsimove[END_REF]. The room was a 4.45 × 3.55 × 2.5 m shoebox with a reverberation time of 128 ms. The microphone spacing was set to 1 m, as the distance between the source and the center of the microphone pair. The direction of arrival (DoA) and the type of the sources are described in Table I for each mixture, along with the duration in seconds. As the MASS dataset provides stereo sources, each one is first converted to mono, downsampled to 16 kHz and filtered with the associated RIRs to create a source image. We finally sum all the source images to create a mixture.

B. Evaluation Criteria

We evaluate the source separation performance in terms of reconstructed source images. We use standard energy ratios: the Signal-to-Distortion Ratio (SDR), Signal-to-Interference Ratio (SIR), Signal-to-Artifact Ratio (SAR) and source Imageto-Spatial distortion Ratio (ISR). These criteria expressed in decibels (dB) are defined in [START_REF] Vincent | The Signal Separation Evaluation Campaign (2007-2010): Achievements and Remaining Challenges[END_REF]. We used the BSS Eval Toolbox available at [START_REF] Vincent | BSS Eval Toolbox Version 3.0 for Matlab[END_REF] to compute these measures.

C. Short Mixing Filters Assumption

The aim of this section is to discuss the short mixing filters assumption used to express the convolutive mixture in the STFT domain. We assume that the mixing filters length is equal to the reverberation time. On the one hand, for rewriting the time-domain convolution as a multiplication in the STFT measures setting average scores SDR A f is blindly estimated 0.0 A e,f is oracle and A l,f is zero 0.9 A f is oracle 7.9 domain as in (3), we theoretically have to use an STFT analysis window longer than the reverberation time [START_REF] Parra | Convolutive blind separation of non-stationary sources[END_REF]. On the other hand, we cannot choose a too long analysis window because of the local stationarity assumption on the source signals. Moreover, the true mixing filters are still tractable in the limit case where the analysis window length equals the reverberation time. For example, considering the dataset presented above, we can choose a 128 ms long window and think of the mixing matrix A f in (4) as containing the frequency response of the mixing filters. Even if this approach is theoretically arguable, it can be shown to be effective in terms of oracle source separation results 4 . This may explain why this mixing model is still widely used in audio source separation. Indeed, even if it is not accurate for representing the true mixing process (time-domain convolution), its simplicity is useful for blind audio source separation and it allows us to achieve satisfactory oracle performance.

ISR
To validate this argument, we conducted a preliminary source separation experiment with the baseline method summarized in Algorithm 1 (for the ML estimation case). The STFT analysis window length is equal to the reverberation time (128 ms). We consider three settings: (1) A f is blindly estimated in the M-step according to equation ( 14); (2) the early part A e,f is fixed to oracle values (using the true mixing filters truncated at the mixing time) and the late part A l,f is fixed to zero; (3) the whole mixing matrix A f = A e,f + A l,f is fixed to oracle values. All other parameters are blindly estimated. The results averaged on the 29 sources of the dataset are shown in Table II. We first observe that even if the length of the STFT analysis window equals the length of the mixing filters, we achieve good separation results when A f is oracle. We also observe that when we neglect the late reverberation by setting A l,f to zero, the separation performance drastically decreases compared with the case where the whole mixing matrix is oracle. This shows the importance of achieving a good estimation for both the early and late parts of the mixing filters.

The conclusion of this preliminary experiment is the following one: if we were able to estimate the true mixing filters, which correspond to room responses, we would achieve good source separation performance, even if the mixing filters length equals the STFT analysis window length. The aim of this work is precisely to incorporate constraints on the estimation of the mixing filters to go towards this objective.

D. Blind Initialization

The GEM algorithm is very sensitive to the parameter initialization. In order to obtain satisfactory separation results, we have to provide a "good initialization". As in [START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF], [START_REF]Multichannel audio source separation with probabilistic reverberation modeling[END_REF] we first tried the initialization procedure based on a hierarchical clustering algorithm to estimate the mixing filters [START_REF] Winter | MAP-based underdetermined blind source separation of convolutive mixtures by hierarchical clustering and 1-norm minimization[END_REF] and a permutation solving method [START_REF] Sawada | Grouping separated frequency components by estimating propagation model parameters in frequency-domain blind source separation[END_REF]. From the estimated mixing filters the sources can be recovered via projection of the mixture over the source directions and binary masking in the TF plan. The NMF parameters can then be estimated from the separated sources. However with this approach we obtained less satisfactory initialization results than with a more empirical method inspired from the one described in [21, section IV.H]. The initialization procedure we used is restricted to stereo mixtures and is detailed in the Appendix.

E. Algorithm Parameters 1) STFT parameters:

We use half-overlapping sine windows. The STFTs are computed using the MATLAB code provided in the context of the SiSEC challenge at [START_REF]Signal Separation Evaluation Campaign (SiSEC 2008)[END_REF]. According to the discussion in section V-C, we fix the length of the analysis window to 128 ms or 2048 points at a sampling rate of 16 kHz.

2) NMF order: We arbitrarily choose an NMF of order 10 for all the sources.

3) Stopping criterion for the descent algorithms: It was sufficient in terms of convergence to repeat the coordinate descent for the update of {A e,f } f for 2 iterations. The stopping criterion for the PCG method used in the update of {A l,f } f is set to 20 iterations.

4) Prior hyperparameters: The AR coefficients for the early reverberation prior are estimated in the M-step of the GEM algorithm. The ARMA(7, 2) coefficients for the late reverberation prior are learned from the theoretical ACVF [START_REF] Lindau | Perceptual evaluation of physical predictors of the mixing time in binaural room impulse responses[END_REF] which is defined according to some known room parameters. The variances of both priors are critical hyperparameters because they have a strong influence on the separation results. Indeed, we see in Fig. 4 that the average SDR for the mixture 4 of the dataset varies from 0.2 to 4.2 dB according to the value of the variances. These variances control the strength of the priors, the smallest they are, the strongest the priors will be and conversely. In particular, if their values are too high the priors have no effect on the results compared with an ML estimation of the mixing filters; their estimation will indeed be dominated by the contribution from the log-likelihood.

F. Source Separation Results

For comparing ML and MAP estimations we run 500 iterations of the GEM algorithms from the same initialization. As can be seen in Fig. 4, the source separation results for MAP estimation are highly dependent on the prior variances.

Moreover the results we observe in Fig. 4 are considerably different if we take another mixture from our dataset. This justifies the use of a grid-search on the variances to assess the performance of MAP estimation. 1) Results with best variances on the grid: For all the variances in the grid (σ 2 κ , σ 2 ) ∈ {10 0 , 10 -1 , ..., 10 -9 } × {10 0 , 10 -1 , ..., 10 -9 } we first compute the source separation performance that leads to the best average SDR individually for each mixture. ML and MAP results with this strategy are shown on average for each mixture in Table III, in columns ML and MAP BEST respectively. The mean results in the last line of Table III are computed by averaging over all the sources of the dataset. We see that over all the mixtures, the SDR improvement due to MAP estimation ranges from 0.1 to 3.6 dB. On average, MAP estimation leads to an improvement of 0.9 dB of ISR, 1.3 dB of SIR, 1.1 dB of SAR and 1.2 dB of SDR. We also represent in Fig. 5a a scatter plot of the SDR for the 29 sources of our dataset. On y and x axes we represent respectively the MAP and ML results. Red circles (resp. blue stars) above (resp. below) the line x = y represent better separation results with MAP (resp. ML) estimation. We observe that globally the improvement due to MAP estimation is significant as it leads to an increase of the SDR on almost all the sources. When it is not the case, the decrease is small as blue stars are very close to the line x = y. Those results show that for each mixture there exist values for the variances that lead to an improvement of the source separation performance.

2) Results with variances chosen by cross-validation: We now present the results where the variances are chosen by cross-validation on the same grid as before. For each fold we consider 7 mixtures out of 8 as a training set and we keep the last one for testing. The variances are chosen as the ones that lead to the best average SDR over all the sources in the training set. Interestingly, for every fold the chosen variances are (σ 2 κ , σ 2 ) = (10 -6 , 10 -7 ). The MAP results for this cross-validation are shown in column MAP CV of Table III. Compared with ML estimation we observe a decrease of the SDR on mixtures 1 and 8 by 1.5 and 1.3 dB respectively. However for all the other mixtures, the SDR is increased by taking the priors into account. The improvement ranges from 0.1 to 3.4 dB. On average MAP estimation with crossvalidation leads to an improvement of 0.7 dB of ISR, 0.6 dB of SIR, 0.2 dB of SAR and 0.7 dB of SDR. We observe from the scatter plot in Fig. 5b that with MAP estimation the quality of the separated sources never degrades much, but on some sources the improvement is important. Indeed, red circles tend to be further from line x = y than blue stars.

VI. CONCLUSION

In this paper we introduced two probabilistic priors for constraining the estimation of the mixing filters in multichannel audio source separation. The priors aim to capture the temporal characteristics of the mixing filters by modeling the frequency correlations of their frequency responses. Experiments showed that these new priors lead to better source separation results than a standard approach with unconstrained mixing filters. However the results are highly dependent on the variances of the priors. In future works we could use a conjugate prior on the variances so that their values are constrained to be close to the ones we obtained with cross-validation, but not exactly equal.

This work targeted a semi-informed scenario where the reverberation time, the volume and the total wall area of the room were assumed to be known to define the late reverberation prior. In future works we could try to blindly estimate the reverberation time (see, e.g, [START_REF] Ratnam | Blind estimation of reverberation time[END_REF]). Moreover it could be interesting to conduct a sensitivity analysis of the method with respect to errors on these room parameters.

Future works will also include the extension of our model to non-punctual sources using spatial covariance matrices. These matrices can be factorized as the outer product of sub-sources mixing matrices [START_REF] Ozerov | A general flexible framework for the handling of prior information in audio source separation[END_REF] on which we could consider the same priors as the ones presented in this paper. The source with a SDR around -20 dB corresponds to the brushes in the mixture 1 of the dataset. This source is not accurately modeled by an NMF.

Finally we could consider not only modeling the frequency correlations of the mixing filters but also the spatial correlations. As used in [START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF], [START_REF] Duong | Spatial location priors for gaussian model based reverberant audio source separation[END_REF] for stereo audio source separation, the spatial correlation functions are theoretically defined according to the distance between the microphones. For that purpose we could consider a full noise covariance matrix in equation [START_REF] Gu | Parametric modeling of wideband and ultrawideband channels in frequency domain[END_REF] or use a vector ARMA model to represent late reverberation in the frequency domain.

APPENDIX BLIND INITIALIZATION FOR STEREO MIXTURES

The initialization procedure we used is as follows: 1) Stack left and right mixture STFTs so as to create a 2F × N complex-valued matrix. 2) Compute a K init -component IS-NMF and split the resulting matrix containing the frequency atoms in two matrices associated to left and right channels. K init is chosen between 3 and 5 depending on the mixture. 3) Reconstruct K init left (i = 1) and right (i = 2) components c ik,f n by Wiener filtering (see, e.g., [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence: With application to music analysis[END_REF]).

By considering an anechoic model we can write c ik,f n = ρ ik e -j2πf τ ik /La |c k,f n |e j arg(c k,f n ) . 4) Compute cik,fn = c ik,f n /e j arg(c 1k,f n ) such that we can write c1k,fn = ρ 1k |c k,f n | and c2k,fn = ρ 2k |c k,f n |e -j2πf (τ 2k -τ 1k )/La . 5) Compute cik,f = 1 N n cik,fn . 6) Compute and unwrap ξ k,f = arg( c2k,f c1k,f ) such that we can write ξ k,f = -2πf (τ 2k -τ 1k )/L a . 7) Assuming punctual and spatially disjoint sources, ξ k,f forms J clusters as multiple components are spatially associated with the same source. Use the K-means algorithm to cluster the ξ k,f . 8) By denoting K j the set of indices k associated to source j, initialize the mixing filters as a ij,f = 1 #Kj k∈Kj cik,f . 9) Compute a pre-separation by projection of the mixture over the source direction and TF masking. We used the MATLAB code provided in the context of the SiSEC challenge at [START_REF]Signal Separation Evaluation Campaign (SiSEC 2008)[END_REF]. 10) Initialize the NMF parameters from the separated sources using the Kullback-Leibler divergence. individually. For that purpose we simply split the initial mixing filters in the time domain according to the mixing time in [START_REF] Duong | Spatial location priors for gaussian model based reverberant audio source separation[END_REF].

Fig. 2 :

 2 Fig. 2: Phase (left figure) and magnitude (right figure) relations between H e (f ) and H e (f -1) illustrating (23) and[START_REF] Févotte | Maximum likelihood approach for blind audio source separation using time-frequency gaussian source models[END_REF]. Blue dots: data from the early part of the RIR represented in Fig.1, red lines: line x = y.

Fig. 3 :

 3 Fig. 3: Real (top figure) and imaginary (bottom figure) parts of the RFR ACVFs for late reverberation. Blue solid line: empirical ACVF, red dash-dot line: theoretical ACVF from (30), green dashed line: ARMA(7, 2) parametrization.

Fig. 4 :

 4 Fig. 4: Average SDR over a grid-search on the prior variances for mixture 4 of the dataset.

Fig. 5 :

 5 Fig.5: Scatter plot of the SDR (in dB) for the 29 sources of the dataset. MAP and ML results are on y and x axis respectively. The source with a SDR around -20 dB corresponds to the brushes in the mixture 1 of the dataset. This source is not accurately modeled by an NMF.

11 )

 11 Initialize σ 2 b,f = 10 -2 i,n |x i,f n | 2 /(IN ). MAP estimation requires the initialization of a e ij,f and a l ij,f

TABLE I :

 I Description of the dataset: duration (in seconds) and source properties, i.e. instrument and DoA (in degrees), for each mixture.

TABLE II :

 II Preliminary source separation results (in dB) averaged over the 29 sources of the dataset.

TABLE III :

 III Source separation results in dB: without prior (ML), with prior and best variances on the grid-search (MAP BEST), with prior and variances chosen by cross-validation (MAP CV). The mean is computed over the 29 sources of the dataset.

The reverberation time in seconds is defined as the time it takes for the sound energy to decrease by 60 dB after extinction of the source.

http://perso.telecom-paristech.fr/leglaive/demoMASSwPRP.html

Moreover, one could wonder if the residual noise in the mixture model (4) does not become overwhelmingly large due to the short mixing filters assumption which is not satisfied. Experimental results which are not detailed here show that most of the mixture power is represented by the estimated source images and the residual noise is negligible.

ACKNOWLEDGMENT

The authors would like to thank the anonymous reviewers for their valuable comments.


Alaa Saade

Florent Krzakala

Marc Lelarge

Lenka Zdeborová

Clustering from sparse pairwise measurements

come

I. INTRODUCTION

A. Problem and model

Similarity-based clustering is a standard approach to label items in a dataset based on some measure of their resemblance. In general, given a dataset {x i } i∈[n] ∈ X n , and a symmetric measurement function s : X 2 → R quantifying the similarity between two items, the aim is to cluster the dataset from the knowledge of the pairwise measurements s ij := s(x i , x j ), for 1 ≤ i < j ≤ n. This information is conveniently encoded in a similarity graph, which vertices represent items in the dataset, and the weighted edges carry the pairwise similarities. Typical choices for this similarity graph are the complete graph and the nearest neighbor graph (see e.g. [START_REF] Luxburg | A tutorial on spectral clustering[END_REF] for a discussion in the context of spectral clustering).

Here however, we will not assume the measurement function s to quantify the similarity between items, but more generally ask that the measurements be typically different depending on the cluster memberships of the items, in a way that will be made quantitative in the following. For instance, s could be a distance in an Euclidean space or could take values in a set of colors (i.e. s does not need to be realvalued). Additionally, we will not assume knowledge of the measurements for all pairs of items in the dataset, but only for O(n) of them chosen uniformly at random. Sampling is a well-known technique to speed up computations by reducing the number of non-zero entries [START_REF] Achlioptas | Fast computation of low rank matrix approximations[END_REF]. The main challenge is to choose the lowest possible sampling rate while still being able to detect the signal of interest. In this paper, we compute explicitly this fundamental limit for a simple probabilistic model and present three algorithms allowing partial recovery of the signal above this limit. Below the limit, in the case of two clusters, no algorithm can give an output positively correlated with the true clusters. Our three algorithms are respectively a belief propagation algorithm and two spectral algorithms based on the non-backtracking operator and the Bethe Hessian. Although these three algorithms are intimately related, so far, a sketch of rigorous analysis is available only for the spectral properties of the non-backtracking matrix. From a practical perspective however, belief propagation and the Bethe Hessian are much simpler to implement and show even better numerical performance.

To evaluate the performance of our proposed algorithms, we construct a model with n items in k predefined clusters of same average size n/k, by assigning to each item i ∈ [n] a cluster label c i ∈ [k] with uniform probability 1/k. We assume that the pairwise measurement between an item in cluster a and another item in cluster b is a random variable with density p a,b . We choose the observed pairwise measurements uniformly at random, by generating an Erdős-Rényi random graph G = (V = [n], E) ∈ G(n, α/n). The average degree α corresponds to the sampling rate: pairwise measurements are observed only on the edges of G, and α therefore controls the difficulty of the problem. From the base graph G, we build a measurement graph by weighting each edge (ij) ∈ E with the measurement s ij , drawn from the probability density p ci,cj . The aim is to recover the cluster assignments c i for i ∈ [n] from the measurement graph thus constructed.

We consider the sparse regime α = O(1), and the limit n → ∞ with fixed number of clusters k. With high probability, the graph G is disconnected, so that exact recovery of the clusters, as considered e.g. in [START_REF] Jog | Information-theoretic bounds for exact recovery in weighted stochastic block models using the renyi divergence[END_REF], [START_REF] Chen | Information recovery from pairwise measurements: A shannon-theoretic approach[END_REF], is impossible. In this paper, we address instead the question of how many measurements are needed to partially recover the cluster assignments, i.e. to infer cluster assignments ĉi such that the following quantity, called overlap, is strictly positive:

max σ∈S k 1 n i 1(σ(ĉ i ) = c i ) -1 k 1 -1 k , (1)
where S k is the set of permutations of [k]. This quantity is monotonously increasing with the number of correctly classified items. In the limit n → ∞, it vanishes for a random guess, and equals unity if the recovery is perfect. Finally, we note an important special case for which analytical results can be derived, which is the case of symmetric clusters:

∀a, b ∈ [k] p a,b (s) := 1(a = b)p in (s) + 1(a = b)p out (s), (2)
where p in (s) (resp. p out (s)) is the probability density of observing a measurement s between items of the same cluster (resp. different clusters). For this particular case, we conjecture that all of the three algorithms we propose achieve partial recovery of the clusters whenever α > α c , where

1 α c = 1 k K ds p in (s) -p out (s) 2 p in (s) + (k -1)p out (s) , (3)
where K is the support of the function p in + (k -1)p out . This expression corresponds to the threshold of a related reconstruction problem on trees [START_REF] Mézard | Reconstruction on trees and spin glass transition[END_REF]. In the following, we substantiate this claim for the case of k = 2 symmetric clusters, and discrete measurement distributions. Note that the model we introduce is a special case of the labeled stochastic block model of [START_REF] Heimlicher | Community detection in the labelled stochastic block model[END_REF]. In particular, for the case k = 2, it was proven in [START_REF] Lelarge | Reconstruction in the labeled stochastic block model[END_REF] that partial recovery is information theoretically impossible if α < α c . In this contribution, we argue that this bound is tight, namely that partial recovery is possible whenever α > α c , and that the algorithms we propose are optimal, in that they achieve this threshold. Note also that the symmetric model ( 2) contains the censored block model of [START_REF] Abbe | Decoding binary node labels from censored edge measurements: Phase transition and efficient recovery[END_REF]. More precisely, if p in and p out are discrete distributions on {±1} with p in (+1) = p out (-1) = 1 -, then α c = (1 -2 ) -2 . In this case, the claimed result is known [START_REF] Saade | Spectral detection in the censored block model[END_REF], and to the best of our knowledge, this is the only case where our result is known.

B. Motivation and related work

The ability to cluster data from as few pairwise comparisons as possible is of broad practical interest [START_REF] Chen | Information recovery from pairwise measurements: A shannon-theoretic approach[END_REF]. First, there are situations where all the pairwise comparisons are simply not available. This is particularly the case if a comparison is the result of a human-based experiment. For instance, in crowdclustering [START_REF] Gomes | Crowdclustering[END_REF], [START_REF] Yi | Crowdclustering with sparse pairwise labels: A matrix completion approach[END_REF], people are asked to compare a subset of the items in a dataset, and the aim is to cluster the whole dataset based on these comparisons. Clearly, for a large dataset of size n, we can't expect to have all O(n 2 ) measurements. Second, even if these comparisons can be automated, the typical cost of computing all pairwise measurements is O(n 2 d) where d is the dimension of the data. For large datasets with n in the millions or billions, or large dimensional data, like high resolution images, this cost is often prohibitive. Storing all O(n 2 ) measurements is also problematic. Our work supports the idea that if the measurements between different classes of items are sufficiently different, a random subsampling of O(n) measurements might be enough to accurately cluster the data. This work is inspired by recent progress in the problem of detecting communities in the sparse stochastic block model (SBM) where partial recovery is possible only when the average degree α is larger than a threshold value, first conjectured in [START_REF] Decelle | Asymptotic analysis of the stochastic block model for modular networks and its algorithmic applications[END_REF], and proved in [START_REF] Mossel | Stochastic block models and reconstruction[END_REF]- [START_REF] Mossel | A proof of the block model threshold conjecture[END_REF]. A belief propagation (BP) algorithm similar to the one presented here is introduced in [START_REF] Decelle | Asymptotic analysis of the stochastic block model for modular networks and its algorithmic applications[END_REF], and argued to be optimal in the SBM. Spectral algorithms that match the performance of BP were later introduced in [START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF], [START_REF] Saade | Spectral clustering of graphs with the bethe hessian[END_REF]. The spectral algorithms presented here are based on a generalization of the operators that they introduce.

C. Outline and main results

In Sec. II, we describe three closely related algorithms to solve the partial recovery problem of Sec. I-A. The first one is a belief propagation (BP) algorithm approximating the Bayes optimal solution. The other two are spectral methods derived from BP. We show numerically that all three methods achieve the threshold [START_REF] Jog | Information-theoretic bounds for exact recovery in weighted stochastic block models using the renyi divergence[END_REF]. Next in Sec. III we substantiate this claim for the spectral method based on the non-backtracking operator.

II. ALGORITHMS

A. Belief propagation

We consider a measurement graph generated from the model of Section I-A. From Bayes' rule, we have:

P({c i }|{s ij }) = 1 Z (ij)∈E p ci,cj (s ij ) , (4)
where Z is a normalization. The Bayes optimal assignment, maximizing the overlap (1), is ĉi = argmax P i , the mode of the marginal of node i. We approximate this marginal using belief propagation (BP):

P i (c i ) ≈ 1 Z i l∈∂i k c l =1 p ci,c l (s il )P l→i (c l ) , (5)
where ∂i denotes the neighbors of node i in the measurement graph G, Z i is a normalization, and the P i→j (c i ) are the fixed point of the recursion:

P i→j (c i ) = 1 Z i→j l∈∂i\j k c l =1 p ci,c l (s il )P l→i (c l ) . (6)
In practice, starting from a random initial condition, we iterate (6) until convergence, and estimate the marginals from (5). On sparse tree-like random graphs generated by our model, BP is widely believed to give asymptotically accurate results, though a rigorous proof is still lacking. This algorithm is general and applies to any model parameters p ab . For now on, however, we restrict our theoretical discussion to the symmetric model [START_REF] Achlioptas | Fast computation of low rank matrix approximations[END_REF]. Eq. ( 6) can be written in the compact form P = F (P), where P ∈ R 2mk and m = |E| is the number of edges in G.

The first step in understanding the behavior of BP is to note that in the case of symmetric clusters (2), there exists a trivial fixed point of the recursion (6), namely P i→j (c i ) = 1/k. This fixed point is uninformative, yielding a vanishing overlap. If this fixed point is stable, then starting from an initial condition close to it will cause BP to fail to recover the clusters. We therefore investigate the linearization of (6) around this fixed point, given by the Jacobian J F .

B. The non-backtracking operator

A simple computation yields

J F = B ⊗ I k - 1 k U k , (7)
where I k is the k × k identity matrix, U k is the k × k matrix with all its entries equal to 1, ⊗ denotes the tensor product, and B is a 2m × 2m matrix called the non-backtracking operator, acting on the directed edges of the graph G, with elements for (ab) and (cd) ∈ E:

B (a→b),(c→d) = w(s cd )1(a = d)1(b = c) , ∀s, w(s) := p in (s) -p out (s) p in (s) + (k -1)p out (s) . (8)
Note that to be consistent with the analysis of BP, our definition of the non-backtracking operator is the transpose of the definition of [START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF]. This matrix generalizes the nonbacktracking operators of [START_REF] Saade | Spectral detection in the censored block model[END_REF], [START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF] to arbitrary edge weights. More precisely, for the censored block model [START_REF] Abbe | Decoding binary node labels from censored edge measurements: Phase transition and efficient recovery[END_REF], we have s = ±1 and w(s) = (1 -2 )s so that B is simply a scaled version of the matrix introduced in [START_REF] Saade | Spectral detection in the censored block model[END_REF]. We also introduce an operator C ∈ R n×2m defined as

C i,j→l = w(s jl )1(i = l) . (9)
This operator follows from the linearization of eq. ( 5) for small P l→i . Based on these operators, we propose the following spectral algorithm. First, compute the real eigenvalues of B with modulus greater than 1. Let r be their number, and denote by v 1 , ..., v r ∈ R 2m the corresponding eigenvectors. If r = 0, raise an error. Otherwise, form the matrix

Y = [v 1 • • • v r ] ∈ R 2m×r
by stacking the eigenvectors in columns, and let X = CY ∈ R n×r . Finally, regarding each item i as a vector in R r specified by the i-th line of X, cluster the items, using e.g. the k-means algorithm.

Theoretical guarantees for the case of k = 2 clusters are sketched in the next section, stating that this simple algorithm succeeds in partially recovering the true clusters all the way down to the transition (3). Intuitively, this algorithm can be thought of as a spectral relaxation of belief propagation. Indeed, for the particular case of k = 2 symmetric clusters, we will argue that the spectral radius of B is larger than 1 if and only if α > α c . As a simple consequence, whenever α < α c , the trivial fixed point of BP is stable, and BP fails to recover the clusters. On the other hand, when α > α c , a small perturbation of the trivial fixed point grows when iterating BP. Our spectral algorithm approximates the evolution of this perturbation by replacing the non-linear operator F by its Jacobian J F . In practice, as shown on figure 1, the non-linearity of the function F allows BP to achieve a better overlap than the spectral method based on B, but a rigorous proof that BP is asymptotically optimal is still lacking.

C. The Bethe Hessian

The non-backtracking operator B of the last section is a large, non-symmetric matrix, making the implementation of the previous algorithm numerically challenging. A much smaller, closely related symmetric matrix can be defined that empirically performs as well in recovering the clusters, and in fact slightly better than B. For a real parameter x ≥ 1, define a matrix H(x) ∈ R n×n with non-zero elements:

H ij (x) =    1 + l∈∂i w(s il ) 2 x 2 -w(s il ) 2 if i = j - xw(sij ) x 2 -w(sij ) 2 if (ij) ∈ E , ( 10
)
where ∂i denotes the set of neighbors of node i in the graph G, and w is defined in [START_REF] Abbe | Decoding binary node labels from censored edge measurements: Phase transition and efficient recovery[END_REF]. A simple computation, analogous to [START_REF] Saade | Spectral detection in the censored block model[END_REF], allows to show that (λ ≥ 1, v) is an eigenpair of B, if and only H(λ)v = 0. This property justifies the following picture [START_REF] Saade | Spectral clustering of graphs with the bethe hessian[END_REF]. For x large enough, H(x) is positive definite and has no negative eigenvalue. As we decrease x, H(x) gains a new negative eigenvalue whenever x becomes smaller than an eigenvalue of B. Finally, at x = 1, there is a one to one correspondence between the negative eigenvalues of H(x) and the real eigenvalues of B that are larger than 1. We call Bethe Hessian the matrix H( 1), and propose the following spectral algorithm, by analogy with Sec. II-B. First, compute all the negative eigenvalues of H(1). Let r be their number. If r = 0, raise an error. Otherwise, denoting v i , ..., v r ∈ R n the corresponding eigenvectors, form the matrix

X = [v 1 • • • v r ] ∈
R n×r by stacking them in columns. Finally, regarding each item i as a vector in R r specified by the i-th line of X, cluster the items, using e.g. the k-means algorithm.

In the case of two symmetric clusters, the results of the next section imply that if α > α c , denoting by λ 1 > 1 the largest eigenvalue of B, the smallest eigenvalue of H(λ 1 ) is 0, and the corresponding eigenvector allows partial recovery of the clusters. While the present algorithm replaces the matrix H(λ 1 ) by the matrix H(1) and is therefore beyond the scope of this theoretical guarantee, we find empirically that the eigenvectors with negative eigenvalues of H(1) are also positively correlated with the hidden clusters, and in fact allow better recovery (see figure 1), without the need to build the non-backtracking operator B and to compute its leading eigenvalue. This last algorithm also has an intuitive justification. It is well known [START_REF] Yedidia | Bethe free energy, kikuchi approximations, and belief propagation algorithms[END_REF] that BP tries to optimize the so-called Bethe free energy. In the same way B can be seen as a spectral relaxation of BP, H(1) can be seen as a spectral relaxation of the direct optimization of the Bethe free energy. In fact, it corresponds to the Hessian of the Bethe free energy around a trivial stationary point (see e.g. [START_REF] Saade | Spectral clustering of graphs with the bethe hessian[END_REF], [START_REF] Saade | Matrix completion from fewer entries: Spectral detectability and rank estimation[END_REF]).

D. Numerical results

Figure 1 shows the performance of all three algorithms on model-generated problems. We consider the symmetric problem defined by (2) with k = 2, 3, fixed p in and p out , chosen to be Gaussian with a strong overlap, and we vary α. All three algorithms achieve the theoretical threshold. While all the algorithms presented in this work assume the knowledge of the parameters of the model, namely the functions p a,b for a, b ∈ [k], we argue that the belief propagation Fig. 2. Clustering of toy datasets using belief propagation. Each dataset is composed of 20000 points, 200 of which come labeled and constitute the training set. We used the Euclidean distance as the measurement function s, and estimated the probability densities p ab on the training set using kernel density estimation (middle row). Although these estimates are very noisy and overlapping, belief propagation is able to achieve a very high accuracy using a random measurement graph G of small average degree α (top row). For comparison, we show in the third row the result of spectral clustering with the normalized Laplacian, using a 3-nearest neighbors similarity graph (see e.g. [START_REF] Luxburg | A tutorial on spectral clustering[END_REF]) built from G, i.e. using only the available measurements.

algorithm is robust to large imprecisions on the estimation of these parameters. To support this claim, we show on figure 2 the result of the belief propagation algorithm on standard toy datasets where the parameters were estimated on a small fraction of labeled data.

III. PROPERTIES OF THE NON-BACKTRACKING OPERATOR

We now state our claims concerning the spectrum of B. We restrict ourselves to the case where k = 2 and p in and p out are distributions on a finite alphabet.

Claim 1: Consider an Erdős-Rényi random graph on n vertices with average degree α, with variables assigned to vertices c i ∈ {1, 2} uniformly at random independently from the graph and measurements s ij between any two neighboring vertices drawn according to the probability density: p ci,cj (s) = 1(c i = c j )p in (s) + 1(c i = c j )p out (s) for two fixed (i.e. independent of n) discrete distributions p in = p out on S. Let B be the matrix defined by ( 8) and denote by

|λ 1 | ≥ |λ 2 | ≥ • • • ≥ |λ 2m |
the eigenvalues of B in order of decreasing magnitude, where m is the number of edges in the graph. Recall that α c is defined by (3). Then, with probability tending to 1 as n → ∞:

(i) If α < α c , then |λ 1 | ≤ α αc + o(1). (ii) If α > α c , then λ 1 = α αc + o(1) and |λ 2 | ≤ α αc + o(1).
Additionally, denoting by v the eigenvector associated with λ 1 , Cv is positively correlated with the planted variables (c i ) i∈ [n] , where C is defined in [START_REF] Saade | Spectral detection in the censored block model[END_REF].

Note that for the censored block model, our claim implies Theorem 1 in [START_REF] Saade | Spectral detection in the censored block model[END_REF]. The main idea which substantiates our claim is to introduce a new non-backtracking operator with spectral properties close to those of B and then apply the techniques developed in [START_REF] Bordenave | Non-backtracking spectrum of random graphs: community detection and non-regular ramanujan graphs[END_REF] to it. We try to use notations consistent with [START_REF] Bordenave | Non-backtracking spectrum of random graphs: community detection and non-regular ramanujan graphs[END_REF]: for an oriented edge e = u → v = (u, v) from node u to node v, we set e 1 = u, e 2 = v and e -1 = (v, u). For a matrice M , its transpose is denoted by M * . We also define

σ i = 2c i -3 for each i ∈ [n].
We start by a simple transformation: if t is the vector in R E defined by t e = σ e1 and is the Hadamard product, i.e. (t x) e = σ e1 x e , then we have

B * x = λx ⇔ B X (t x) = λ(t x) , (11)
with B X defined by B X ef = B f e σ e1 σ e2 . In particular, B X and B have the same spectrum and there is a trivial relation between their eigenvectors. With X e = σ e1 w(s e )σ e2 , we have:

B X ef = X e 1(e 2 = f 1 )1(e 1 = f 2 )
. Moreover, note that the random variables (A f = σ f1 σ f2 ) f ∈E are random signs with P(A f = 1) = 1/2 and the random variables (X f ) f ∈E are such that [START_REF] Bordenave | Non-backtracking spectrum of random graphs: community detection and non-regular ramanujan graphs[END_REF]. More precisely, we define P the linear mapping on R E defined by (P x) e = Y e x e -1 (i.e. in matrix form P ef = Y e 1(f = e -1 )). Note that P * = P and since Y 2 e = 1, P is an involution so that P is an orthogonal matrix. A simple computation shows that (B Y ) k P = P (B Y ) * k , hence (B Y ) k P is a symmetric matrix. This symmetry corresponds to the oriented path symmetry in [START_REF] Bordenave | Non-backtracking spectrum of random graphs: community detection and non-regular ramanujan graphs[END_REF] and is crucial to our analysis. If (τ j,k ), 1 ≤ j ≤ 2m are the eigenvalues of (B Y ) k P and (x j,k ) is an orthonormal basis of eigenvectors, we deduce that

EX f = EX 2 f = 1
(B Y ) k = 2m j=1 τ j,k x j,k (P x j,k ) * . ( 12
)
Since P is an orthogonal matrix (P x j,k ), 1 ≤ j ≤ 2m is also an orthonormal basis of R E . In particular, [START_REF] Decelle | Asymptotic analysis of the stochastic block model for modular networks and its algorithmic applications[END_REF] gives the singular value decomposition of (B Y ) k . Indeed, if t j,k = |τ j,k | and y j,k = sign(τ j,k )P x j,k , then we get (B Y ) k = 2m j=1 t j,k x j,k y * j,k , which is precisely the singular value decomposition of (B Y ) k . As shown in [START_REF] Bordenave | Non-backtracking spectrum of random graphs: community detection and non-regular ramanujan graphs[END_REF], for large k, the decomposition (12) carries structural information on the graph.

A crucial element in the proof of [START_REF] Bordenave | Non-backtracking spectrum of random graphs: community detection and non-regular ramanujan graphs[END_REF] is the result of Kesten and Stigum [START_REF] Kesten | A limit theorem for multidimensional galton-watson processes[END_REF], [START_REF]Additional limit theorems for indecomposable multidimensional galton-watson processes[END_REF] and we give now its extension required here which can be seen as a version of Kesten and Stigum's

Fig. 1 .

 1 Fig.1. Performance in clustering model-generated measurement graphs in the symmetric case[START_REF] Achlioptas | Fast computation of low rank matrix approximations[END_REF]. The overlap is averaged over 20 realizations of graphs of size n = 10 5 , with k = 2, 3 clusters, and Gaussian p in , pout with mean respectively 1.5 and 0, and unit variance. The theoretical transition (3) is at αc ≈ 2.63 for k = 2, and conjectured to be at αc ≈ 5.5 for k = 3. All three methods achieve the theoretical transition, although the Bethe Hessian (H) and belief propagation (BP) achieve a higher overlap than the non-backtracking operator (B).

2 s

 2 (p in (s) -p out (s))2 p in (s) + p out (s) = 1 α c .We now define another non-backtracking operator B Y . First, letting(s)= pout(s) pin(s)+pout(s) ∈ [0, 1], we define the sequence of independent random variables { Ỹe } e∈E with P( Ỹe = +1|s e ) = 1 -P( Ỹe = -1|s e ) = 1 -(s e ), so that E[ Ỹe |s e ] = w(s e ). We define Y e = Ỹe σ e1 σ e2 and finally B Y ef = Y f 1(e 2 = f 1 )1(e 1 = f 2 ) , so that E[B Y |G, {s e } e∈E ] = B X . It turns out that the analysis of the matrix B Y can be done with the techniques developped in

ACKNOWLEDGMENT

This work has been supported by the ERC under the European Union's FP7 Grant Agreement 307087-SPARCS and by the French Agence Nationale de la Recherche under reference ANR-11-JS02-005-01 (GAP project).

work in a random environment. We write N * = {1, 2, . . . } and U = ∪ n≥0 (N * ) n the set of finite sequences composed by N * , where (N * ) 0 contains the null sequence ∅. For u, v ∈ U , we note |u| = n for the lenght of u and uv for the sequence obtained by the juxtaposition of u and v. Suppose that {(N u , A u1 , A u2 , . . . )} u∈U is a sequence of i.i.d. random variables with value in N × R N * such that N u is a Poisson random variable with mean α and the A ui are independent i.i.d. random signs with P(A u1 = 1) = P(A u1 = -1) = 1 2 . We then define the following random variables: first s u such that P(s u |A u = 1) = p in (s u ) and P(s u |A u = -1) = p out (s u ), then X u = A u w(s u ) and Y u = A u Ỹu where P( Ỹu = 1|s u ) = 1 -P(Y u = -1|s u ) = 1 -(s u ). We assume that for all u ∈ U and i > N u , A ui = s ui = X ui = Y ui = 0. N u will be the number of children of node u and the sequence (s u1 , . . . , s uNu ) the measurements on edges between u and its children. We set for u = u 1 u 2 . . . u n ∈ U ,

. . Y u1...un . We define (with the convention 0 0 = 0) ,

Then conditionnaly on the variables (s u ) u∈U , M t is a martingale converging almost surely and in L 2 as soon as α > α c . The fact that this martingale is bounded in L 2 follows from an argument given in the proof of Theorem 3 in [START_REF] Heimlicher | Community detection in the labelled stochastic block model[END_REF].

In order to apply the technique of [START_REF] Bordenave | Non-backtracking spectrum of random graphs: community detection and non-regular ramanujan graphs[END_REF], we need to deal with the -th power of the non-backtracking operators. For not too large, the local struture of the graph (up to depth ) can be coupled to a Poisson Galton-Watson branching process, so that the computations done for M above provide a good approximation of the -th power of the non-backtracking operator and we can use the algebraic tools about perturbation of eigenvalues and eigenvectors, see the Bauer-Fike theorem in Section 4 in [START_REF] Bordenave | Non-backtracking spectrum of random graphs: community detection and non-regular ramanujan graphs[END_REF].

IV. CONCLUSION

We have considered the problem of clustering a dataset from as few measurements as possible. On a reasonable model, we have made a precise prediction on the number of measurements needed to cluster better than chance, and have substantiated this prediction on an interesting particular case. We have also introduced three efficient and optimal algorithms, based on belief propagation, to cluster modelgenerated data starting from this transition. Our results suggest that clustering can be significantly sped up by using a number of measurements linear in the size of the dataset, instead of quadratic. These algorithms, however, require an estimate of the distribution of the measurements between objects depending on their cluster membership. On toy datasets, we have demonstrated the robustness of the belief propagation algorithm to large imprecisions on these estimates, paving the way for broad applications in real world, large-scale data analysis. A natural avenue for future work is the unsupervised estimation of these distributions, through e.g. an expectationmaximization approach.


Laetitia Canabady-Rochelle
email: laetitia.canabady-rochelle@univ-lorraine.fr

Cédric Paris

Katalin Selmeczi

Caroline Gaucher

Igor Clarot

Pierre Leroy

Halima Alem

Raphaël Schneider

Patrick Chaimbault

Laurence Muhr

Sandrine Boschi-Muller

Screening and Purification of Metal Chelating Peptides in Hydrolysates for their Antioxidant Properties

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction Reactive Oxygen Species (ROS) induced oxidation

Ageing and various pathologies in human related to oxidative stress, in particular inflammatory bowel diseases with iron anemia complication Degradation of consumption products (food, cosmetic matrix)

Living organisms developped their own antioxidant defenses against these ROS

Imbalance

Oxidative stress

Importance to discover natural antioxidants and evaluation of their bioactivity

An interdisciplinary approach

Some results

SPR using a M 2+ -chip:

-a miniaturized approach of IMAC-type affinity chromatography -screening of the most interesting sources of metal-chelating peptides.

Positive correlations evidenced between SPR screening and metal chelation tests using UV-Vis spectrometry on five sources of peptides.

The affinity parameters are determined by SPR on some synthetic model peptides and then used to predict the separation of peptides on IMAC chromatographic column.

Novel chromatography-MS method currently developed:

-Screening of metal-chelating peptides sequences in a given complex mixture -Methodology successfully set-up on a mixture of synthetic peptides -Further tests on real hydrolysates Identified peptides currently assessed for their antioxidant activity by chemical metal chelation and by in cellulo tests.

Potential health application

-Treatment of ferriprive anemia related to MICI, or inversely to hematochromatosis -Anchoring of these chelating peptides on nanoparticles for bioimaging or therapeutic applications.

Metal-chelating peptides obtained after protein hydrolysis -Antioxidant properties related to the inhibition of Fenton reaction, responsible of ROS production -Used in pharmaceutics, foods and cosmetic Novel strategy for the screening and purification of the peptides to target their antioxidant properties. Novel approach developed in order to save time. Mass shift of +56 u and isotopic pattern = ironcontaining ion. Bioactivity tests

Metal chelating peptide

Health applications Coordination chemistry

Nanoparticules for bioimaging

Processes and Modelling

17th International conference on Oxidative Stress, Reduction and Redox Homeostasis Paris, June 13-15, 2016 Empirical approach Fractioning SEPARATION Complex mixture of peptides

	PROTEINES PROTEINES		HYDROLYSAT HYDROLYSAT		HYDROLYSAT
	VEGETALES VEGETALES		INTERMEDIAIRE INTERMEDIAIRE	FINAL
	H 2 N H 2 N	COOH COOH	H 2 N H 2 N	COOH COOH	H 2 N H 2 N	COOH
			Peptide bioactif	Peptide bioactif
			intermédiaire