Xiaofei Li

Yutong Ban

Laurent Girin

Xavier Alameda-Pineda

Radu Horaud

A CASCADED MULTIPLE-SPEAKER LOCALIZATION AND TRACKING SYSTEM

Keywords: sound-source localization, multiple moving speakers, tracking, reverberant environments, LOCATA Challenge

This paper presents an online multiple-speaker localization and tracking method, as the INRIA-Perception contribution to the LO-CATA Challenge 2018. First, the recursive least-square method is used to adaptively estimate the direct-path relative transfer function as an interchannel localization feature. The feature is assumed to associate with a single speaker at each time-frequency bin. Second, a complex Gaussian mixture model (CGMM) is used as a generative model of the features. The weight of each CGMM component represents the probability that this component corresponds to an active speaker, and is adaptively estimated with an online optimization algorithm. Finally, taking the CGMM component weights as observations, a Bayesian multiple-speaker tracking method based on the variational expectation maximization algorithm is used. The tracker accounts for the variation of active speakers and the localization miss measurements, by introducing speaker birth and sleeping processes. The experiments carried out on the development dataset of the challenge are reported.

INTRODUCTION

For multiple-speaker localization, the W-disjoint orthogonality (WDO) [START_REF] Yilmaz | Blind separation of speech mixtures via timefrequency masking[END_REF] assumption is widely used. It assumes that the audio signal is dominated by only one speaker in each small region of the time-frequency (TF) domain, because of the natural sparsity of speech signals in this domain. After applying the short-time Fourier transform (STFT), interchannel localization features, e.g. interaural phase difference, can be extracted. To assign the interchannel features to multiple speakers, a mixture of Gaussian mixture models (GMMs) is used as a generative model of the interchannel features of multiple speakers in [START_REF] Mandel | Model-based expectationmaximization source separation and localization[END_REF] with each GMM representing one speaker, and then feature assignment was done based on the maximum likelihood criteria. In [START_REF] Dorfan | Tree-based recursive expectation-maximization algorithm for localization of acoustic sources[END_REF], instead of setting one GMM for each speaker, a single complex GMM (CGMM) is used for all speakers with each component representing one candidate speaker location. After maximizing the likelihood of the features, with an expectation-maximization (EM) algorithm, the weight of each component represents the probability that there is an active speaker at the corresponding candidate location. To localize moving speakers, in [START_REF] Schwartz | Speaker tracking using recursive EM algorithms[END_REF], based on a CGMM model similar to [START_REF] Dorfan | Tree-based recursive expectation-maximization algorithm for localization of acoustic sources[END_REF], a recursive EM algorithm was proposed to update online the CGMM component weights. Counting and localization of active speakers can be jointly carried out by selecting the components with large weights. Taking This work was supported by the ERC Advanced Grant VHIA #340113.

the instantaneous outputs of a localization method as observations and using a speaker dynamic model, Bayesian tracking techniques estimate the posterior distribution of source locations, e.g. [START_REF] Vermaak | Nonlinear filtering for speaker tracking in noisy and reverberant environments[END_REF][START_REF] Valin | Robust localization and tracking of simultaneous moving sound sources using beamforming and particle filtering[END_REF]. To tackle the tracking problem with an unknown and time-varying number of speakers, additional model features such as observationto-speaker assignments, speaker track birth/death processes and a model of speech activity can be included [START_REF] Evers | Bearing-only acoustic tracking of moving speakers for robot audition[END_REF][START_REF] Ba | An on-line variational bayesian model for multi-person tracking from cluttered scenes[END_REF][START_REF] Gebru | Audio-visual speaker diarization based on spatiotemporal Bayesian fusion[END_REF][START_REF] Ban | Exploiting the complementarity of audio and visual data in multi-speaker tracking[END_REF].

In the present paper, we present our online multiple-speaker localization and tracking method contributed to the LOCATA Challenge, which is composed of three modules:

• A recursive DP-RTF (direct-path relative transfer function) estimation module. In the STFT domain, the room impulse response (RIR) can be approximated by the convolutive transfer function (CTF) [START_REF] Avargel | System identification in the short-time Fourier transform domain with crossband filtering[END_REF][START_REF] Talmon | Relative transfer function identification using convolutive transfer function approximation[END_REF]. DP-RTF is defined by the ratio between the first taps of the CTF of two microphones, thus encodes the direct-path information and is used as an interchannel feature being robust against reverberation. The CTF estimation used for DP-RTF extraction was formulated in batch mode in [START_REF] Li | Estimation of the direct-path relative transfer function for supervised sound-source localization[END_REF][START_REF]Multiple-speaker localization based on direct-path features and likelihood maximization with spatial sparsity regularization[END_REF] with speakers considered as static. Based on recursive least-square (RLS), the online CTF estimation method was proposed in [START_REF] Li | Online localization of multiple moving speakers in reverberant environments[END_REF][START_REF] Li | Online Localization and Tracking of Multiple Moving Speakers in Reverberant Environment[END_REF], and will be briefly presented in this paper, for moving speaker localization.

• An online multiple-speaker localization module. In [START_REF]Multiple-speaker localization based on direct-path features and likelihood maximization with spatial sparsity regularization[END_REF], we adopt the above-mentioned CGMM model [START_REF] Dorfan | Tree-based recursive expectation-maximization algorithm for localization of acoustic sources[END_REF] to assign the DP-RTF features to speakers, in addition, an entropy-based regularization term was used to impose the spatial sparsity of the estimated component weights. However, [START_REF]Multiple-speaker localization based on direct-path features and likelihood maximization with spatial sparsity regularization[END_REF] only considered the batch mode and static speakers. The recursive EM algorithm proposed in [START_REF] Schwartz | Speaker tracking using recursive EM algorithms[END_REF] was adopted in [START_REF] Li | Online localization of multiple moving speakers in reverberant environments[END_REF] for online liklihood maximization without using the entropy regularization. Furtherly, an online optimization method, i.e. exponentiated gradient (EG) [START_REF] Kivinen | Exponentiated gradient versus gradient descent for linear predictors[END_REF], was used in [START_REF] Li | Online Localization and Tracking of Multiple Moving Speakers in Reverberant Environment[END_REF] for simultaneously online liklihood maximization and entropy minimization. This method will be briefly presented in this paper.

• A multiple-speaker tracking module. In [START_REF] Li | Online Localization and Tracking of Multiple Moving Speakers in Reverberant Environment[END_REF], a multiple-speaker tracking method was proposed. The results of the above localization module are taken as inputs and efficiently exploited in a Bayesian framework. The problem is efficiently solved by a variational expectation maximization (VEM) algorithm. Speaker birth and sleeping processes are included in the tracking process. The sleeping process is efficient to tackle a missed detection by the localization procedure. This paper will briefly present this multiple-speaker tracker, please refer to [START_REF] Ba | An on-line variational bayesian model for multi-person tracking from cluttered scenes[END_REF][START_REF] Ban | Exploiting the complementarity of audio and visual data in multi-speaker tracking[END_REF][START_REF] Ban | Tracking a varying number of people with a visually-controlled robotic head[END_REF][START_REF] Li | Online Localization and Tracking of Multiple Moving Speakers in Reverberant Environment[END_REF] for more detailed algorithmic derivation and description.

In Fig. 1, the three layers of temporal evolutions depicts these three modules, respectively. In the following, we present them one by one from Section 2 to 4, and then gives the experiments on the development dataset of the LOCATA challenge [START_REF] Löllmann | The LOCATA challenge data corpus for acoustic source localization and tracking[END_REF] in Section 5.

RECURSIVE DP-RTF ESTIMATION

To simplify the presentation, let us first consider the noise-free single-speaker case. In the time domain, the i-th, i = 1, . . . , I, microphone signal is:

x i (n) = a i (n) s(n),
where n is the time index, s(n) is the source signal, a i (n) is the RIR from the source to the i-th microphone, and denotes the convolution. Applying the STFT, and using the CTF approximation, we have for each frequency index f = 0, . . . , F -1:

x i t,f = a i t,f s t,f
, where x i t,f and s t,f are the STFT coefficients of the corresponding signals, and the CTF a i t,f is a subband representation of a i (n). Here, the convolution is executed w.r.t the frame index t. The first CTF coefficient a i 0,f mainly consists of the direct-path information, and DP-RTF is defined as the ratio between the first CTF coefficients of two channels: a i 0,f /a r 0,f , where channel r is the reference channel. Based on the cross-relation method [START_REF] Xu | A least-squares approach to blind channel identification[END_REF], for one microphone pair (i, j), we have x i t,f a j f = x j t,f a i f , where denotes matrix/vector transpose, the convolution vectors are a i f = [a i 0,f , . . . , a i Q-1,f ] and x i t,f = [x i t,f , . . . , x i t-Q+1,f ] , with Q denoting the CTF length. We concatenate the CTF vector of all channels as a f = [a 1 f , . . . , a I f ] . For each microphone pair (i, j), we construct a cross-relation equation in terms of a f . For this aim, define:

x ij t,f = [0, . . . , 0 (i-1)Q , x j t,f , 0, . . . , 0 (j-i-1)Q , -x i t,f , 0, . . . , 0 (I-j)Q ] . (1)
Then we have x ij t,f a f = 0. The CTF vector a f can be estimated by solving this equation. To avoid a trivial solution, i.e. a f = 0, we constrain the first CTF coefficient of the reference channel, say r = 1, to be equal to 1. This leads to the following equation

xij t,f ãf = y ij t,f , (2)
where -y ij t,f is the first entry of

x ij t,f , xij t,f is x ij t,f
with the first entry removed, and ãf is the relative CTF vector:

ãf = ã1 f a 1 0,f , a 2 f a 1 0,f , . . . , a I f a 1 0,f . (3)
where

ã1 f = [a 1 1,f , . . . , a 1 Q-1,f ]
, is a 1 f with the first entry removed. For i = 2, . . . , I, the DP-RTFs appear in (3) as the first entries of

a i f a 1 0,f .
The DP-RTF estimation amounts to solving the linear problem Eq. ( 2). Eq. ( 2) is defined for one microphone pair at one frame. For the online case, we would like to update the estimate of ãf using the current frame, say t. There is a total of M = I(I -1)/2 distinct microphone pairs, For notational convenience, instead of using ij , we use m = 1, . . . , M denote the index of microphone pair. The fitting error of ( 2) is e m t,f = y m t,f -xm t,f ãf . At the current frame t, for the microphone pair m, RLS aims to minimize the error

J m t,f = t-1 t =1 M m =1 λ t-t |e m t ,f | 2 + m m =1 |e m t,f | 2 , ( 4
)
which sums up the fitting error of all the microphone pairs for the past frames and the microphone pairs up to m for the current frame.

The forgetting factor λ ∈ (0, 1] gives exponentially lower weight to older frames, whereas at one given frame, all microphone pairs have the same weight. In RLS, this minimization problem can be recursively solved. For the detailed recursion procedure please refer to Algorithm 1 in [START_REF] Li | Online Localization and Tracking of Multiple Moving Speakers in Reverberant Environment[END_REF]. For each frame t, let ãt,f denote the CTF estimate, and ci t,f , i = 2, . . . , I denote the DP-RTF estimates extracted from ãt,f . Note that implicitly we have c1 t,f = 1.

We now introduce how to extend the above method to the noisy multiple-speaker case. We assume that, over a short time, the speakers are static and only one source is active at each frequency bin. Therefore, the CTF can be estimated using a small number of recent frames. This can be done by adjusting the forgetting factor λ: To approximately have a memory of P frames, we can set λ = P -1 P +1 . P should be empirically set to achieve a good tradeoff between the validity of the above assumptions and a robust CTF estimate. To suppress the noise, we use the inter-frame spectral subtraction algorithm proposed in [START_REF] Li | Estimation of relative transfer function in the presence of stationary noise based on segmental power spectral density matrix subtraction[END_REF][START_REF] Li | Online localization of multiple moving speakers in reverberant environments[END_REF]. At each frequency bin f , frames are first classified into speech frames and noise frames, then interframe spectral subtraction is applied between them. In the RLS process, only the speech frames (after spectral subtraction) are used, and the noise frames are skipped. In practice, a DP-RTF estimate can sometimes be unreliable due to the imperfect performance of the above-mentioned methods. We use the consistency test method proposed in [START_REF]Multiple-speaker localization based on direct-path features and likelihood maximization with spatial sparsity regularization[END_REF][START_REF] Li | Online Localization and Tracking of Multiple Moving Speakers in Reverberant Environment[END_REF] to detect the unreliable estimates. Finally, at frame t, we obtain a set of features Ct = {{ĉ i t,f }i∈I f } F -1 f =0 , where I f ⊆ {2, . . . , I} denotes the set of microphone indices that pass the consistency test. Each of the features is assumed to be associated with a single speaker.

ONLINE MULTIPLE-SPEAKER LOCALIZATION

In order to assign the DP-RTF features in Ct to speakers, the CGMM generative model proposed in [START_REF] Dorfan | Tree-based recursive expectation-maximization algorithm for localization of acoustic sources[END_REF] is adopted. We define a set D of D candidate source locations. Let d = 1, . . . , D denote the location index. The probability, that an observed feature ĉi t,f is emitted by candidate locations, is modelled by a CGMM:

P (ĉ i t,f |D) = D d=1 w d Nc(ĉ i t,f ; c i,d f , σ 2 ), (5)
where the mean c i,d f is precomputed based on the direct-path propagation model from the d-th candidate location to the microphone i. The variance σ 2 is empirically set as a constant value. The component weight w d ≥ 0, with D d=1 w d = 1, is the only free model parameter. Let us denote the vector of weights with w = [w 1 , ..., w D ] .

Let Lt denote the log-likelihood (normalized by the number of features) of the features in Ct, as a function of w d . Once Lt is maximized, the component weight w d represents the probability that there exist an active speaker at the d-th candidate location. In addition, taking into account the fact that the number of actual active speakers is much lower than the number of candidate locations, an entropy minimization is taken as a regularization term to impose a sparse distribution of w d . The entropy is defined as H = -D d=1 w d log(w d ). Then the overall cost function is -Lt + γH. In the present online framework, we proceed to a recursive update of the weight vector at current frame t, hence now denoted wt, based on the estimate at previous frame wt-1 and on the DP-RTF features of the current frame. This can be formulated as an online optimization problem [START_REF] Kivinen | Exponentiated gradient versus gradient descent for linear predictors[END_REF]:

wt = argmin w χ(w, wt-1) + η(-Lt + γH), (6)
with the constraints that w d , d = 1, . . . , D are positive and sum to 1. χ(w, wt-1) is some distance measure between w and wt-1.

The positive constant η controls the parameter update rate. To exploit the fact that w d are probability masses, we use the Kullback-Leibler divergence, i.e. χ(w, wt-1)

= D d=1 w d log w d w d t-1
, which results in the exponentiated gradient algorithm [START_REF] Kivinen | Exponentiated gradient versus gradient descent for linear predictors[END_REF]. Let ∆t-1 =

∂(η(-L t +γH)) ∂w d w d t-1
denote the partial derivatives of η(-Lt + γH) w.r.t w d at the point w d t-1 . Then, the exponentiated gradient r d t-1 = e -∆ t-1 is used to update the weights

w d t = r d t-1 w d t-1 D d =1 r d t-1 w d t-1 , ∀d ∈ D. (7)
It is obvious from [START_REF] Evers | Bearing-only acoustic tracking of moving speakers for robot audition[END_REF] that the parameter constraints, namely positivity and summation to 1, are automatically satisfied.

MULTIPLE-SPEAKER TRACKING

In the following, upper case letters denote random variables while lower case letters denote their realizations. Let N be the maximum number of tracks (speakers), and let n be the speaker index. Moreover, let n = 0 denote no speaker, or background noise. Let Stn be a latent (or state) variable associated with speaker n at frame t, and let St = (St1, . . . , Stn, . . . , StN ). Stn is composed of two parts: the speaker direction and the speaker velocity. In this work, speaker direction is defined by an azimuth θtn. To avoid phase (circular) ambiguity we describe the direction with the unit vector U tn = (cos(θtn), sin(θtn)) . Moreover, let Vtn ∈ R be the angular velocity. Altogether we define a realization of the state variable as stn = [utn; vtn] where the notation [•; •] stands for vertical vector concatenation.

As mentioned above, the CGMM component weight w d represents the probability that there exist an active speaker at the d-th candidate location. The frame-wise localization of active speakers can be carried out by peak picking over w d . However, to fully use the weight information, without applying peak picking, all the candidate locations and their associated weights are used. Formally, let Ot = (Ot1, . . . , O td , . . . , OtD) be the observed variables at frame t. Each realization o td of O td is composed of a candidate location, or azimuth θtd ∈ D, and a weight w td . As above, let the azimuth be described by a unit vector b td = (cos( θtd ), sin( θtd )) . In summary we have o td = [b td ; w td ]. Moreover, let Z td be a (latent) assignment variable associated with each observed variable O td , such that Z td = n means that the observation indexed by d at frame t is assigned to speaker n ∈ {0, . . . , N }.

Bayesian Tracking Model

The problem at hand can be cast into the estimation of the filtering distribution p(st, zt|o1:t), and further inference of st and zt. By applying the Bayes rule, the filtering distribution is proportional to: p(st, zt|o1:t) ∝ p(ot|st, zt)p(zt)p(st|o1:t-1), [START_REF] Ba | An on-line variational bayesian model for multi-person tracking from cluttered scenes[END_REF] which contains the following three terms.

The audio observation model p(ot|st, zt) describes the distribution of the observations given speakers state and assignment. We assume the different observations are independent conditionally to speakers state and assignment. For each observation, we adopt the weighted-data GMM model of [START_REF] Gebru | Em algorithms for weighted-data clustering with application to audio-visual scene analysis[END_REF]: p(b td |Z td = n, stn; w td ) = N (b td ; Mstn, 1 w td Σ) for n ∈ {1, . . . , N }, where the matrix M = [I2×2, 02×1] projects the state variable onto the space of source directions and Σ is a covariance matrix (set empirically to a fixed value). Note that the weight plays the role of a precision: The higher the weight w td , the more reliable the source direction b td . The case Z td = 0 follows a uniform distribution over the volume of the observation space, i.e. p(b td |Z td = 0) = U(vol(G)).

The prior distribution of the assignment variable, i.e. p(zt), is independent over observations and is assumed to be uniformly distributed over all the speakers, i.e. p(Z td = n) = π dn = 1 N +1 .

To calculate the state predictive distribution p(st|o1:t-1), we marginalize it over st-1: p(st|o1:t-1) = p(st|st-1)p(st-1|o1:t-1)dst-1.

(

) 9
We model the state dynamics p(st|st-1) as a linear-Gaussian first-order Markov process, independent over the speakers, i.e. p(st,n|st-1,n) = N (stn; Dt-1,nst-1,n, Λtn), where Λtn is the dynamics' covariance matrix and Dt-1,n is the state transition matrix. Importantly, since the state position subvector utn lies on the unit circle, the dynamic model is designed for circular motion.

Given the estimated azimuth angle θt-1,n at frame t -1, the state transition matrix can be written as:

Dt-1,n =   1 0 -sin(θt-1,n) 0 1 cos(θt-1,n) 0 0 1   . ( 10
)
In the following Dt-1,n is written as D for notational simplicity.

Variational Expectation Maximization Algorithm

To estimate the model parameters Θ, and infer the posteriori distribution of (st, zt), we adpot an EM algorithm that alternates between computing and maximizing the expected complete-data loglikelihood:

J(Θ, Θ o ) = E p(z t ,s t |o 1:t ,Θ o ) [log p(zt, st, o1:t|Θ)] , (11)
where E denotes expectation, Θ o are the "old" model parameter estimates (obtained at previous iteration). Given the hybrid combinatorial-continuous nature of the latent space, it is computationally heavy to operate with the exact a posterior distribution p(zt, st|o1:t). We thus use a variational approximation to solve the problem [START_REF] Ba | An on-line variational bayesian model for multi-person tracking from cluttered scenes[END_REF], which factorizes the posterior distribution as:

p(zt, st|o1:t) ≈ q(zt, st) = q(zt) N n=0 q(stn). (12)
It is seen that the posterior distribution factorizes across speakers. This principle is also valid for time t -1. As a result, the predictive distribution (9) also factorizes across speakers. If the posterior distribution q(st-1,n) is assumed to be a Gaussian with mean µ t-1,n and variance Γt-1,n, then the per-speaker predictive distribution p(stn|o1:t-1) is a Gaussian:

p(stn|o1:t-1) = N (stn; Dµ t-1,n , DΓt-1,nD + Λtn). (13)
In E-S step, it can be derived that the variational posterior distribution of q(stn) is a Gaussian distribution, with variance and mean respectively as

Γtn = D d=1 α tdn w d t H Σ -1 H + Λtn + DΓt-1,nD -1 -1
,

µ tn = Γtn H Σ -1 D d=1 α tdn w d t b td + Λtn + DΓt-1,nD -1 Dµ t-1,n ,
where α tdn = q(Z td = n) is the variational posterior distribution of the assignment variable, which is derived in E-Z step, the model parameters Λtn is udated in M-step, please refer to [START_REF] Li | Online Localization and Tracking of Multiple Moving Speakers in Reverberant Environment[END_REF] for detailed descriptions. At time t, VEM converges after a few iterations, e.g. 5 iterations used in this work. After convergency, the posterior mean of stn, i.e. µ t = (µ t1 , . . . , µ tn , . . . , µ tN ), is output as the result of speaker tracking, and together with the posterior variance Γt = (Γt1, . . . , Γtn, . . . , ΓtN ) are transmitted to time t + 1.

Speaker Track Birth Process and Activity Detection

A track birth process is used to initialize new tracks, i.e. new speakers that enter the scenario. The general principle is the following. In a short period of time, say from frame t -L to frame t, we assume that at most one new (yet untracked) speaker appears. For each frame from t -L to t, we select the observation assigned to the background noise with the highest weight, and obtain an observation sequence õt-L:t. We then calculate the marginal likelihood of this sequence according to our model, τ0 = p(õt-L:t). If these observations have been generated by a new speaker, they exhibit smooth trajectories, and τ0 will be high. Therefore, the birth process is conducted by comparing τ0 with a threshold. The posterior distribution of the assignment variable, i.e. α idn , can be used for multiple-speaker activity detection. This can be formalized as testing, for each frame t and each speaker n, whether the weighted assignments D d=1 α idn w d i (averaged over a small number of frames) is larger than a threshold.

EXPERIMENTS ON LOCATA DEVELOPMENT DATA

We report the results of on the LOCATA development corpus for tasks #3 and #5 with a single moving speaker, and tasks #4 and #6 with two moving speakers, each task comprising three recorded sequences. In this work, we use four microphones with indices {5, 8, 11, 12} out of the twelve microphones of a spherical array built in the head of a humanoid robot, i.e. NAO, to perform azimuth localization and tracking. These four microphones are mounted on the top of the robot head, and they approximately lie in a horizontal plane. We perform 360 • -wide azimuth estimation and tracking: D = 72 azimuth directions at every 5 • in [-175 • , 180 • ] are used as candidate directions. The TDOAs are computed based on the coordinate of microphones, which are then used to compute the phase of the CGMM means, while the magnitude of the CGMM means are set to a constant for all the frequencies. All the recorded signals are resampled to 16 kHz. The STFT uses the Hamming window with length of 16 ms and shift of 8 ms. The CTF length is Q = 8 frames. The RLS forgetting factor λ is computed using ρ = 1. The exponentiated gradient update factor is η = 0.07. The entropy regularization factor is γ = 0.1. For the tracker, the covariance matrix is set to be isotropic Σ = 0.03I2. A localization and tracking example Fig. 2 shows an example for a LOCATA sequence. Two speakers are moving and continuously speaking with short pauses. Fig. 2 (b) shows that the localization method achieves a heatmap with negligible interferences and smooth peak evolution. From Fig. 2 (c), it is seen that the tracking method further smooth the speaker moving trajectories. Even when the observations have a low weight, the tracker is still able to give the correct speaker trajectories. This is ensured by exploiting the source dynamics model. As a result, the tracker is able to preserve the identity of speakers in spite of the (short) speech pauses. In the presented sequence example, the estimated speaker identities are quite consistent with the ground truth.

Quantitative results

The following metrics are used for quantitative evaluation. The detected speaker is considered to be successfully localized if the azimuth difference is not larger than 15 • . The absolute error is calculated for the successfully localized sources. The mean absolute error (MAE) is computed by averaging the absolute error of all speakers and frames. For the unsuccessful localizations, we count the miss detections (MD) (speaker active but not detected) and false alarms (FA) (speaker detected but not active). Then the MD and FA rates are computed, using all the frames, as the percentage of the total MDs and FAs out of the total number of actual speakers, respectively. In addition to these localization metrics, we also count the identity switches (IDs) to evaluate the tracking continuity. IDs represents the number of the identity changes in the tracks for a whole test sequence. Table 1 gives the quantitative localization and tracking results. The localization results are obtained by applying peak picking on the GMM weights. The tracker slightly reduces the FA rate compared to the localization method alone mainly by eliminating some spurious peaks that are present in the localization outputs. It also reduces the MD rate since some correct speaker trajectories can be recovered even when the observations have (very) weak weights, as explained above. The identity switches are mainly due to the crossing of speaker trajectories, a hard case for the source dynamics model.

CONCLUSION

In this paper, we presented the INRIA-Perception contribution to the LOCATA Challenge 2018. We combined i) a recursive DP-RTF feature estimation method, ii) an online multiple-speaker localization method, and iii) an multiple-speaker tracking method. The resulting framework provides online speaker counting, localization and consistent tracking (i.e. preserving speaker identity over a track in spite of intermittent speech production).

Figure 1 :

 1 Figure 1: Flowchart of the three localization and tracking modules.

Figure 2 :

 2 Figure 2: Results of speaker localization and tracking for Recording #1 / Task #6. (a) Ground truth trajectory and voice activity (red for speaker 1, black for speaker 2). Intervals in the trajectories are speaking pauses. (b) Result for localization, i.e. GMM Weights. (f) Result for tracker. Black and red colors demonstrate a succesful tracking, i.e. continuity of the tracks despite of speech pauses.

Table 1 :

 1 Localization and tracking results.

	Localization	24.1	12.7	4.0	-
	Tracking	22.7	12.4	4.1	10

MD rate (%) FA rate (%) MAE ( • ) IDs

September 17-20, 2018, Tokyo, Japan


C Shen

D Brito

J Diaz

D Zhang

C Bordes

F Faucher

S Garambois

On the use of a pulsed-laser source in laboratory seismic experiments

Reproduction of large-scale seismic exploration at laboratory-scale with controllable sources is a promising approach that could not only be applied to study small-scale physical properties of the medium, but also contribute to significant progress in wave-propagation understanding and complex media imaging at exploration scale via upscaling methods. We seek to characterize the properties of a laser-generated seismic source for new geophysical experiments at laboratory scale. This consists in generating seismic waves by pulsed-laser impacts and measuring the displacement wavefield by laser vibrometry. Parallel 2D/3D simulations using Discontinuous Galerkin discretization method and analytic predictions have been done to match the experimental data.

Pulsed-laser source : General

Source stability and reproducibility

Distribution of first arrival seismic wave amplitudes u z measured by LDV at the epicenter of the 10 mm thick aluminum block. P d denotes the incident power density.

Pulsed-laser source : Seismogram All regime

Two major parameters control the regime :

1) The incident energy ;

2) The laser spot size Regime evolution with d under constant input energy. Since the spot size depends on d, P d varies. However, the observable regime cannot be determined unilaterally by P d . Considerable diffractions caused by the cavity disturbs heavily the first arrival time (2 nd Tomo). The first break based inversion algorithm doesn't take into account this complicated phenomenon, leading to noisy, unstable and hardly convergent results.

Seismogram and radiation patterns

Toward Full Waveform Inversion (FWI) ?

The idea of Full Waveform Inversion (FWI) is to perform a quantitative reconstruction of the physical parameters, [START_REF] Tarantola | Inversion of seismic reflection data in the acoustic approximation[END_REF][START_REF] Virieux | An overview of full-waveform inversion in exploration geophysics[END_REF]. It is based on an iterative minimization of the residuals, defined as the difference between the observations and simulations, in order to recover the medium parameters (i.e. velocity, density). The problem writes as

min m J (m) = 1 2 d -F (m) 2 , ( 1
)
where J is the cost function, d the observed seismogram and F (m) are simulations using an initial model m. We wish to employ the code that has been developed in [START_REF] Faucher | Contributions to seismic full waveform inversion for time harmonic wave equations : stability estimates, convergence analysis, numerical experiments involving large scale optimization algorithms[END_REF], in particular starting with an initial model obtained from tomography. We illustrate in the figure below the observed and simulated traces for a single source. The deployment of iterative minimization would increase the accuracy of the velocity reconstruction, and allow a more precise characterization of the core interior.

Conclusion :

This laser-generated seismic source opens new perspectives on various applications such as precise Non-Destructive Test on metals under high temperatures [START_REF] Scruby | Non-contact ultrasonic measurements on steel at elevated temperatures[END_REF], micro-seismic exploration on small and intermediate scale samples of random shapes in the laboratory etc. We are especially interested in its eventual geophysical applications in rock mechanics [START_REF] Dupuy | Estimation of rock physics properties from seismic attributes ?part 2 : Applications[END_REF], rocks or digital rocks imaging for geological reservoirs explorations. The laser-impulsed source appears to be well controllable, flexible and reproducible under some precautions. The combination of the pulsed-laser source and the LDV is particularly adapted to generate broadband seismic full wavefields in heterogeneous natural rocks. The punctual feature of the source makes it convenient to model and implement into numerical schemes. We aim at performing FWI with these data and further more, doing amplitude/frequency related analyses for anisotropy [START_REF] Audoin | Measurement of stiffness coefficients of anisotropic materials from pointlike generation and detection of acoustic waves[END_REF], attenuation and poroelasticity etc.

  Lab set-up for pulsed-laser source characterization Two different experimental tools were mounted to investigate in Cartesian coordinates or in cylinder coordinates. 1 : Q-Switched laser generator ; 2 : convergent lens ; 3 : Aluminium cuboid samples of various thickness (10, 50, 100 mm), V p ≈6350 m/s, V s ≈V p /2 ; 4 : single point Laser Doppler Vibrometer (LDV) ; 5 core sample ; 6 piezoelectric source. laser-generated seismic source under thermoelastic regime with displacements computed by combining the wave equation and the thermoelastic equations [2, 3, 5]. b) : laser-generated seismic source under ablation regime which is modelled as a point-source [1, 2]. Zoom on the measured ablation-regime pulse Temporal laser-generated source S 0 retrieved from a record on the 50 mm thick aluminum block. b) : Frequency spectrum S 0 of the same source. c) First order time derivative of S 0 , thus showing the surface vibration velocities. d) Frequency spectrum of S ′ 0 .

  Seismogram measured along linear receivers on the 50 mm thick aluminum block with different sources, accompanied by simulated/modelled radiation patterns. a) : radiation pattern of a point source (laser ablation) ; b) : radiation pattern of a piezoelectric source (Φ ≈ 10mm) ; c) : radiation pattern of a typical thermoelastic source[START_REF] Bernstein | Line source representation for laser-generated ultrasound in aluminum[END_REF] under laser irritation.Epicentral records under the ablation regimeOverview on single experimental/simulated seismic traces (normal components). 2D/3D simulations are done by FE schemes featuring Discontinuous Galerkin (DG) with Interior Penalty (IP).

First

  ( el (f the b(tt(m secti(n (T(m(TV by J.Virie-. et al.)


Jessica Vial
email: jessicavial@mpimet.mpg.de

Christophe Cassou

Francis Codron

Sandrine Bony

Yohan Ruprich-Robert

Influence of the Atlantic meridional overturning circulation on the tropical climate response to CO 2 forcing

Internal multidecadal variability is an important modulator of the tropical climate response to CO 2 forcing • By impacting global atmospheric energetics, AMOC fluctuations affect the strength of the mean atmospheric tropical circulation • The dependence of the CO 2 -forced circulation response to AMOCrelated background state follows a "strong-gets-weaker" pattern

Introduction

Predicting the precipitation response to increasing greenhouse gases (GHG) atmospheric concentration remains a great scientific challenge. The need to better anticipate future precipitation changes is particularly strong in the tropics where nearly 50% of the Earth's population is influenced by rainfall and monsoon variability [START_REF] De Carvalho | The monsoons and climate change: observations and modeling[END_REF]. Past investigations have led to considerable progress in our understanding of the key factors controlling the long-term (centuries ahead) precipitation forced-response assessed from coupled models under specific GHG emission scenarios [START_REF] Chou | Mechanisms of global warming impacts on regional tropical precipitation[END_REF]Held and Soden, 2006;Chou et al., 2009;Seager et al., 2010;[START_REF] Xie | Global warming pattern formation: Sea surface temperature and rainfall[END_REF]Bony et al., 2013;[START_REF] Chadwick | Spatial patterns of precipitation change in CMIP5: Why the rich do not get richer in the tropics[END_REF]. Physical mechanisms have also been proposed to explain the considerable spread among the models' projections [START_REF] Ma | Regional patterns of sea surface temperature change: A source of uncertainty in future projections of precipitation and atmospheric circulation[END_REF][START_REF] Xie | Towards predictive understanding of regional climate change[END_REF]Long et al., 2016;[START_REF] Oueslati | Interpreting the inter-model spread in regional precipitation projections in the tropics: role of surface evaporation and cloud radiative effects[END_REF]. On the other hand, few studies so far have evaluated the importance of the decadal-to-multidecadal internal variability in modulating the long-term response of regional rainfall to anthropogenic forcing [START_REF] Xie | Towards predictive understanding of regional climate change[END_REF][START_REF] Cassou | Decadal climate variability and predictability: challenges and opportunities[END_REF], despite the dominant contribution of natural climate variability to uncertainties in regional precipitation projections [START_REF] Hawkins | The potential to narrow uncertainty in projections of regional precipitation change[END_REF].

This study aims at exploring and understanding the weight/fingerprint of the multidecadal internal variability in the estimation of GHG-forced precipitation changes in the tropics using the CNRM-CM5 atmosphere-ocean coupled model. Abrupt doubling CO 2 experiments have been carried out to extract the forced response. Multi-member ensembles are built here to evaluate the sensitivity (or interaction) of the forced response to (with) the variability of the Atlantic Meridional Overturning Circulation (AMOC). In CNRM-CM5, AMOC fluctuations are associated with surface ocean fin-gerprint, which projects on the so-called Atlantic Multidecadal Variability (AMV) pattern. In this model, the AMV corresponds to the leading mode of surface temperature internal low-frequency variability at global scale. Spectral analysis shows maximum energy concentration within the 80-120yr band as assessed from the long 850-yr preindustrial control (piControl) experiment [START_REF] Ruprich-Robert | Combined influences of seasonal east atlantic pattern and north atlantic oscillation to excite atlantic multidecadal variability in a climate model[END_REF] produced for the fifth phase of the Coupled Model Intercomparison Project (CMIP5, [START_REF] Taylor | An overview of CMIP5 and the experiment design[END_REF]).

In the following, abrupt CO 2 -doubling experiments are initialized from two extreme phases of the model AMOC and the impact of the initial ocean state on the tropical precipitation forced-response is analyzed by contrasting the two distinct ensembles. Each of them includes 10 members differing by their first day atmospheric states so that the full range of uncertainties associated with internal variability (ocean + atmosphere) is covered -because of the combination of both macro and micro perturbations for the ensemble generation, following [START_REF] Hawkins | Irreducible uncertainty in near-term climate projections[END_REF] terminology. This protocol contrasts to most previous studies, which either combine single model realization from several models [START_REF] Vecchi | Global warming and the weakening of the tropical circulation[END_REF]Long et al., 2016;[START_REF] Oueslati | Interpreting the inter-model spread in regional precipitation projections in the tropics: role of surface evaporation and cloud radiative effects[END_REF] or which assess uncertainties through the sole atmospheric state sensitivity [START_REF] Selten | Tropical origins for recent and future Northern Hemisphere climate change[END_REF][START_REF] Deser | Uncertainty in climate change projections: the role of internal variability[END_REF][START_REF] Kang | Uncertainty in climate change projections of the Hadley circulation[END_REF].

The rationale for the choice of selecting two contrasted states of AMOC internal variability from a given model relies on the fact that observational records over the historical period (on which anthropogenic forcing is active) are too short to know the actual stage of AMOC with respect to its own long-term variability. This is an irreducible source of uncertainties for the forced response, which should be quantified and better understood.

The rest of the paper is outlined as follows. The model simulations and methodology are described in Section 2. Following an overview of the AMOC fingerprints in the CNRM-CM5 Model (Section 3), we analyze the sensitivity of the projected changes in tropical precipitation and circulation to the strength of the AMOC (Section 4), and propose a physical interpretation (Section 5). Conclusions and discussion are given in Section 6.

Materials and method

Climate model simulations

We use the CNRM-CM5.1 Coupled General Circulation Model developed for CMIP5. The atmospheric component is the spectral ARPEGE-Climat-5.2 model, with a Gaussian grid of ∼1.4 • horizontal resolution and a vertical discretization over 31 levels (including 26 levels in the troposphere). The ocean model is NEMO (Nucleus for European Modeling of the Ocean) on the grid of ∼1 • horizontal resolution and a vertical discretization over 42 levels (including 18 levels in the upper 250-m mixed layer ocean). See [START_REF] Voldoire | The CNRM-CM5.1 global climate model: description and basic evaluation[END_REF] for a complete description of the model.

Oceanic conditions of Jan 1 st of year 303 and 141 are selected from the 850-year piControl experiment of CNRM-CM5 for the production of two unforced ensembles of 30 years and 10 members each. These two dates cor-respond to extreme opposite states of the AMOC in piControl (-4.1 and +2.8 standard deviation, respectively) associated with very pronounced negative and positive phases of AMV, respectively [see Fig. 4 in Ruprich-Robert and Cassou, 2015]. Individual members of each ensemble (hereafter referred to as Pcold and Pwarm) differ only by their atmospheric initial conditions, chosen randomly [see [START_REF] Ménégoz | Role of the Atlantic Multidecadal Variability in modulating the climate response to a Pinatubo-like volcanic eruption[END_REF], for a complete description of the protocol]. Two twin ensemble experiments (i.e., identical initial conditions) are conducted (hereafter referred to as Fcold and Fwarm) with an instantaneous doubling of atmospheric CO 2 (abrupt2xCO2 relative to preindustrial conditions using CMIP5 nomenclature) imposed at the beginning of the simulations and then hold fixed.

Calculations and notations

Here are the calculations and notations of the different metrics applied on annually-averaged fields from the ensembles:

• ∆X is the projected climate response to CO 2 abrupt doubling for a given ensemble (X F -X P ). When applicable, we evaluate both the 30-year averaged change and the contribution arising from the 1 styear only [as an approximation of the direct CO 2 -forced fast response, Bony et al., 2013].

• δX stands for the anomaly of a given member with respect to the two-ensemble mean. ∂X is the anomaly expressed in percentage, or fractional difference (relative to the two-ensemble mean).

• The statistical significance of AMOC-related ensemble mean differences (based upon p < 0.10, unless otherwise stated) is evaluated against a null hypothesis of zero change using a 2-sided Student's ttest.

• I refers to as the tropical atmospheric overturning circulation index defined by Bony et al. [2013]. I ≡ ω 500 ↓ -ω 500 ↑ with upward (ω 500 ↑< 0) and downward (ω 500 ↓> 0) pressure vertical velocities at 500 hPa averaged over 35 When applicable, we compute W over land + ocean or just ocean areas.

• The top of the atmosphere (TOA) instantaneous radiative forcing (RF ) quantifies the direct radiative impact of the doubling CO 2 concentration alone. It is commonly diagnosed from offline radiative transfer simulations [e.g., [START_REF] Collins | Radiative forcing by well-mixed greenhouse gases: Estimates from climate models in the Intergovernmental Panel on Climate Change (IPCC) Fourth Assessment Report (AR4)[END_REF], but due to the large number of ensemble experiments considered, we use the simple linear (least squares) regression model developed by [START_REF] Huang | Inhomogeneous radiative forcing of homogeneous greenhouse gases[END_REF]:

RF = RF + i A i y i -y i y i
, where A i is the regression coefficient of a dependent variable y i with global mean multiyear mean y i , and RF , the intercept of the regression model, represents the global mean of RF . This regression model accounts for several dependent variables i, which are taken from the P-ensemble simulations, including the temperature at the surface and at 10 hPa, the vertically-integrated water vapor and the cloud radiative forcing (difference between all-sky and clear-sky radiation fluxes at the TOA). See [START_REF] Huang | Inhomogeneous radiative forcing of homogeneous greenhouse gases[END_REF] for more details. • R a refers to the column-integrated atmospheric radiative fluxes (R a = R T OA -R surf ace , positive fluxes upward). The 1 st -year change in R a is used as an approximation of the direct CO 2 -forced radiative heating in the atmosphere.

AMOC fingerprints in CNRM-CM5

Figure 1 illustrates the yearly evolution of AMOC and Global Mean Surface Temperature (GMST) for the two 10-member P and F ensembles over the 30 years of integration. The model exhibits significant and persistent initial value predictability for both quantities as expected from spectral properties and physical mechanisms of AMOC variability assessed from piControl [START_REF] Ruprich-Robert | Combined influences of seasonal east atlantic pattern and north atlantic oscillation to excite atlantic multidecadal variability in a climate model[END_REF]. The 10-member Pcold and Pwarm ensembles for the AMOC time-series are clearly disjoint throughout the entire 30 years of simulations (Fig. 1a). This is a fortiori true for AMOC ensemble means, as well as for GMST, with a difference of ∼2.5 Sv and 0.11 • C, respectively.

Regionally, the AMOC fingerprints in CNRM-CM5 are consistent with previous inferences from observations and model-based studies [see Buckley and Marshall, 2016, for a review]. Difference between Pcold and Pwarm ensembles show a significant (marginal) surface and tropospheric cooling (warming) over the northern (southern) hemisphere when AMOC is weak due to reduced northward ocean heat transport. A southward displacement of the Hadley circulation and ITCZ-related precipitation from their climatological position is obtained consistently with the modification of the mean inter-hemispheric gradient (Supplementary Fig. 1 and [START_REF] Ménégoz | Role of the Atlantic Multidecadal Variability in modulating the climate response to a Pinatubo-like volcanic eruption[END_REF] for further description). A contraction and reinforcement of the Walker cell branches in the Pacific is also noticeable. Therefore, the AMOC-related anomalies are planetary and persistent so that the two ensembles can be treated as samples for two contrasted mean background states set by pure internal variability.

AMOC declines steadily in response to abrupt CO 2 forcing, but the two ensembles clearly remain disjoint over the 30-yr integration (Fig. 1c). Rapid warming occurs in both F-ensembles, which stay statistically distinguishable in GMST up to ∼15 years (Fig. 1d). These emphasize on the very strong persistence properties of AMOC in the model despite extreme forcing perturbation.

Sensitivity of precipitation and circulation forced-responses to AMOC states

The tropical precipitation and tropospheric circulation responses to CO 2 forcing and their modulation by different AMOC-related background states are presented in Figure 2. The CNRM-CM5 model exhibits increased rainfall in most of the deep tropics and enhanced drying over subtropical regions, similar to previous CMIP3/CMIP5 assessments [e.g., Seager et al., 2010;[START_REF] Xie | Towards predictive understanding of regional climate change[END_REF]Long et al., 2016;[START_REF] Oueslati | Interpreting the inter-model spread in regional precipitation projections in the tropics: role of surface evaporation and cloud radiative effects[END_REF]. Focus is laid on the tropical Pacific region, as it is where the fingerprint of multidecadal internal variability in the forced responses is the most important. Differences in rainfall changes between the cold and warm AMOC-related states are in the same order of magnitude (Fig. 2c) as the individual CO 2induced response assessed from ensemble-mean changes (Fig. 2a,b). This clearly indicates that AMOC is an important modulator of the projected forced-changes in this region. Comparison of Figures 2c and2f shows that the impact of different AMOC states on precipitation forced-responses is strongly related to circulation changes (pattern correlation equals to 0.95). A stepwise analysis of the vertically-integrated water budget equation further demonstrates that the mean circulation changes alone can provide a full accounting of the difference in the equatorial Pacific rainfall changes between the cold and warm ensembles (Supplementary Fig. 2 to 4). More- (df) same but for pressure vertical velocity at 500 hPa. (g -i) fractional forced changes (relative to P-ensembles) in Walker circulation as a function of its mean intensity in P-ensembles for warm (red) and cold (blue) ensembles: (g) 30-year mean change per degree of global surface warming, (h) 1 st -year change and (i) last 10-year mean change per degree of global surface warming. Each individual member is represented by a filled circle and the ensemble means by asterisks.

over, a decomposition of Figure 2f into its zonally averaged and zonally asymmetric component suggests that most of the modulation in rainfall changes originates from differences in the Walker circulation changes rather than in the Hadley circulation changes (Supplementary Fig. 5). The overall weakening of the Walker circulation (defined in Section 2.2) is amplified by a factor of two in the cold ensemble (Fig. 2g). This occurs almost everywhere across equatorial Pacific, but more particularly over the ocean (Fig. 2d -f). We further show that the direct CO 2 -forced Walker circulation response (Fig. 2h), rather than the temperature-mediated transient circulation response (Fig. 2i), explains to a large extent the sensitivity of precipi-tation forced-changes to AMOC states. Figure 2h shows a significant initial value predictability of the direct CO 2 -forced Walker circulation responses. The large range of projected changes in Walker circulation depends on its climatological strength: the Walker circulation is more strongly reduced when it is stronger a priori in the preindustrial control climate state (Fig. 2h). This relationship is even stronger when W is computed over ocean areas only (correlation R = 0.7 -not shown).

In the next section, we investigate the physical mechanism controlling these contrasted circulation responses to CO 2 forcing. To this end, we address two questions: (i) How do different AMOC-related states influence the tropical mean circulation in the control climate state? (Section 5.1); (ii) Why does the direct circulation response to increased CO 2 depend on the mean background state? (Section 5.2).

Physical understanding

AMOC-related influence on tropical circulation in control climate

Global constraints on precipitation and water vapor are commonly used to understand the tropical circulation response to global warming [e.g., [START_REF] Mitchell | On CO2 climate sensitivity and model dependence of results[END_REF]Held and Soden, 2006;[START_REF] Vecchi | Global warming and the weakening of the tropical circulation[END_REF]. We here apply the same energetic constraints to explain the sensitivity of the tropical mean overturning circulation to the AMOC state, since the latter has a global influence on the modeled climate (Fig. 1b and Supplementary Fig. 1) and thus justifies the use of such a conceptual framework.

When the AMOC is weaker in control experiments, the surface temperature decreases significantly by 0.11K (globally) and by 0.04K (in the tropics). As shown in Figure 3, this decrease in temperature is also associated with (i) a decrease in precipitation of 1.6%/K (Fig. 3a) that results from a weaker radiative cooling and surface evaporation (not shown) and (ii) a decrease in column water vapor amount following the Clausius-Clapeyron constraint of 7%/K (Fig. 3b). Following [START_REF] Mitchell | On CO2 climate sensitivity and model dependence of results[END_REF], this differential response of the tropical (or equivalently global) mean decrease in precipitation and humidity to surface cooling can be used to estimate the fractional increase of the tropical overturning circulation in response to surface cooling. Using δM/M = δP/P -0.07δT s /T s [Held and Soden, 2006;[START_REF] Vecchi | Global warming and the weakening of the tropical circulation[END_REF], where M is the convective mass flux, P is the precipitation and T s the surface temperature, we expect an overall strengthening of the convective mass fluxes by 5.4%/K in the cold ensemble (Fig. 3c). The relationship between this theoretical strengthening in convective mass flux and the actual change in the overturning circulation simulated by the model is fairly good with ensemble means lying right on top of the expected line set by the theory (Fig. 3d). This confirms that the influence of AMOC on the tropical mean circulation corresponding to pure by-products of internal climate processes can be very well explained by the latter conceptual framework. The AMOC-related change in the tropical overturning circulation is also reflected in the Hadley and Walker circulations, which are significantly reinforced when AMOC is weaker (Fig. 3e,f). Note that other factors control the Hadley and Walker circulations, influencing the exact magnitude of their AMOC-forced intensification and their relationships to the tropical mean overturning circulation. These factors could result from a redistribution of the convective and subsiding areas, which often accompanies the change in the intensity of the tropical overturning circulation [e.g., [START_REF] Vecchi | Global warming and the weakening of the tropical circulation[END_REF] or in the northward ocean heat transport [e.g., L'hévéder et al., 2015]. When AMOC is weak, the overall strengthening of the Walker and Hadley circulations is here associated with an eastward shift of convection over the central equatorial Pacific (Supplementary Fig. 1d) and a southward shift of the Hadley circulation, resulting in a strengthening (weakening) of the northern (southern) cell (Supplementary Fig. 1e). De-spite this phenomenon, our analysis suggests that the global constraints on precipitation and water vapor are not only useful to understand the mean circulation weakening with global surface warming, but also the relationship between the multidecadal variations of the global mean surface temperature and of the tropical overturning circulation (Hadley/Walker) associated with the AMOC internally generated fluctuations. These results are robust across a large ensemble of CMIP5 models, except maybe for the AMOCrelated variations in the Hadley circulation, which appear more model-dependent (Supplementary Text 2 and Fig. 6).

Dependence of the direct CO 2 -forced response on the background state

We now explore the reasons why the direct effect of CO 2 on the Walker circulation depends on the climatological strength of the Walker circulation in the preindustrial climate state, and thus on the state of the AMOC (cf. Fig. 2h). One mechanism, proposed by [START_REF] Merlis | Direct weakening of tropical circulations from masked CO 2 radiative forcing[END_REF], relates the direct weakening of the overturning circulation to the spatial structure of the instantaneous radiative forcing to increased CO 2 . The argument is based on the moist energetics of the tropical atmosphere, which constrain the mean overturning circulation to export energy from energy-rich (warm/moist) regions toward relatively energy-poor (cold/dry) regions [START_REF] Neelin | Modeling tropical convergence based on the moist static energy budget[END_REF][START_REF] Trenberth | The flow of energy through the Earth's climate system[END_REF].

As shown by [START_REF] Huang | Inhomogeneous radiative forcing of homogeneous greenhouse gases[END_REF], the spatial structure of the instantaneous radiative forcing strongly depends on the horizontal gradients of atmospheric and surface variables of the control climate state, especially surface and atmospheric temperature, atmospheric humidity and clouds. Figure 4a shows that the CO 2 radiative forcing (as estimated by Huang et al.

[2016] -see also Section 2.2) is stronger in large-scale subsidence areas and weaker in regions of large-scale rising motions (compare RF with black contours in Fig. 4a). This dependence of the forcing on the mean circulation is explained by the presence of deep convective clouds and strong water vapor mixing ratios in regions of large-scale ascent, which masks the effect of increased CO 2 on TOA radiative fluxes [START_REF] Merlis | Direct weakening of tropical circulations from masked CO 2 radiative forcing[END_REF][START_REF] Huang | Inhomogeneous radiative forcing of homogeneous greenhouse gases[END_REF][START_REF] Huang | On the pattern of CO2 radiative forcing and poleward energy transport[END_REF][START_REF] Xia | Differential Radiative Heating Drives Tropical Atmospheric Circulation Weakening[END_REF]. [START_REF] Trenberth | The flow of energy through the Earth's climate system[END_REF]. Over the tropical oceans, much of the instantaneous radiative forcing is confined to the atmosphere [e.g., [START_REF] Huang | On the pattern of CO2 radiative forcing and poleward energy transport[END_REF][START_REF] Xia | Differential Radiative Heating Drives Tropical Atmospheric Circulation Weakening[END_REF]. Therefore, the effect of the spatially inhomogeneous radiative forcing, shown in Fig. 4a, is to radiatively heat the atmosphere more strongly in the subsiding than ascending region (as suggested in Fig. 4b by the 1 st -year change in R a ), which then requires the circulation to weaken to restore local atmospheric energetics. This also suggests that when the tropical overturning circulation is stronger in the control climate state, cloud and water vapor gradients between convective and subsidence areas are reinforced, and the larger spatial heterogeneity of the CO 2 forcing potentially leads to a stronger weakening of the circulation.

Consistent with this interpretation and because the Walker circulation is stronger when the AMOC is weaker (in the cold ensemble), the TOA ra- diative forcing and 1 st -year atmospheric radiative heating are found to be (i) stronger in the eastern equatorial Pacific, where the subsiding branch of circulation is stronger, and (ii) weaker in the central equatorial Pacific region, where the rising branch of circulation is stronger (Fig. 4 c,d and4 e). Therefore, the stronger Pacific equatorial zonal gradient in radiative forcing that is produced in the cold ensemble is expected to induce a stronger weakening of the Walker circulation by the direct effect of CO 2 on the atmospheric radiative fluxes (as suggested in Figure 4 f).

Conclusion

In this study, the tropical precipitation and tropospheric circulation responses to CO 2 forcing and its modulation by different states of the Atlantic Multidecadal Overturning Circulation (AMOC) have been investigated using multi-member ensembles from the fully coupled CNRM-CM5 General Circulation Model. A very large range of possible precipitation projections is obtained in the equatorial Pacific region due to pure internal variability, assessed here through an ensemblist approach and a combination of macro and micro perturbations for its generation. The associated uncertainties are comparable to -or even larger than -the forced response estimated through ensemble-mean changes. This underscores the crucial need for realizing multi-member ensembles (by sampling both atmospheric and oceanic initial conditions) to better characterize the forced model's climate change projections, even for extreme forcing perturbations as is the case here.

Contrasted precipitation and circulation forced-responses arise depending on the strength of the AMOC in the preindustrial climate state, used as initial states. When the AMOC is weaker, the slowdown of the Walker circulation in response to increased CO 2 is amplified (compared to when the AMOC is stronger), which substantially influences the amplitude of rainfall changes over Indonesia, southeast Asia and the equatorial Pacific ocean. The direct CO 2 -forced effect alone on circulation changes explains the overall sensitivity of precipitation and circulation forced-changes to AMOC states. A mechanism is proposed to explain how the direct weakening of the tropical circulation in response to CO 2 depends on AMOC-related background climate states. We show that:

1. Global energetic constraints on precipitation and water vapor explain the impact of different AMOC states on the mean tropical circulation in the control climate. When the AMOC is weak, a strengthening of the tropical overturning circulation (Hadley and Walker circulations) is obtained consistently with a significant decrease in global mean surface temperature. 2. The direct weakening of the circulation in response to increased CO 2 results from the dependence of the CO 2 radiative forcing on the spatial distribution of clouds and water vapor [START_REF] Merlis | Direct weakening of tropical circulations from masked CO 2 radiative forcing[END_REF]. Therefore, the stronger overturning circulation (when AMOC is weaker), which enhances the cloud and water vapor contrasts between convective and subsidence areas and thus the spatial heterogeneity of the CO 2 forcing, leads to a stronger direct weakening of the circulation.

Our interpretation provides an example of a new mechanism by which the response of tropical precipitation to CO 2 forcing depends on the background climate state. We call it the "strong-gets-weaker" mechanism, as it essentially operates through the circulation changes: the stronger the circulation in the control climate state, the weaker it gets in response to CO 2 forcing, enhancing the amplitude of rainfall changes at the regional scales. Although the example given here explains the AMOC influence on Walker circulation forced changes, this mechanism might also provide the backdrop for any factor that influences the intensity of the mean tropical circulation in the control climate state. A point that should motivate further investigation for understanding the inter-model spread in CO 2 -forced circulation changes.

These insights are based on the analysis of abrupt doubling CO 2 experiments; the question arises as to how relevant they are to interpret the AMOC influence on the climate response to gradually increasing CO 2 concentration. We analyzed twin ensemble transient climate change simulations using the exact same initial conditions, in which CO 2 increases at 2% per year until doubling the initial CO 2 concentration. We found very similar ensemble-mean sensitivities of precipitation and circulation forced-changes to AMOC states (Supplementary Fig. S7 -S14). Moreover, strong differences between respective members (from abrupt and transient CO 2 forcing experiments) are noticeable, despite the fact that they have the exact same initial conditions. This confirms the robustness of our results, which emphasize the very important role of the background climate state for the projected tropical and circulation changes, at least in the CNRM-CM5 model.

Contents of this file

1. Text S1 to S2 2. Figures S1 to S14 3. Table S1 D R A F T July 4, 2018, 2:42pm D R A F T
Text S1. Origin of AMOC-related differences in tropical precipitation changes.

The impact of different AMOC-related background states on the 30-year mean CO 2 -forced changes in tropical precipitation is here analysed using the monthly-mean vertically integrated water budget equation [e.g., Chou et al., 2009;Seager et al., 2010;Bony et al., 2013;Long et al., 2016] :

∆P = ∆E -∆ [ωΓ q ] -∆ V . ∇q + res ( 1
)
where P is the surface precipitation (mm d -1 ), E the surface evaporation (mm d -1 ), ω the vertical pressure velocity profile (Pa s -1 ), V the horizontal wind profile (m s -1 ), q the specific humidity profile (×L v , in J kg -1 ), Γ q its vertical gradient (J kg -1 Pa -1 ) and ∇q its horizontal gradient (J kg -1 m -1 ) and res is the residual term owing to sub-monthly transient eddy convergence. Square brackets refer to mass-weighted vertical integrals and ∆ refers to 30-year CO 2 -forced changes.

Figure S2 shows that the change in vertical moisture advection (∆ [ωΓ q ] in Fig. S2d) is clearly the dominant source of AMOC-related differences in the tropical precipitation change, with negligible contributions from the other three terms. We further decompose the vertical advection term into a thermodynamic and a dynamic component to diagnose the relative contributions of moisture and circulation changes, respectively [e.g., Seager et al., 2010;Long et al., 2016] : ∆ [ωΓ q ] = ∆P ther + ∆P dyn , where ∆P ther = [ω∆Γ q ] and ∆P dyn = [Γ q ∆ω].

Figure S3 shows that all of the difference in tropical precipitation changes between the cold and warm ensembles is due to ∆P dyn (Fig. S3b).

The separate influence of the AMOC state on the preindustrial moisture profile and on the D R A F T July 4, 2018, 2:42pm D R A F T forced changes in the mean circulation is then assessed by expressing the dynamical component as the product of surface specific humidity in unperturbed (P-ensemble) climate simulations (q sf c ) and CO 2 -forced change in mid-tropospheric vertical motion (∆ω 500 ) :

∆P dyn -1 g q sf c ∆ω 500 . This has been shown to be a good approximation for the Tropics [Held and Soden, 2006;[START_REF] Huang | Patterns of the seasonal response of tropical rainfall to global warming[END_REF][START_REF] Huang | Regional response of annual-mean tropical rainfall to global warming[END_REF]Long et al., 2016;[START_REF] Wills | How stationary eddies shape changes in the hydrological cycle: Zonally asymmetric experiments in an idealized GCM[END_REF]; here, the correlation across all experiments between ∆P dyn and -1 g q sf c ∆ω 500 in the two main areas where the difference between the cold and warm ensembles is strong (the Niño 3 and Niño 4 regions) is close to unity.

Figures S4 shows that ∆ω 500 accounts for most of the difference in the tropical rainfall changes between the cold and warm ensembles (Figures S4b) and of the variance in ∆P across the warm and cold ensembles (Fig. S4d andS4f). Therefore, it is clear that the influence of the AMOC on the change in tropical precipitation operates essentially through the modulation of the mean circulation changes, and especially of its zonally asymmetric component (the Walker circulation, Fig. S5c).

Text S2. Robustness of the relationship between the low-frequency variations of the AMOC and of the tropical atmospheric circulation

In the CNRM-CM5 multi-member ensembles that we use for this study, we compare two extreme states of the AMOC that are associated with pronounced and persistent negative/positive phases of the AMV, which in turn affect the climate at the global scale. However, the spatio-temporal properties of the simulated AMOC and AMV, the relationship between D R A F T July 4, 2018, 2:42pm D R A F T them and their influence on climate at the larger scale can vary substantially between models [e.g., [START_REF] Zhang | Multidecadal North Atlantic sea surface temperature and Atlantic meridional overturning circulation variability in CMIP5 historical simulations[END_REF][START_REF] Ba | A multi-model comparison of Atlantic multidecadal variability[END_REF][START_REF] Qasmi | Teleconnection Between Atlantic Multidecadal Variability and European Temperature: Diversity and Evaluation of the Coupled Model Intercomparison Project Phase 5 Models[END_REF][START_REF] Yan | Underestimated AMOC variability and implications for AMV and predictability in CMIP models[END_REF]. Moreover, in most models the persistence, amplitude and large-scale spatial coherence of the AMOC and AMV are underestimated with respect to observations [e.g., [START_REF] Qasmi | Teleconnection Between Atlantic Multidecadal Variability and European Temperature: Diversity and Evaluation of the Coupled Model Intercomparison Project Phase 5 Models[END_REF][START_REF] Yan | Underestimated AMOC variability and implications for AMV and predictability in CMIP models[END_REF] -keeping in mind the large uncertainties also in the observations due to sampling and the different existing methodologies for extracting externally forced variations [START_REF] Qasmi | Teleconnection Between Atlantic Multidecadal Variability and European Temperature: Diversity and Evaluation of the Coupled Model Intercomparison Project Phase 5 Models[END_REF].

The CNRM-CM5 model is no exception to these systematic biases, although it lies above the multi-model mean for the persistence and amplitude of the AMV [START_REF] Qasmi | Teleconnection Between Atlantic Multidecadal Variability and European Temperature: Diversity and Evaluation of the Coupled Model Intercomparison Project Phase 5 Models[END_REF] Fig. 1 and Supplementary Fig. S6] and tends to have a stronger hemispheric-scale temperature response to AMOC fluctuations compared to other models [Yan et al., 2018, their Fig. 3].

Given the large diversity in model behaviors, we here aim to provide some insights into the robustness of our findings across different coupled models, and especially regarding the relationship between the AMOC and the atmospheric tropical circulation in the control climate state (cf. Section 5.1), which is determining in the sensitivity of the CO 2 -forced tropical climate response to the AMOC state in the CNRM-CM5 model. More specifically, here we address the question of whether or not the global energetic framework (described in Section 5.1) can robustly constrain the response of the atmospheric tropical overturning circulation to the AMOC/AMV state in the control climate. To do so, we perform a similar analysis as in Figure 3, but by compositing the tropically-averaged climate variables and atmospheric circulation indices with respect to positive/negative phases of the AMV across 28 coupled model simulations using the long piControl experiments of the CMIP5 archive (cf. Table S1).

D R A F T July 4, 2018, 2:42pm D R A F T Among those models, upon availability of the ocean meridional overturning circulation data, 8

were also retained to assess the relationship with respect to the AMOC phases. The AMV index is defined as the detrended, low-pass filtered annual North Atlantic SST anomalies averaged over 0 • -60 • N and 0 • -80 • W [e.g., [START_REF] Qasmi | Teleconnection Between Atlantic Multidecadal Variability and European Temperature: Diversity and Evaluation of the Coupled Model Intercomparison Project Phase 5 Models[END_REF], using a 30-year low-pass filter (instead of the more traditional 10-year filter) to be more comparable with the 30-year averaged periods used in our main analysis. The AMOC index is similarly defined as detrended, 30-year low-pass filtered annual anomalies, but from the maximum meridional stream function deeper than 500 m at 26 • N [START_REF] Qasmi | Teleconnection Between Atlantic Multidecadal Variability and European Temperature: Diversity and Evaluation of the Coupled Model Intercomparison Project Phase 5 Models[END_REF].

The results presented in Figure S6 clearly show that the relationship between the phase of the AMOC/AMV and the strength of the atmospheric tropical circulation (I) is robust across CMIP5 models and well explained by the global constraints on precipitation and water vapor (Fig. S6a -d). This is also true for the Walker circulation, which is strongly related to the multidecadal variations of the tropical mean overturning circulation (Fig. S6f). However, the AMOC/AMV-related variations in the Hadley circulation appear to be more model-dependent (Fig. S6a -e). One possible explanation could be related to different sensitivities of the

Figure 1 :

 1 Figure 1: Yearly evolution of (a) AMOC and (b) GMST for Pwarm (orange) and Pcold (green) ensembles. Thin lines stand for individual members and thick lines for ensemble means.(c, d) same but for Fwarm (red) and Fcold (blue) ensembles, with the P-ensemble mean quantities reproduced for reference.

Figure 2 :

 2 Figure 2: (a -c) 30-year mean response to CO 2 forcing in precipitation for (a) cold and (b) warm ensemble means over tropical Pacific (30 • S-30 • N, 90 • E-70 • W) and (c) AMOC-related difference between cold and warm ensembles (hatching indicating statistical significance).(df) same but for pressure vertical velocity at 500 hPa. (g -i) fractional forced changes (relative to P-ensembles) in Walker circulation as a function of its mean intensity in P-ensembles for warm (red) and cold (blue) ensembles: (g) 30-year mean change per degree of global surface warming, (h) 1 st -year change and (i) last 10-year mean change per degree of global surface warming. Each individual member is represented by a filled circle and the ensemble means by asterisks.

Figure 3 :

 3 Figure3: (a -c) AMOC-related fractional differences (∂X in Section 2.2) for tropicallyaveraged (a) precipitation, (b) vertically-integrated water vapor and (c) convective mass flux -estimated as the difference between (a) and (b) as in[START_REF] Vecchi | Global warming and the weakening of the tropical circulation[END_REF] -as a function of AMOC-related absolute differences in tropically-averaged surface temperature (δT s , in K). (d -e) AMOC-related fractional differences for the tropically-averaged (d) overturning circulation, (e) Hadley circulation and (f) Walker circulation as a function of AMOC-related fractional differences in the estimated tropical convective mass flux (in d) and in the simulated overturning circulation (in e, f). All differences are significant at p<0.03. The color code is the same as Figure1for P-ensembles; warm (orange) and cold (green) members are shown by filled circles and their respective ensemble means by asterisks. Also reported in (a -c) is the slope of the least-square regression fits (solid lines). In (d) the dashed line highlights the one-to-one relationship between ∂I and ∂M .

Figure 4 :

 4 Figure 4: (a -b) The two-ensemble mean (a) instantaneous CO 2 radiative forcing (RF in Section 2.2, shaded) and (b) 1 st -year change in atmospheric radiative cooling (∆R a in Section 2.2, shaded); overlayed contours represent the two-ensemble mean 30-year mean atmospheric circulation (ω 500 ) for the P-ensembles. (c -d) AMOC-related difference for (c) RF and (d) ∆R a ; hatching indicating statistical significance and overlayed contours show the AMOC-related difference in ω 500 . (e) Zonal gradient in RF (∂ x RF in Section 2.2) as a function of Walker circulation in the P-ensembles. (f) 1 st -year fractional change (compared with P-ensembles) in Walker circulation as a function of ∂ x RF . W is computed over ocean areas only to ensure that the relationships are not influenced by land-sea contrasts. The legend is similar to Figure 2: warm (red) and cold (blue) members are shown by filled circles and their respective ensemble means by asterisks.

Figure S1 :Figure S3 :Figure S5 :Figure S6 :

 S1S3S5S6 Figure S1: AMOC-related difference between the cold and warm ensembles for (a) surface tempera-

Figure S7 :

 S7 Figure S7: 30-year abrupt CO 2 -forced change in precipitation (in mm d -1 ) over tropical Pacific for

Figure S8 :

 S8 Figure S8: Same as Figure S7 but for the warm ensemble: abrupt CO 2 -forced change in precipitation

Figure S9 :

 S9 Figure S9: Same as Figure S7 but for abrupt CO 2 -forced change in mid-tropospheric circulation (assessed with ω 500 , in hPa d -1 ) in the cold ensemble.

Figure S10 :

 S10 Figure S10: Same as Figure S8 but for the warm ensemble : abrupt CO 2 -forced change in ω 500 in

Figure S11 :

 S11 Figure S11: Same as Figure S7 but for transient CO 2 -forced experiments: change in precipitation

Figure S12 :

 S12 Figure S12: Same as Figure S8 but for transient CO 2 -forced experiments: change in precipitation

Figure S13 :

 S13 Figure S13: Same as Figure S9 but for transient CO 2 -forced experiments: change in ω 500 in the

Figure S14 :

 S14 Figure S14: Same as Figure S10 but for transient CO 2 -forced experiments: change in ω 500 in the

is the Hadley circulation index defined as I, but on zonally aver

  aged values of ω 500 over 35 • S-35 • N. • W is the Pacific Walker circulation index defined as the difference between ω 500 ↑ over the western/central Pacific (15 • S-15 • N, 110 • E-160 • W) and ω 500 ↓ over the eastern basin (15 • S-15 • N, 80 • -140 • W).

• S-35 • N. • H

  • δ x RF and δ x R a are zonal gradients of RF and R a over the equatorial Pacific estimated by the difference between (15 • S-15 • N, 110 • E-160 • W) and (15 • S-15 • N, 80 • -140 • W), as in W calculation.

Acknowledgments

This work is supported by the French National Research Agency (ANR) project MORDICUS (ANR-13-SENV-0002-01). Primary data and scripts used in the analysis and other supplementary information are archived at ftp://ftp.cerfacs.fr/pub/globc/exchanges/cassou/Data4Papers/ GRL_Vial_etal_2018/. We acknowledge the World Climate Research Programme's Working Group on Coupled Modelling, which is responsible for CMIP, and we thank the climate modeling groups (listed in Table S1 of this paper) for producing and making available their model output on the CEDA portal (http://browse.ceda.ac.uk). Additionally, we are grateful to the two reviewers for their helpful comments.

should motivate further research to better understand the differences in the Hadley/Walker circulation responses to low-frequency AMOC/AMV fluctuations.

X -23

Table S1: List of the CMIP5 coupled models used to produce the AMOC/AMV compositing analysis in Figure S6, together with the respective length of their picontrol simulations (in years


Mohamad Jouni

Mauro Dalla Mura

Pierre Comon

Some issues in computing the CP decomposition of NonNegative Tensors

Tensor decompositions are still in the process of study and development. In this paper, we point out a problem existing in nonnegative tensor decompositions, stemming from the representation of decomposable tensors by outer products of vectors, and propose approaches to solve it. In fact, a scaling indeterminacy appears whereas it is not inherent in the decomposition, and the choice of scaling factors has an impact during the execution of iterative algorithms and should not be overlooked. Computer experiments support the interest in the greedy algorithm proposed, in the case of the CP decomposition.

Introduction

Tensors of order d are represented by data arrays with d indices, (d = 2 for matrices). They provide unique features as they are a suitable data structure for representing multimodal or multisource data, in which each diversity is represented by one of the ways of the tensor. One of the most interesting applications of tensors is the Canonical Polyadic (CP) decomposition defined below, which aims at representing a tensor as a sum of decomposable rank one tensors, revealing relationships among its d ways.

CP decomposition. In this paper, we shall focus our attention on the CP decomposition of third order tensors. To begin with, a tensor D is decomposable if it can be expressed as the outer product of vectors, i.e.: D ijk = a i b j c k , which will be denoted compactly as D = a ⊗ b ⊗ c, where ⊗ is the outer (tensor) product. Next, every real tensor T of order d = 3 and dimensions I × J × K, admits a CP Decomposition of the following form:

T = R r=1 λ r D(r), (1)
where D(r) def = a(r) ⊗ b(r) ⊗ c(r), a(r), b(r) and c(r) being real vectors, which can be stored in the so-called factor matrices, A, B and C respectively, of size I × R, J × R, and K × R respectively, and λ r are real positive scalars. The CP decomposition reveals tensor rank when R is minimal, which will be assumed from now on; for instance, tensors D(r) are of rank 1. Note that another writing of (1) in terms of factor matrices is T ijk = R r=1 λ r A ir B jr C kr . In addition, because of the presence of λ r , the columns of factor matrices may be normalized to 1.

At this stage, it is important to stress that there is no scaling ambiguity in the CP decomposition [START_REF] Ruiz-Tolosa | From Vectors to Tensors[END_REF], contrary to what is sometimes claimed in the literature. Only the representation of tensors D(r) by triplets of vectors is subject to this indeterminacy. In fact, by definition, tensors are precisely equivalence classes with respect to scaling [START_REF] Ruiz-Tolosa | From Vectors to Tensors[END_REF][START_REF] Landsberg | Tensors: Geometry and Applications[END_REF][START_REF] Hackbusch | Tensor Spaces and Numerical Tensor Calculus[END_REF][START_REF] Comon | Tensors: a brief introduction[END_REF]: the triplets (a, b, c) and (αa, βb, γc) represent the same tensor provided that αβγ = 1.

The rank R of the CP decomposition ( 1) is of particular interest in applications since it is related to the intrinsic dimensionality of multilinear data. Furthermore, the CP decomposition, contrary to other tensor decompositions, e.g., Tucker's or High-Order Singular Value Decomposition (HOSVD), enjoys uniqueness if the rank is not too large [START_REF] Sidiropoulos | On the uniqueness of multilinear decomposition of N-way arrays[END_REF][START_REF]Three-way arrays: Rank and uniqueness of trilinear decompositions[END_REF][START_REF] Chiantini | An algorithm for generic and low-rank specific identifiability of complex tensors[END_REF]. Uniqueness is of utmost importance since it eventually allows physical interpretation of relationships among the ways of a tensor.

Motivation

Nonnegativity. When the observation tensor T contains only real nonnegative entries, it is suitable to impose decomposable tensors D(r) to also be nonnegative. By doing this, we define a nonnegative rank, R + , which may be larger than R. This is actually already true for matrices (tensors of order 2). In fact, Herbert E. Robbins exhibited a simple example of a 5 × 5 matrix having rank 3 but nonnegative rank 4; see [START_REF] Cohen | Nonnegative ranks, decompositions and factorizations of nonnegative matrices[END_REF][START_REF] Comon | Tensors: a brief introduction[END_REF] for its expression. It is thus necessary to define the nonnegative CP decomposition of a nonnegative tensor as:

T = R + r=1 λ r a(r) ⊗ b(r) ⊗ c(r), (2)
where a i (r) ∈ R + , b j (r) ∈ R + and c k (r) ∈ R + , ∀(i, j, k, r).

There are many applications where nonnegativity is relevant, as to provide better interpretable results when dealing with variables related to physical quantities such as luminance in images, spectra or chemical concentrations [START_REF] Cichocki | Nonnegative Matrix and Tensor Factorizations[END_REF][START_REF] Comon | Handbook of Blind Source Separation, Independent Component Analysis and Applications[END_REF]. There exist many algorithms aiming at computing the CP decomposition of nonnegative tensors [START_REF] Cichocki | Nonnegative Matrix and Tensor Factorizations[END_REF][START_REF] Huang | A flexible and efficient algorithmic framework for constrained matrix and tensor factorization[END_REF]. However, due to measurement noise or modeling errors, the tensor to decompose may not be nonnegative or may have a too large rank, hence requiring to be approximated. It turns out that, given any real tensor T of rank R, it is fortunately always possible to find a best nonnegative approximation of T of given nonnegative rank R + . This problem is indeed well-posed [START_REF] Lim | Nonnegative approximations of nonnegative tensors[END_REF][START_REF] Qi | Uniqueness of nonnegative tensor approximations[END_REF] (which would not be the case in R instead of R + ).

Projection onto the nonnegative orthant: In the nonnegative CP decomposition (2), all quantities are nonnegative. For instance, vector a(r) belongs to the nonnegative orthant (R + ) I . In iterative algorithms, this constraint is ensured at each iteration by projecting a computed value onto the nonnegative orthant. This is where the problem shows up. In fact, projecting D(r) or its building vectors {a(r), b(r), c(r)} do not yield the same result. Since this observation is already true for matrices, a simple example will be most convincing.

Example. Take the matrix M below, of rank 1. Now its projection M + has rank 2. So it is preferred to project its supporting vectors {a, b} instead. The obtained vectors are {a + , b + } and yield a matrix of nonnegative rank equal to 1:

M = 4 -2 -2 1 = 2 -1 ⊗ 2 -1 = a⊗b, M + = 4 0 0 1 , a + ⊗b + = 4 0 0 0 .
The problem is that vectors {a, b} are not uniquely defined. We could have taken {-a, -b} without changing M . Should we do that, we obtain instead:

M = 4 -2 -2 1 = -2 1 ⊗ -2 1 = a⊗b, M + = 4 0 0 1 , a + ⊗b + = 0 0 0 1 .
We see that the projected matrix of nonnegative rank 1 is not the same. This issue comes from the fact that no care has been taken of the scaling indeterminacies (which reduce to sign indeterminacies thanks to the use of factors λr) inherent in the representation of a rank-1 tensor by a triplet of vectors.

Algorithms resorting to projection include Alternating Nonnegative Least Squares (ANLS) [START_REF] Cichocki | Nonnegative Matrix and Tensor Factorizations[END_REF], Projected and Compressed ANLS (ProCo) [START_REF] Cohen | Fast decomposition of large nonnegative tensors[END_REF], or Alternating Direction Method of Multipliers (ADMM) [START_REF] Huang | A flexible and efficient algorithmic framework for constrained matrix and tensor factorization[END_REF], among others. Hard thresholding is the procedure in which it is the easiest to illustrate the occurring of the problem.

Algorithm 1 Alternating Nonnegative Least Squares (ANLS)

Require: T , B = B[0], C = C[0] t = 0; while stopping criterion is not met, do t = t + 1 compute A from B[t -1] and C[t -1]; A[t] ← A + compute B from C[t -1] and A[t]; B[t] ← B + compute C from A[t] and B[t]; C[t] ← C + normalize columns: A[t] ← A[t]Λ -1 A ; B[t] ← B[t]Λ -1 B ; C[t] ← C[t]ΛAΛB; end while normalize columns: C[t] ← C[t]Λ -1 return A[t], B[t], C[t], Λ
ANLS. One algorithm that has been widely used to compute CP decomposition (1) is the Alternating Least Squares (ALS) algorithm. ALS minimizes with respect to matrices A, B, C in an alternating fashion, the loss:

Φ = ijk [T ijk - R r=1 λ r A ir B jr C kr ] 2 .
(

) 3
Factor matrices are updated in turns during each iteration until a certain condition is attained (e.g. the number of iterations or a certain threshold on the reconstruction error). When a nonnegative decomposition is sought, each factor matrix can be projected onto the nonnegative orthant right after its calculation; this is the ANLS algorithm [9, p.47]. The pseudo-code is given in Alg.1.

Proposed approach

We illustrate the problem with hard thresholding (cf. Section 4), but our solution could also reveal useful in soft thresholding as well. The problem is worse when all entries in a column vector are set to zero; this prevents its normalization (as it would lead to a division by zero) or imposes an erroneous reduction of the rank (due to the arbitrary removal of the null columns). The solution we describe overcomes these two difficulties most of the time, up to negligible extraneous computation load. We propose to implement this in a procedure to be executed before projection. The concept goes as follows. Because of normalization, the scaling indeterminacy reduces merely to signs. In fact, in every decomposable tensor D(r), we have two variables, , η ∈ {-1, +1}, which are to be used as sign flippers for the columns a(r), b(r) and c(r) that are together involved in an outer product term, without changing the result of the outer product given by:

a(
Vectors a + (r), b + (r) and c + (r) are defined in a similar manner, with positive entries. In particular, a + (r) + a -(r) = a (r). Given a triplet of vectors, (a, b, c), there are 4 possibilities to construct a nonnegative decomposable tensor D [ ] by just flipping their signs without changing the value of (a ⊗ b ⊗ c) and by setting to zero negative values in each vector (refer to Table.1).

( ,η) (+,+) (-,-) (+,-) (-,+) 1 2 3 4 D [ ]
D [1] = D [2] = D [3] = D [4] where for the sake of convenience, a + stands for vector a [ ] (r) + , and similarly for a -, b + , b -, c + and c -.

= a + ⊗ b + ⊗ c + a + ⊗ b -⊗ c -a -⊗ b + ⊗ c -a -⊗ b -⊗ c + Table 1.
We are interested to know which combination would yield the minimal number of resets. Ultimately, we are concerned about (i) avoiding to set a whole vector to zero, which would lead to decrease the rank. This goal can mean "set as few entries to zero as possible". And we also aim at (ii) minimizing the distance between the original tensor and its nonnegative approximation.

We explored several criteria. The first is to minimize

Φ 0 = T -r D [ ] (r) 2 .
This criterion is very costly to optimize, due to the large number of combinations. In fact, for every r, there are 4 possibilities to assign ( , η), and this assignment can be different for each r. This would result in 4 R possibilities to explore. This is why we propose two greedy algorithms searching for the optimal solution D [ ] (r) independently for every r. One possibility is to minimize w.r.t. the following product for every r independently, and for the L 2 norm:

Φ 1 ( , r) = D(r) -D [ ] (r) 2 2 . (8)
Let us express this criterion for = 1, without loss of generality. We have for any fixed r (that we drop for the sake of convenience):

Φ 1 (1, r) = D(r) 2 2 + D [1] (r) 2 2 -2 ijk a i a + i b j b + j c k c + k . ( 9
)
The last term can be rewritten as 2 (a T a + )(b T b + )(c T c + ). Next, it is also equal to 2 a + 2 b + 2 c + 2 , since a + and a -are orthogonal and a = a + -a -. This suggests another criterion to minimize w.r.t. :

Φ 2 ( , r) = a -• b -• c - (10)
Criteria Φ 1 and Φ 2 are easy to optimize w.r.t. ( , η), i.e w.r.t. , and need negligible extraneous computation load.

A Toy Example

Consider the factor matrices:

; Φ1( ) ← D -D [ ] 2 2 ; end while Find o = arg min Φ1( ); a(r) ← ( o) η( o) a(r); b(r) ← ( o) b(r); c(r) ← η( o) c(r); return a(r), b(r), c(r); Algorithm 3 Modified ANLS Require: T , B = B[0], C = C[0] t = 0; while stopping criterion is not met, do t = t + 1 compute A from B[t -1] and C[t -1]; A[t] compute B from C[t -1] and A[t]; B[t] compute C from A[t] and B[t]; C[t] r = 0; while r < R do update a(r), b (r)
, and c(r) using Alg.2; r = r + 1; end while

A[t] ← A + ; B[t] ← B + ; C[t] ← C + ;
replace null columns in A, B or C by random values using e.g. absolute value of standard Gaussian distribution; (see Sect.5) normalize columns:

A[t] ← A[t]Λ -1 A ; B[t] ← B[t]Λ -1 B ; C[t] ← C[t]ΛAΛB; end while normalize columns: C[t] ← C[t]Λ -1 return A[t], B[t], C[t], Λ
When computing the CP Decomposition with Alg.1, after one update of A, one of its columns became negative, and hence one of its columns got discarded, and the rank was decreased by 1. Note that, for the sake of conciseness, during the loop of updates only the columns of A and B are normalized and their norms multiply C; after the loop ends, C is normalized and its column norms (containing that of A and B) form the values of Λ (cf. Alg.1). As can be seen in Fig. 1, 107 realizations out of 500 are unsuccessful with ANLS, that is, 107 realizations generate one fully negative column in a factor matrix which is then zeroed due to hard thresholding. This eventually leads to a decrease of the rank down to 4 and hence to a large reconstruction error (close to 10 -1 ). Among those 107 pathological cases, our simple function described by Alg.2 could cope with 92 of them without a significant increase in complexity. However, 15 realizations remain unsolved, because they correspond to either one of two particular cases: (i) either one column, say a(r o ), is fully negative, and the two others, namely b(r o ) and c(r o ) are fully positive, or (ii) all the three columns are fully negative. In order to cope with the latter cases, a straightforward improvement was brought in Alg. 3, by drawing a fresh column vector (also the absolute value of i.i.d drawn from a standard Gaussian distribution) to replace null vectors when generated in the unsolved pathological cases, before normalizing the columns. The results can be seen in Fig. 1 (bottom) and Fig. 2, where all 15 cases were solved and the rank was preserved.

Concluding remarks

In this paper, we have emphasized the fact that rank-1 tensors should not be treated as a collection of vectors without care, and showed an illustration in the case of ANLS using hard thresholding. In the latter case, two modifications have been proposed to fix the problem. In future works, we plan to investigate applications to other algorithms such as ADMM, and/or using soft thresholding. The influence of noise would also deserve to be further addressed.

  r) ⊗ b(r) ⊗ c(r) = ( ηa(r)) ⊗ ( b(r)) ⊗ (ηc(r)), ∀( , η) ∈ {-1, +1}. (4) This formula covers all 4 combinations of sign flipping of vectors, without affecting the result of the original outer product. Now denote by a (r) = η a(r), b (r) = b(r), and c (r) = η c(r), and: a -(r) = a (r) where a (r) < 0, and 0 elsewhere (5) b -(r) = b (r) where b (r) < 0, and 0 elsewhere (6) c -(r) = c (r) where c (r) < 0, and 0 elsewhere.

  The 4 possibilities created by sign flipping.

Algorithm 2

 2 Minimization of Φ 1 w.r.t. to and η for fixed r Require: a(r), b(r), c(r) initialize ( ), η( ) for each possible D [ ] ; (according to Table.1) normalize: a ← a(r) a(r) ; b ← b(r) b(r) ; c ← c(r) c(r) ; = 0; while ≤ 4, do = + 1; a ← ( ) η( ) a; b ← ( ) b; c ← η( ) c; compute a + , b + , and c +

  At the end the output of standard ANLS results in:At the end the output of Modified ANLS results in: 500 realizations of 10 × 5 matrices {A, B, C} are drawn. The rank of the tensor that is tested is hence R = 5. Entries of factor matrices are the absolute value of i.i.d. drawn from a standard Gaussian distribution. On each realization, both ANLS and a modified version based on the minimization of Φ 1 are run.

Fig. 1 .

 1 Fig. 1. Histograms of the error, in log10 scale, obtained after 500 iterations. Top: ANLS. Middle: ANLS modified with Φ1. Bottom: ANLS modified with Φ1 and with column reinitialization.

Fig. 2 .

 2 Fig.2. The average error obtained after 500 iterations, as a function of the number of iterations in ANLS (blue crosses), ANLS modified with Φ1 (red circles), and ANLS modified with Φ1 and with column reinitialization (yellow triangles).

This work was supported in part by ERC Advanced Grant 2013-320594 "Decoda".