Simon Leglaive

Umut Şimşekli

Antoine Liutkus

Laurent Girin

Radu Horaud

Speech enhancement with variational autoencoders and alpha-stable distributions

Keywords: Speech enhancement, variational autoencoders, alpha-stable distribution, Monte Carlo expectation-maximization

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

Speech enhancement is one of the central problems in audio signal processing [START_REF] Loizou | Speech enhancement: theory and practice[END_REF]. The goal is to recover a clean speech signal after observing a noisy mixture. In this work, we address single-channel speech enhancement, which can be seen as an under-determined source separation problem, where the sources are of different nature.

One popular statistical approach for source separation combines a local Gaussian model of the time-frequency signal coefficients with a variance model [START_REF] Vincent | Probabilistic modeling paradigms for audio source separation[END_REF]. In this framework, non-negative matrix factorization (NMF) techniques have been used to model the timefrequency-dependent signal variance [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence: With application to music analysis[END_REF][START_REF] Mohammadiha | Supervised and unsupervised speech enhancement using nonnegative matrix factorization[END_REF]. Recently, discriminative approaches based on deep neural networks (DNNs) have also been investigated for speech enhancement, with the aim of estimating either clean spectrograms or time-frequency masks, given noisy spectrograms [START_REF] Wang | Supervised speech separation based on deep learning: An overview[END_REF][START_REF] Xu | A regression approach to speech enhancement based on deep neural networks[END_REF][START_REF] Weninger | Speech enhancement with LSTM recurrent neural networks and its application to noise-robust ASR[END_REF]]. As a representative example, a DNN is used in [START_REF] Xu | A regression approach to speech enhancement based on deep neural networks[END_REF] to map noisy speech log-power spectrograms into clean speech log-power spectrograms.

Even more recently, generative models based on deep learning, and in particular variational autoencoders (VAEs) [START_REF] Kingma | Auto-encoding variational Bayes[END_REF], have been used for single-channel [START_REF] Bando | Statistical speech enhancement based on probabilistic integration of variational autoencoder and non-negative matrix factorization[END_REF][START_REF] Leglaive | A variance modeling framework based on variational autoencoders for speech enhancement[END_REF] and multi-channel speech enhancement [START_REF] Sekiguchi | Bayesian multichannel speech enhancement with a deep speech prior[END_REF][START_REF] Leglaive | Semi-supervised multichannel speech enhancement with variational autoencoders and non-negative matrix factorization[END_REF]. These generative model-based approaches provide important advantages and justify the interest of semi-supervised methods for speech enhancement. Indeed, as shown in [START_REF] Bando | Statistical speech enhancement based on probabilistic integration of variational autoencoder and non-negative matrix factorization[END_REF][START_REF] Leglaive | A variance modeling framework based on variational autoencoders for speech enhancement[END_REF], fully-supervised methods such as [START_REF] Xu | A regression approach to speech enhancement based on deep neural networks[END_REF] may have issues for generalizing to unseen noise types. The method proposed in [START_REF] Leglaive | A variance modeling framework based on variational autoencoders for speech enhancement[END_REF] was shown to outperform both a semi-supervised NMF baseline and the fully-supervised deep learning approach [START_REF] Xu | A regression approach to speech enhancement based on deep neural networks[END_REF].

In most cases, probabilistic models for source separation or speech enhancement rely on a Gaussianity assumption, which turns This work is supported by the ERC Advanced Grant VHIA #340113.

out to be restrictive for audio signals [START_REF] Liutkus | Generalized Wiener filtering with fractional power spectrograms[END_REF]. As a result, heavytailed distributions have started receiving increasing attention in the audio processing community [START_REF] Şimşekli | Alpha-stable matrix factorization[END_REF][START_REF] Liutkus | Cauchy nonnegative matrix factorization[END_REF][START_REF] Yoshii | Student's t nonnegative matrix factorization and positive semidefinite tensor factorization for singlechannel audio source separation[END_REF]. In particular, α-stable distributions (cf. Section 3) are becoming popular heavy-tailed models for audio modeling due to their nice theoretical properties [START_REF] Kuruoglu | Signal processing in alpha-stable noise environments: a least lp-norm approach[END_REF][START_REF] Liutkus | Generalized Wiener filtering with fractional power spectrograms[END_REF][START_REF] Şimşekli | Alpha-stable matrix factorization[END_REF][START_REF] Leglaive | Alpha-stable multichannel audio source separation[END_REF][START_REF] Fontaine | Explaining the parameterized Wiener filter with alpha-stable processes[END_REF][START_REF] Fontaine | Multichannel Audio Modeling with Elliptically Stable Tensor Decomposition[END_REF][START_REF] Şimşekli | Alpha-stable low-rank plus residual decomposition for speech enhancement[END_REF].

In this work, we investigate the combination of a deep learningbased generative speech model with a heavy-tailed α-stable noise model. The rationale for introducing a noise model based on heavytailed distributions as opposed to a structured NMF approach as in [START_REF] Leglaive | A variance modeling framework based on variational autoencoders for speech enhancement[END_REF] is to avoid relying on restricting assumptions regarding stationarity or temporal redundancy of the noisy environment, that may be violated in practice, leading to errors in the estimates. In addition, we let the noise model remain unsupervised in order to avoid the aforementioned generalization issues regarding the noisy recording environment. We develop a Monte Carlo expectation-maximization algorithm [START_REF] Wei | A Monte Carlo implementation of the EM algorithm and the poor man's data augmentation algorithms[END_REF] for performing maximum likelihood estimation at test time. Experiments performed under challenging conditions show that the proposed approach outperforms the competing approaches in terms of both perceptual quality and intelligibility.

SPEECH MODEL

We work in the short-term Fourier transform (STFT) domain where B = {0, ..., F -1}×{0, ..., N -1} denotes the set of time-frequency bins. For (f, n) ∈ B, f denotes the frequency index and n the timeframe index. We use s f n , b f n , x f n ∈ C to denote the complex STFT coefficients of the speech, noise, and mixture signals, respectively.

As in [START_REF] Bando | Statistical speech enhancement based on probabilistic integration of variational autoencoder and non-negative matrix factorization[END_REF][START_REF] Leglaive | A variance modeling framework based on variational autoencoders for speech enhancement[END_REF], independently for all (f, n) ∈ B, we consider the following generative speech model involving a latent random vector hn ∈ R L , with L ≪ F : hn ∼ N (0, I);

(1)

s f n hn ∼ Nc(0, σ 2 s,f (hn)), (2)
where N (x; µ, Σ) denotes the multivariate Gaussian distribution for a real-valued random vector, I is the identity matrix of appropriate size, and Nc(x; µ, σ 2 ) denotes the univariate complex proper Gaussian distribution. As represented in Fig. 1a,

{σ 2 s,f ∶ R L ↦ R+} F -1
f =0 is a set of non-linear functions corresponding to a neural network which takes as input hn ∈ R L . This variance term can be understood as a model for the short-term power spectral density of speech [START_REF] Liutkus | Gaussian processes for underdetermined source separation[END_REF]. We denote by θs the parameters of this generative neural network.

An important contribution of VAEs [START_REF] Kingma | Auto-encoding variational Bayes[END_REF] is to provide an efficient way of learning the parameters of such a generative model. Let s = {sn ∈ C F } N -1 n=0 be a training dataset of clean-speech STFT time frames and h = {hn ∈ R L } N -1 n=0 the set of associated latent random vectors. Taking ideas from variational inference, VAEs estimate the parameters θs by maximizing a lower bound of the log-likelihood ln p(s; θs) defined by:

L (θs, ψ)=E q(h s;ψ) [ln p (s h; θs)]-DKL (q (h s; ψ) ∥ p(h)) , (3)
where q (h s; ψ) denotes an approximation of the intractable true posterior distribution p(h s; θs), and DKL(q ∥ p) = Eq[ln(q p)] is the Kullback-Leibler divergence. Independently for all the dimensions l ∈ {0, ..., L -1} and all the time frames n ∈ {0, ..., N -1}, q(h s; ψ) is defined by:

h l,n sn ∼ N μl sn ⊙2 , σ2 l sn ⊙2 , (4)
where h l,n = (hn) l and (⋅) ⊙⋅ denotes element-wise exponentiation. As represented in Figure 1b,

{μ l ∶ R F + ↦ R} L-1 l=0 and {σ 2 l ∶ R F + ↦ R+} L-1 l=0
are non-linear functions corresponding to a neural network which takes as input the speech power spectrum at a given time frame. ψ denotes the parameters of this recognition network, which also have to be estimated by maximizing the variational lower bound defined in [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence: With application to music analysis[END_REF]. Using (1), ( 2) and ( 4) we can develop this objective function as follows:

L (θs, ψ) c = - F -1 f =0 N -1 n=0 E q(hn sn;ψ) dIS s f n 2 ; σ 2 f (hn) + 1 2 L-1 l=0 N -1 n=0 ln σ2 l sn ⊙2 -μl sn ⊙2 2 -σ2 l sn ⊙2 , (5)
where dIS(x; y) = x yln(x y) -1 is the Itakura-Saito divergence. Finally, using the so-called "reparametrization trick" [START_REF] Kingma | Auto-encoding variational Bayes[END_REF] to approximate the intractable expectation in (5), we obtain an objective function which is differentiable with respect to both θs and ψ and can be optimized using gradient-ascent-based algorithms. It is important to note that the only reason why the recognition network is introduced is to learn the parameters of the generative network.

NOISE AND MIXTURE MODELS

In the previous section we have seen how to learn the parameters of the generative model ( 1)-( 2). This model can then be used as a speech signal probabilistic prior for a variety of applications. In this paper we are interested in single-channel speech enhancement. We do not assume prior knowledge about the recording environment, so that the noise model remains unsupervised. Independently for all (f, n) ∈ B, the STFT coefficients of noise are modeled as complex circularly symmetric α-stable random variables [START_REF] Samorodnitsky | Stable non-Gaussian random processes: stochastic models with infinite variance[END_REF]:

b f n ∼ SαS(σ b,f ), (6)
where α ∈]0, 2] is the characteristic exponent and σ b,f ∈ R+ is the scale parameter. As proposed in [START_REF] Fontaine | Multichannel Audio Modeling with Elliptically Stable Tensor Decomposition[END_REF] for multichannel speech enhancement, this scale parameter is only frequency-dependent, it does not depend on time. For algorithmic purposes, the noise model ( 6) can be conveniently rewritten in an equivalent scale mixture of Gaussians form [START_REF] Andrews | Scale mixtures of normal distributions[END_REF], by making use of the product property of the symmetric α-stable distribution [START_REF] Samorodnitsky | Stable non-Gaussian random processes: stochastic models with infinite variance[END_REF]:

φ f n ∼ P α 2 S 2 cos(πα 4) 2 α ; (7) b f n φ f n ∼ Nc(0, φ f n σ 2 b,f ), (8)
where φ f n ∈ R+ is called the impulse variable. It locally modulates the variance of the conditional distribution of b f n given in [START_REF] Kingma | Auto-encoding variational Bayes[END_REF]. P α 2 S denotes a positive stable distribution of characteristic exponent α 2. It corresponds to a right-skewed heavy-tailed distribution defined for non-negative random variables [START_REF] Magron | Lévy NMF for robust nonnegative source separation[END_REF]. These impulse variables can be understood as carrying uncertainty about the stationary noise assumption made in the marginal model [START_REF] Xu | A regression approach to speech enhancement based on deep neural networks[END_REF], where the scale parameter does not depend on the time-frame index.

The observed mixture signal is modeled as follows for all (f, n) ∈ B:

x f n = √ gns f n + b f n , (9)
where gn ∈ R+ represents a frame-dependent but frequencyindependent gain. The importance of this parameter was experimentally shown in [START_REF] Leglaive | A variance modeling framework based on variational autoencoders for speech enhancement[END_REF]. We further consider the conditional independence of the speech and noise STFT coefficients so that:

x f n hn, φ f n ∼ Nc 0, gnσ 2 s,f (hn) + φ f n σ 2 b,f . (10)

INFERENCE

Let θu = g = {gn ∈ R+} N -1 n=0 , σ 2 b = {σ 2 b,f ∈ R+} F -1 f =0
be the set of model parameters to be estimated. For maximum likelihood estimation, in this section we develop a Monte-Carlo expectation maximization (MCEM) algorithm [START_REF] Wei | A Monte Carlo implementation of the EM algorithm and the poor man's data augmentation algorithms[END_REF], which iteratively applies the so-called E-and M-steps until convergence, which we detail below. Remember that the speech generative model parameters θs have been learned during a training phase (see Section 2). We denote by x = {x f n } (f,n)∈B the set of observed data while z = hn,

φ n = {φ f n } F -1 f =0 N -1
n=0 is the set of latent variables. We will also use xn = {x f n } F -1 f =0 and zn = {hn, φ n } to respectively denote the set of observed and latent variables at a given time frame n.

Monte Carlo E-Step.

Let θ ⋆ u be the current (or the initial) value of the model parameters. At the E-step of a standard expectation-maximization algorithm, we would compute the following conditional expectation of the complete-data log-likelihood

Q(θu; θ ⋆ u ) = E p(z x;θs,θ ⋆ u ) [ln p(x, z; θs, θu)].
However, this expectation cannot be here computed analytically. We therefore approximate Q(θu; θ ⋆ u ) using an empirical average:

Q(θu; θ ⋆ u ) c = - 1 R R r=1 (f,n)∈B ln gnσ 2 s,f h (r) n + φ (r) f n σ 2 b,f + x f n 2 gnσ 2 s,f h (r) n + φ (r) f n σ 2 b,f -1 , (11)
where c = denotes equality up to an additive constant, and z

(r) n = h (r) n , φ (r) n = φ (r) f n F -1
f =0 , r ∈ {1, ..., R}, is a sample drawn from the posterior p(zn xn; θs, θ ⋆ u ) using a Markov Chain Monte Carlo (MCMC) method. This approach forms the basis of the MCEM algorithm [START_REF] Wei | A Monte Carlo implementation of the EM algorithm and the poor man's data augmentation algorithms[END_REF]. Note that unlike the standard EM algorithm, it does not ensure an improvement in the likelihood at each iteration. Nevertheless, some convergence results in terms of stationary points of the likelihood can be obtained under suitable conditions [START_REF] Chan | Monte Carlo EM estimation for time series models involving counts[END_REF].

In this work we use a (block) Gibbs sampling algorithm [START_REF] Robert | Monte Carlo Statistical Methods[END_REF]. From an initialization z (0) n , it consists in iteratively sampling from the so-called full conditionals. More precisely, at the m-th iteration of the algorithm and independently for all n ∈ {0, ..., N -1}, we first sample h

(m) n ∼ p hn xn, φ (m-1) n ; θs, θ ⋆ u . Then, we sample φ (m) f n ∼ p φ f n x f n , h (m) n ; θs, θ ⋆
u independently for all f ∈ {0, ..., F -1}. Those two full conditionals are unfortunately analytically intractable, but we can use one iteration of the Metropolis-Hastings algorithm to sample from them. This approach corresponds to the Metropolis-within-Gibbs sampling algorithm [28, p. 393]. One iteration of this method is detailed in Algorithm 1. The proposal distributions for hn and φ f n are respectively given in lines 2 and 6. The two acceptance probabilities required in lines 3 and 7 are computed as follows:

α (h) n = min ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ 1, p hn F -1 ∏ f =0 p x f n hn, φ (m-1) f n ; θs, θ ⋆ u p h (m-1) n F -1 ∏ f =0 p x f n h (m-1) n , φ (m-1) f n ; θs, θ ⋆ u ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ ; (12)
α (φ) n = min ⎛ ⎝ 1, p x f n h (m) n , φfn ; θs, θ ⋆ u p x f n h (m) n , φ (m-1) f n ; θs, θ ⋆ u ⎞ ⎠ . (13)
The two distributions involved in the computation of those acceptance probabilities are defined in ( 1) and [START_REF] Leglaive | A variance modeling framework based on variational autoencoders for speech enhancement[END_REF]. Finally, we only keep the last R samples for computing Q(θu; θ ⋆ u ), i.e. we discard the samples drawn during a so-called burn-in period.

M-

Step. At the M-step we want to minimize -Q(θu; θ ⋆ u ) with respect to θu. Let C(θu) denote the cost function associated with this non-convex optimization problem. Similar to [START_REF] Févotte | Algorithms for nonnegative matrix factorization with the β-divergence[END_REF], we adopt a majorization-minimization approach. Let us introduce the two following sets of auxiliary variables c = {c (r) f n ∈ R+} r,f,n and λ = {λ (r) k,f n ∈ R+} k,r,f,n . We can show using standard concave/convex inequalities (see e.g. [START_REF] Févotte | Algorithms for nonnegative matrix factorization with the β-divergence[END_REF]) that C(θu) ≤ G(θu, c, λ), where:

G(θu, c, λ) = 1 R R r=1 (f,n)∈B ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ ln c (r) f n + 1 c (r) f n gnσ 2 s,f h (r) n + φ (r) f n σ 2 b,f -c (r) f n + x f n 2 ⎛ ⎝ λ (r) 1,f n 2 gnσ 2 s,f h (r) n + λ (r) 2,f n 2 φ (r) f n σ 2 b,f ⎞ ⎠ ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ . ( 14
)
Moreover, this upper bound is tight, i.e. C(θu) = G(θu, c, λ), for c (r)

f n = gnσ 2 s,f h (r) n + φ (r) f n σ 2 b,f ; (15)
λ (r) 1,f n = gnσ 2 s,f h (r) n gnσ 2 s,f h (r) n + φ (r) f n σ 2 b,f -1 ; (16)
λ (r) 2,f n = φ (r) f n σ 2 b,f gnσ 2 s,f h (r) n + φ (r) f n σ 2 b,f -1 . ( 17
)
Minimizing G(θu, c, λ) with respect to the model parameters θu is a convex optimization problem. By zeroing the partial derivatives of G with respect to each scalar in θu we obtain update rules that depend on the auxiliary variables. We then replace these auxiliary variables with the formulas given in ( 15)-( 17), which makes the upper Algorithm 1 m-th iteration of the Metropolis-within-Gibbs sampling algorithm 1: independently for all n ∈ {0, ..., N -1} do 2:

Sample hn ∼ N (h

(m-1) n , 2 I) 3:
Compute acceptance probability α (h) n (equation ( 12))

4:

Set h

(m) n = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ hn if α (h) n > u (h) n ∼ U([0, 1]) h (m-1) n otherwise 5:
independently for all f ∈ {0, ..., F -1} do Compute acceptance probability α

(φ)
n (equation ( 13))

8:

Set φ (m)

f n = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ φfn if α (φ) n > u (φ) n ∼ U([0, 1]) φ (m-1) f n otherwise 9:
end for 10: end for bound G tight. This procedure ensures that the cost C(θu) decreases [START_REF] Hunter | A tutorial on MM algorithms[END_REF]. The resulting final updates are given as follows:

σ 2 b,f ← σ 2 b,f ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ N -1 ∑ n=0 x f n 2 R ∑ r=1 φ (r) f n v (r) x,f n -2 N -1 ∑ n=0 R ∑ r=1 φ (r) f n v (r) x,f n -1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ 1 2 ; ( 18
)
gn ← gn ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ F -1 ∑ f =0 x f n 2 R ∑ r=1 σ 2 s,f h (r) n v (r) x,f n -2 F -1 ∑ f =0 R ∑ r=1 σ 2 s,f h (r) n v (r) x,f n -1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ 1 2 , ( 19
)
where we introduced v (r)

x,f n = gnσ 2 s,f h (r) n + φ (r)
f n σ 2 b,f in order to ease the notations. Non-negativity is ensured provided that the parameters are initialized with non-negative values.

Speech reconstruction.

Once the unsupervised model parameters θu are estimated with the MCEM algorithm, we need to estimate the clean speech signal. For all (f, n) ∈ B, let sfn = √ gns f n be the scaled version of the speech STFT coefficients. We estimate these variables according to their posterior mean, given by:

ŝfn = E p(s f n x f n ;θu,θs) [s f n ] = E p(zn xn;θu,θs) E p(s f n zn,xn;θu,θs) [s f n ] = E p(zn xn;θu,θs) gnσ 2 s,f (hn) gnσ 2 s,f (hn) + φ f n σ 2 b,f x f n . (20)
As before, this expectation cannot be computed analytically so it is approximated using the Metropolis-within-Gibbs sampling algorithm detailed in Algorithm 1. This estimate corresponds to a probabilistic Wiener filtering averaged over all possible realizations of the latent variables according to their posterior distribution.

EXPERIMENTS

Reference method: The proposed approach is compared with the recent method [START_REF] Leglaive | A variance modeling framework based on variational autoencoders for speech enhancement[END_REF]. The speech signal in this paper is modeled in the exact same manner as in the current work, only the noise model differs. This latter is a Gaussian model with an NMF parametrization of the variance [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence: With application to music analysis[END_REF]: 0 .8 5 0 .8 5 0 .8 6 0 .8 7 0 .8 7 0 .8 9 0 .8 9 0 .9 0 0 .9 0 0 .9 1 0 .9 1 Fig. 2: Speech enhancement results obtained with the proposed method as a function of the characteristic exponent α in the noise model [START_REF] Xu | A regression approach to speech enhancement based on deep neural networks[END_REF]. α = 2.0 actually corresponds to α = 1.999. The value of the median is indicated within each boxplot. and H b are estimated from the noisy mixture signal. This method also relies on an MCEM algorithm. By comparing the proposed method with [START_REF] Leglaive | A variance modeling framework based on variational autoencoders for speech enhancement[END_REF], we fairly investigate which noise model leads to the best speech enhancement results. Note that the method proposed in [START_REF] Leglaive | A variance modeling framework based on variational autoencoders for speech enhancement[END_REF] was shown to outperform both a semi-supervised NMF baseline and the fully-supervised deep learning approach [START_REF] Xu | A regression approach to speech enhancement based on deep neural networks[END_REF], the latter having difficulties for generalizing to unseen noise types. We do not include here the results obtained with these two other methods.

b f n ∼ Nc 0, (W b H b ) f,n , where W b ∈ R F ×K
Database: The supervised speech model parameters in the proposed and the reference methods are learned from the training set of the TIMIT database [START_REF] Garofolo | TIMIT acoustic phonetic continuous speech corpus[END_REF]. It contains almost 4 hours of 16-kHz speech signals, distributed over 462 speakers. For the evaluation of the speech enhancement algorithms, we mixed clean speech signals from the TIMIT test set and noise signals from the DEMAND database [START_REF] Thiemann | The Diverse Environments Multi-channel Acoustic Noise Database (DEMAND): A database of multichannel environmental noise recordings[END_REF], corresponding to various noisy environments: domestic, nature, office, indoor public spaces, street and transportation. We created 168 mixtures at a 0 dB signal-to-noise ratio (one mixture per speaker in the TIMIT test set). Note that both speakers and sentences are different than in the training set.

Parameter setting: The STFT is computed using a 64-ms sine window (i.e. F = 513) with 75%-overlap. Based on [START_REF] Leglaive | A variance modeling framework based on variational autoencoders for speech enhancement[END_REF], the latent dimension in the speech generative model ( 1)-( 2) is fixed to L = 64. In this reference method, the NMF rank of the noise model is fixed to K = 10. The NMF parameters are randomly initialized. For both the proposed and the reference method, the gain gn is initialized to one for all time frames n. For the proposed method, the noise scale parameter σ b,f is also initialized to one for all frequency bins. We run 200 iterations of the MCEM algorithm. At each Monte-Carlo E-Step, we run 40 iterations of the Metropolis-within-Gibbs algorithm and we discard the first 30 samples as the burn-in period. The parameter2 in line 2 of Algorithm 1 is set to 0.01.

Neural network: The structure of the generative and recognition networks is the same as in [START_REF] Leglaive | A variance modeling framework based on variational autoencoders for speech enhancement[END_REF] and is represented in Fig. 1. Hidden layers use hyperbolic tangent (tanh(⋅)) activation functions and output layers use identity activation functions (I d (⋅)). The output of these last layers is therefore real-valued, which is the reason why we consider logarithm of variances. For learning the parameters θs and φ, we use the Adam optimizer [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] with a step size of 10 -3 , exponential decay rates for the first and second moment estimates of 0.9 and 0.999 respectively, and an epsilon of 10 -7 for preventing division by zero. 20% of the TIMIT training set is kept as a validation set, and early stopping with a patience of 10 epochs is used. Weights are initialized using the uniform initializer described in [START_REF] Glorot | Understanding the difficulty of training deep feedforward neural networks[END_REF].

Results: The estimated speech signal quality is evaluated in terms of standard energy ratios expressed in decibels (dBs) [START_REF] Vincent | Performance measurement in blind audio source separation[END_REF]: the signal-to-distortion (SDR), signal-to-interference (SIR) and signalto-artifact (SAR) ratios. We also consider the perceptual evaluation of speech quality (PESQ) [START_REF] Rix | Perceptual evaluation of speech quality (PESQ)-a new method for speech quality assessment of telephone networks and codecs[END_REF] measure (in [-0.5, 4.5]), and the short-time objective intelligibility (STOI) measure [START_REF] Taal | An algorithm for intelligibility prediction of time-frequency weighted noisy speech[END_REF] (in [0, 1]). For all measures, the higher the better. We first study the performance of the proposed method according to the choice of the characteristic exponent α in the noise model [START_REF] Xu | A regression approach to speech enhancement based on deep neural networks[END_REF]. Results presented in Fig. 2 indicate that according to the PESQ and STOI measures, the best performance is obtained for α = 2 (Gaussian case). 1 The SDR indicates that we should choose α = 1.8. Indeed, for greater values of α, the SIR starts to decrease (∼1dB difference between α = 1.8 and α = 2) while the SAR remains stable (results are not shown here due to space constraints). Therefore, in Fig. 3 we compare the results obtained with the reference [START_REF] Leglaive | A variance modeling framework based on variational autoencoders for speech enhancement[END_REF] and the proposed method using α = 1.8. With the proposed method, the estimated speech signal contains more interferences (SIR is lower) but less artifacts (SAR is higher). According to the SDR both methods are equivalent. But for intelligibility and perceptual quality, artifacts are actually more disturbing than interferences [START_REF] Venkataramani | Performance based cost functions for end-to-end speech separation[END_REF], which is the reason why the proposed method obtains better results in terms of both STOI and PESQ measures. For reproducibility, a Python implementation of our algorithm and audio examples are available online. 2

CONCLUSION

In this work, we proposed a speech enhancement method exploiting a speech model based on VAEs and a noise model based on alphastable distributions. At the expense of more interferences, the proposed α-stable noise model reduces the amount of artifacts in the estimated speech signal, compared to the use of a Gaussian NMFbased noise model as in [START_REF] Leglaive | A variance modeling framework based on variational autoencoders for speech enhancement[END_REF]. Overall, it is shown that the proposed approach improves the intelligibility and perceptual quality of the enhanced speech signal. Future works include extending the proposed approach to a multi-microphone setting using multivariate αstable distributions [START_REF] Leglaive | Alpha-stable multichannel audio source separation[END_REF][START_REF] Fontaine | Multichannel Audio Modeling with Elliptically Stable Tensor Decomposition[END_REF].

Fig. 1 :

 1 Fig. 1: Generative and recognition networks. Beside each layer is indicated its size and the activation function.

+

  and H b ∈ R K×N + . It is also unsupervised in the sense that both W b

Fig. 3 :

 3 Fig.3: Results obtained with the reference and the proposed methods (using α = 1.8). The median is indicated within each boxplot.

Actually α = 1.999 because for α =

the positive α-stable distribution in (6) is degenerate. 2 https://team.inria.fr/perception/icassp2019-asvae/


Yves Le

Yves Le Jan

On Markovian random networks

We investigate the relation to random configurations and combinatorial maps of the Eulerian networks defined by Poissonian emsembles of Markov loops

Introduction

The purpose of this note is to show that random networks, which can be defined as images of Poissonian Markov loop ensembles (also known as "loop soups"), are naturally associated to random configurations and combinatorial maps. The relation with configurations, defined as families of entering and exiting half edges attached to each vertex and coupled to form edges, follows easily from the distribution of the edge occupation field, given in [START_REF] Le | Markov loops, free field and Eulerian networks[END_REF]. The relation with random combinatorial maps is more novel and should be investigated further. In the last section, we collect a few properties of these random maps which can be deduced from [START_REF] Le | Markov paths, loops and fields[END_REF].

We first present briefly the framework of our study, described in [START_REF] Le | Markov paths, loops and fields[END_REF], and recall a few useful results. Consider a graph G, i.e. a set of vertices X together with a set of non oriented edges E. We assume it is connected, and that there is no loop-edges nor multiple edges. The set of oriented edges, denoted E o , is viewed as a subset of X 2 . An oriented edge px, yq is defined by the choice of an ordering in an edge tx, yu.

Given a graph G " pX, Eq, a set of non negative conductances C x,y " C y,x indexed by the set of edges E (C x,y vanishes if tx, yu is not an edge), and a 1 non negative killing measure κ on the set of vertices X, we can associate to them an energy (or Dirichlet form) E, we will assume to be positive definite, which is a transience assumption. For any function f on X, we have:

Epf, f q " 1 2 ÿ x,y C x,y pf pxq ´f pyqq 2 `ÿ x κ x f pxq 2 .
There is a duality measure λ defined by λ x " ř y C x,y `κx . Let G x,y denote the symmetric Green's function associated with E. Its inverse equals M λ ´C with M λ denoting the diagonal matrix defined by λ. The associated λ-symmetric Markov chain is defined by the transition matrix P x y " Cx,y λy . We denote by µ the discrete loop measure associated with this Markov chain. µ is first defined on based loops by µprx 0 , x 1 , ...x n´1 , x 0 sq "

1 n P x 0 x 1 P x 1 x 2 ...P x n´1 x 0 .
A loop is an equivalence class of based loops under time shift. Given any oriented edge px, yq of the graph, we denote by N x,y plq the total number of jumps made from x to y by the loop l and set N x plq " ř y N x,y plq.

We have µptl, N x plq " nq " p1´pxq n n and therefore µptl, N x ą 0uq " ´lnpp x q, with p x " 1 λxG x,x . We refer to chapter 2 of [START_REF] Le | Markov paths, loops and fields[END_REF] for the proof. Note however that in [START_REF] Le | Markov paths, loops and fields[END_REF] we considered mostly loops in continuous time, including one-point loops. From proposition 7 and formula 2-11 in [START_REF] Le | Markov paths, loops and fields[END_REF], we get that µp1 ´e´t lx q " logp1 `tG x,x q ´logp1 `t λx q " log ´λx`tλxG x,x λx`t ¯, lx being the sum of N x plq independent exponential variables of mean λ x . We check that this equals ř 8

1 p n x n
´λx λx`t ¯n . The Poissonian loop ensemble L is defined to be the Poisson process of loops of intensity µ. We refer to [START_REF] Le | Markov paths, loops and fields[END_REF] for more details. Notations here are slightly different. Given any oriented edge px, yq of the graph, we denote by N x,y pLq or simply N x,y the total number of jumps made from x to y by the loops of L and set again N x " ř y N x,y .. We define a network to be a N-valued function defined on oriented edges of the graph. It is given by a matrix k with N-valued coefficients which vanishes on the diagonal and on entries px, yq such that tx, yu is not an edge of the graph. We say that k is Eulerian if

ÿ y k x,y " ÿ y k y,x .
For any Eulerian network k, we define k x to be ř y k x,y " ř y k y,x . It is obvious that the field N defines a random network which verifies the Eulerian property.

For a finite graph, the distribution of the random network defined by L was given in [START_REF] Le | Markov loops, free field and Eulerian networks[END_REF]. See also [START_REF] Werner | On the spatial Markov property of soups of unoriented and oriented loops[END_REF] for related Markov properties.

Theorem 1.1 For any Eulerian network k, P pN " kq " detpI ´P q ś

x k x ! ś

x,y k x,y ! ź

x,y P kx,y x,y .

Networks and configurations

The distribution of N on Eulerian networks can be obtained differently when G is finite. We introduce a specific type of configuration model. Let C be the space of configurations c defined as follows: attach to each vertex x c x entering half edges and c x exiting half edges, numbered from 1 to c x . Then provide a coupling between entering and exiting half-edges in such a way that they define an oriented edge of G.

For each oriented edge px, yq, there are c x,y half edges exiting from x coupled with c x,y half edges entering in y. For each vertex x, ř y c x,y " ř y c y,x is equal to c x and we see that such a configuration c defines a Eulerian network c. Moreover each Eulerian network k is the image by this projection map c ÝÑ c of ś

x pkx!q 2 ś px,yqPE o kx,y! different configurations. Indeed, there are

ś x kx! ś y kx,y!
way of partitioning exiting half edges and ś

x kx! ś y ky,x! ways of partitioning entering half edges according to the oriented edge of G they define, and then ś px,yqPE o k x,y ! ways to couple them.

We can now conclude easily that:

Theorem 2.1 The probability Q defined on C by Qpcq " detpI ´P q 1 ś x c x ! ź x,y P cx,y
x,y is such that, for any Eulerian network k, Qpc " kq " P pN " kq.

We say two configurations are opposite if they are exchanged by reversing the orientation of all half edges. We say they are equivalent if they can be exchanged by a circular permutation at each vertex, acting simultaneously on entering and exiting half-edges. We can deduce from this result an expression for counting configurations with given vertex degrees:

Corollary 2.1 Denote C pix,xPXq the set of configurations with vertex degrees i x : C pix,xPXq " tc P C, @x P X, c x " i x u Then the joint exponential generating function of the cardinalities of these sets,

ř pix,xPXq |C pix,xPXq | ś x s ix x ix! , defined for 0 ď s x ă 1
dx , is given by the inverse of the determinant of the matrix δ x,y ´sx A x,y , A being the adjacency matrix of G.

Proof. Taking unit unit conductances and killing rates 1 ´dx s x , we rewrite ř c Qpcq " 1 in terms of cardinalities.

Remark: A similar construction can be made with the even networks defined by the edge occupation field N p 1 2 q of a loop ensemble L p 1 2 q of intensity

1 2 µ (Cf [5]). N p 1 2 q
tx,yu is the total number of jumps from x to y or from y to x made by the loops of L p 1 2 q . The space of configurations is now the set C ev of numbered G-maps with an even number of incident edges at each vertex. At each vertex x there are 2kx! ś y k tx,yu ! way of partitioning the 2k x incident edges and there are ś tx,yuPE k tx,yu ! ways to couple the numbers assigned to the edges between x and y. Each even network k is the image of

ś x p2kx!q ś tx,yuPE k tx,yu !
different configurations and the distribution of these configurations is uniform conditionally to the network.

Configurations and Wilson algorithm

Let us recall the relation given in section 8-2 of [START_REF] Le | Markov paths, loops and fields[END_REF] between loop ensembles and Wilson algorithm (cf [START_REF] Bruce | Generating random spanning trees more quickly than the cover time[END_REF]), for a finite graph. The choice of an arbitrary order on the vertices x 1 , x 2 ... allows to construct from a finite sequence of random walks a random spanning tree together with an independent set of loops l i , based at x i and included in the complement of tx 1 , ... , x i´1 u. Note that l i can be empty as we are here in discrete time. The random network defined by this set of loops has the same distribution as N. Moreover, all such sets of loops defining the same network have equal probability. In order to construct not only the network N but the Poissonian ensemble L, these based loops l i have to be randomly divided at their base points x i . This partition is defined by a Poisson-Dirichlet distribution if we consider loops indexed by continuous time, with exponential holding time: see remark 21 in section 8-2 of [START_REF] Le | Markov paths, loops and fields[END_REF]. As mentioned in [START_REF] Le | Markov loops topology[END_REF] there is a discrete version of this splitting method: If l i visits x i n x i times, we can view each based loop l i as a set of n x i excursions out of x i and partition it randomly according to the exchangeable partition probability function (EPPF)

ś k 1 pn j ´1q! nx i !
occurring in the so-called "chinese restaurant process" ([8], sections 2-1 and 3-1). The EPPF

ś k 1 pn j ´1q! nx i !
is the probability of a partition into a set of size n 1 , a set of size n 2 etc . The distribution of the block sizes, presented in exchangeable random order is

1 k! ś k 1 n j (cf section 2-1 in [8]
). The proof is based on the fact that for any partition tm 1 , m 2 , ...m k u of m, the probability that the k-tuple m 1 , ... m k is equal to the set tN x plq, l P Lu in random order equals P p|tl P L, N x plq ą 0qu| " kq

ś k j"1 µpNxplq"m j q µpNxplqą0q " px k! ś k j"1 p1´pxq m j m j
. As P pN x pLq " mq " p x p1´p x q m , conditionally to N x pLq " m, this probability equals

1 k! ś k 1 m j
. By exchangeability, performing the partitition in tm 1 , m 2 , ...m k u based loops according to the original order in the excursion set yields the same distribution. We then concatenate the excursions in each part according to their original order and define L to be the associated set of (unbased) discrete time loops. Its distribution is uniform conditionally on the network.

Note that his construction can be related to random configurations. Let us call exit configuration a class of configurations with the same partition of exiting half edges. More precisely, an exit configuration is determined by assigning to each exiting half edge the vertex to which the oriented edge will lead, without specifying the entering half edge of that vertex. All configurations c in such a class c ex have equal probability and the image of Q on the set of exit configurations C ex is Qpc ex q " detpI ´P q ś x,y P cx,y

x,y . All exit configurations inducing the same network have equal probability. Such an exit configuration c ex , with the choice of an order on X, allows to construct a sample of the family of based loops l i considered above. The first loop is based at the first vertex x 1 and visits it c ex x 1 times. It starts with the first exiting half-edge, continues to the vertex to which it is associated, then with the first half edge exiting this vertex and so on, using and taking at each vertex exiting half edges in increasing order until all have been used at x 1 . Then iterate this procedure with the remaining configuration, starting at the first vertex it contains, according to the order initially defined on X. This correspondence is a bijection, and if we start with a random exit configuration, the distribution obtained on the sequence of based loops is the same as in Wilson's algorithm. Note for the following that we get the same distribution if exiting half edges are used in random order. Partitioning these based loops as above yields a set of unbased loops L ex distributed as L. An entrance configuration c in is defined in an analogous way, in terms of entering half edges, assigning to each of them the vertex from which the oriented edge starts. Similarly, we can define L in and we get in that way two sets of loops which are clearly independent conditionally on the network c.

Networks and combinatorial maps

The relations between configurations and Poissonian loop ensembles appears to be deeper, as shown in the following.

-Recall that a combinatorial map can be defined as a graph, with possibly multiple edges between vertices, equipped with a combinatorial embedding, i.e. a cyclic ordering of edges around each vertex (see [START_REF] Mohar | Graphs on Surfaces[END_REF], [START_REF] Lando | Graphs on Surfaces and their Applications[END_REF]). This ordering allows to define combinatorial faces (termed as facial walks in section 4-1 of [START_REF] Mohar | Graphs on Surfaces[END_REF]) as cycles of the permutation φ on oriented edges obtained by composing the edge orientation reversal ρ by the shift defined by the cyclic order. Note that the facial walks may visit a vertex several times. The map can be drawn on a surface with possibly several connected components (see [START_REF] Mohar | Graphs on Surfaces[END_REF], p.85). Let us say it is a G-map if its vertex set is contained in X and if its edges are multiple copies of elements of E (if G is complete, the second condition is redundant). Let us say a map is numbered if a first edge is chosen at every vertex.

-An element c of C defines a numbered G-map Mpcq with an even number of incident edges at each vertex. The first edge at each vertex is determined by the first exiting half edge and the half edge it is coupled to, the second by the first entering half edge and the half-edge it is coupled to, the third by the second exiting half edge and the half-edge it is coupled to, and so on, alternating entering and exiting half edges according to their original cyclic order. We actually get in this way, an orientation on the edges of the map. At each vertex, the cyclic order on edges alternates exiting and entering edges. The number of oriented edges carried by px, yq is equal to cx,y . Two equivalent configurations define the same G-map.

-An equivalent class of configurations defines also alternating signs on the faces of its associated map. The positive faces F `are defined by the cycles of the permutation on edges obtained by mapping each exiting half-edge with the exiting half-edge following, (in the above defined cyclic order at each vertex), the entering half-edge coupled to it. Oriented face contours are defined by the orientation of their as explained at the end of the previous section (the exiting half edges here are used in random order). Conditionally on these based loops, the random ordered partitions of them we have obtained at each vertex are exactly the size-biased presentation of the random partition defined by the "chinese restaurant process' (see section 2-1 in [START_REF] Pitman | Combinatorial Stochastic Processes. 32th St Flour Summer School[END_REF]). This observation completes the proof of the theorem as the case of L ´can be treated in the same way.

Expectations calculations

The study of the geometrical properties of the random combinatorial maps defined by random configurations should be pushed further. In this last section, we present a few results in this direction. Using the results of [START_REF] Le | Markov paths, loops and fields[END_REF], expressions for the expected number of vertices, edges and faces of Mpcq can be given. These are respectively: Ep|tx P X, c x ą 0u|q " ř x p1 ´1 λxG x,x q (as N x follows a geometric distribution of parameter 1 λxG x,x (see the remark following formula 4-3 in [START_REF] Le | Markov paths, loops and fields[END_REF]), Ep ř x,y cx,y q " ř x,y C x,y G x,y (see formula 2-10 in [START_REF] Le | Markov paths, loops and fields[END_REF]) and Ep|L `pcq|`|L ´pcqq " ´2 ln detpI ´P q (see formula 2-5 in [START_REF] Le | Markov paths, loops and fields[END_REF]) . Similar expressions can be given for the variances of the vertices and edges numbers. Much more intricate expressions can also be given for the number and the sizes of the connected components of the map (see [START_REF] Le | Markovian loop clusters on graphs[END_REF]).

All these results apply in particular to the case of complete graphs, which is of special interest in relation with the theory of combinatorial maps. For d vertices and constant killing rate κ, λ " d ´1 `κ and G " 1 d`κ pI `1 κ Jq, J denoting the pd, dq matrix with all entries equal to 1. Setting p " 1 d´1`κ , detpI ´P q " p1 ´pd ´1qpqp1 `pq d´1 . Inserting this in the expressions recalled hereinabove, it appears that if we let d and κ increase to infinity, with κ d converging to zero, the expected number of vertices and the expected number of edges are both equivalent to d κ . More precisely, they are respectively equal to d pκ`1q `1 `op1q and d κ ´1 `op1q. The expected number of faces is 2 lnpd{κq ´2 `op1q.

Following section 4-1 in [START_REF] Mohar | Graphs on Surfaces[END_REF], the Euler characteristic χpcq of the combinatorial map Mpcq is defined as the number of vertices + the number of faces -the number of edges, i.e |tx P X, c x ą 0u| `|L `pcq| `|L ´pcq| ´Npcq with Npcq " ř x,y cx,y .

Its expectation can be expressed as

Epχq " ÿ x p1 ´1 λ x G x,x q ´2 ln detpI ´P q ´ÿ x,y C x,y G x,y .

In particular, for the complete graph, Epχq " dp1 ´p1 ´1 κ`1 qp1 `1 d´1`κ qq dpd´1q

κpd`κq ´2 lnp κ d´1`κ q ´2pd ´1q lnp1 `1 d´1`κ q. The expectation of the characteristic χ can converge to any value v by taking κpd, vq " c d upd, vq with upd, vq ą 1 defined as the unique solution of the equation upd, vq ĺnpupd, vqq " lnpdq ´v. Note that upd, vq is equivalent to lnpdq as d increases to infinity. If we take κpdq " c d lnpdq , the expectation of χ is equivalent to lnplnpdqq as d increases to infinity. Note however that the number of faces in the 2-cell embedding in surfaces defined by the combinatorial map is larger than the number of facial walks. Hence the sum of the Euler characteristics of the surfaces is smaller or equal to χpcq (see proof of 3-4-1 in [START_REF] Mohar | Graphs on Surfaces[END_REF]). The average number of vertices visited by a facial walk is

ř x µpNxą0q µp1q
" ´řx lnpP pNx"0qq ´ln detpI´P q " ř x lnpλxG x,x q ´ln detpI´P q It is smaller from the average length of the walk which is ř

x µpNxq µp1q " ř x pλxG x,x ´1q ´ln detpI´P q . The difference is the average number of multiple points on the loop, counted with their multiplicities minus 1. For the complete graph, with d and κ increasing to infinity, with κ d converging to zero, we get that the is the average number of multiple points is equivalent to d 2κ 2 lnp d κ q which converges to 1 2 if we take κpdq " c d lnpdq .

Consider finally essential vertices, i.e. such that N x ą 1. Their expected number in general is ř x p1 ´1 λxG x,x q 2 which equals dpd`2κ´1q 2 pκ`1q 2 pd´1`κq 2 in the complete graph case. If we take κpdq "

  faces and the expected number of essential vertices are both equivalent to lnpdq as d increases to infinity.

Key words and phrases: Markov loops, networks, maps 0 AMS 2010 subject classification: 60C05, 60J27, 60G60.

Acknowledgement: Thanks are due to the referees for suggesting several improvements of this paper.

edges. The negative faces F ´are defined by the cycles of the permutation on edges obtained by mapping an entering half-edge with the entering half-edge following the exiting half-edge coupled to it. Oriented face contours are defined by reversing the orientation of their edges. Taking the opposite configuration exchanges F `and F ´. Each edge of the map is in this way adjacent to two faces of opposite signs, and giving it opposite orientations as part of a face contour.. Note that F determines the configuration, up to equivalence. The same holds for F ´. The image on G of the sets of oriented face contours of F `and F ´are multisets of loops denoted pL `and L ´q.

-We can now state our main result: Theorem 4.1 The loop ensembles L `and L ´defined by a random configuration under the probability Q have both the same distribution as L .

Proof. Note first that (in contrast with the original version of the algorithm presented in the previous section), this construction does not require additional independent partitions. We can again choose an arbitrary order on the vertices to run a different version of Wilson algorithm to construct L `. This algorithm can be viewed as a way, starting from c ex (i.e. the partition of exiting half edges) to sample the random coupling with entering edges progressively. This determines a configuration c sampled uniformly in c ex . We start a based loop at the first vertex with the first exiting half edge, then move to the vertex attached to it by the partition, and choose an entering half edge to couple it. Then follow the exiting half edge following this entering edge in the alternating cyclic order, and so on until we reach again the base point. Then, a based loop is created if the first exiting half-edge follows the last entering one, in the cyclic order. We iterate the process until this happens, and create a based loop by concatenating the excursions we constructed. Note that conditioning on the initially given partition of exiting half edges and on the partial coupling with entering half edges done until a return time (this includes the excursions out of x 1 up to this return time), the creation of a based loop occurs with probability one over the number of unused entering half edges at the base point. Then restart at the same vertex with the next free exiting half-edge. Iterate until all half edges at the first vertex have been used.

-We can then restart the algorithm at the next vertex with unused exiting half-edge, starting with the first one. Iterate until all half edges of all vertices have been used. If we start with a random initial exit configuration, the concatenation of the excursions obtained at each vertex has the same distribution as the family of based loops tl i u defined by Wilson algorithm