Hassan Mortada

Vincent Mazet

Charles Soussen

Christophe Collet

Lionel Poisson

Parameterized Source Separation For Delayed Spectroscopic Signals

Keywords: Source separation, anechoic unmixing, correlated sources, delays

This paper addresses delayed (also known as anechoic) source separation when the source shape can be modeled by parameterized waveforms. An Alternating Least Squares (ALS) scheme is proposed to estimate the source shape parameters in a first step and both the mixing coefficients and the delays in a second step.

For the challenging delay parameter estimation step, we adopt a strategy inspired by greedy algorithms. For highly correlated sources, the separation becomes ambiguous, and a second algorithm is proposed: a regularization term is added to favor slow delay evolution within each source. Results on synthetic and real data demonstrate the effectiveness of both algorithms compared to state-of-the-art methods for highly correlated Gaussian waveforms.

Introduction

We consider the detection and estimation of patterns with varying characteristics in a temporal sequence of signals. This problem arises when dealing with decomposition of spectra, i.e., the estimation of emission lines (the so-called peaks) in a series of spectroscopic signals. The peak characteristics are assumed to evolve slowly through the sequence. In optical spectroscopy, acquisitions are obtained when varying a physical parameter. For instance, in time-resolved (b) The peaks and their labels. Each color represents the peaks sharing the same label.

photoelectron spectroscopy, several spectra are acquired at different times [START_REF] Stolow | Femtosecond time-resolved photoelectron spectroscopy[END_REF].

In astrophysics, the study of galaxy kinematics [START_REF] Gilmore | Kinematics, chemistry, and structure of the galaxy[END_REF][START_REF] Villeneuve | Nonlinear deconvolution of hyperspectral data with MCMC for studying the kinematics of galaxies[END_REF] leads to the acquisition of multispectral images where each pixel identifies with a spectrum whose peaks undergo varying delays due to the redshift and internal gas motions (Doppler effect). In chemometrics as well [START_REF] Malinowski | Factor analysis in chemistry[END_REF], external variations such as temperature or viscosity induce a spectrum shift and variation.

In many applications, prior knowledge or experiments allow the practitioner to model the peaks by a parameterized waveform, e.g., Gaussian [START_REF] Villeneuve | Nonlinear deconvolution of hyperspectral data with MCMC for studying the kinematics of galaxies[END_REF][START_REF] Mazet | Unsupervised joint decomposition of a spectroscopic signal sequence[END_REF] or Lorentzian [START_REF] Hollas | Modern spectroscopy[END_REF] functions. The first goal of this work is to estimate the amplitudes, delays and shape parameters of the peaks in each spectrum. The second goal is to match together the estimated peaks occurring in different spectra, that is to assign a distinct label to each peak arising in a given spectrum. We will assume that the shape parameters of the peaks having the same label are invariant from one spectrum to another. Fig. 1 shows an example of synthetic spectra and the expected output.

The problem was addressed in the Bayesian framework for a one-dimensional sequence of spectra [START_REF] Mazet | Unsupervised joint decomposition of a spectroscopic signal sequence[END_REF], modeled as the sum of Gaussian waveforms with Markovian priors on their parameters to promote a slow evolution. The parameters and labels were then estimated using a reversible jump Monte Carlo Markov chain algorithm. This method can deal with varying shape parameters of the peaks having the same label but suffers from high computation time.

Another statistical approach has been proposed in [START_REF] Vimond | Efficient estimation for a subclass of shape invariant models[END_REF], which models each signal in the sequence as a shifted version of a periodic, unknown waveform. The shift and the amplitude of the function are also unknown. In contrast to our problem, the author considers a unique waveform in each signal and there is no assumption of slow evolution of the unknowns through the sequence.

The problem can also be seen as a sparse approximation problem using an overcomplete dictionary, constructed by sampling the peak parameters in the case of parameterized peaks [START_REF] Austin | On the relation between sparse reconstruction and parameter estimation with model order selection[END_REF][START_REF] Bourguignon | Restoration of astrophysical spectra with sparsity constraints: Models and algorithms[END_REF][START_REF] Trigano | Cosa: An accelerated ISTA algorithm for dictionaries based on translated waveforms[END_REF][START_REF] Yaghoobi | Parametric dictionary design for sparse coding[END_REF]. Processing multiple spectra can be addressed by estimating the peak parameters separately for each spectrum, then matching the peaks together. However, this approach does not take the slow evolution knowledge into account during peak estimation. Alternatively, joint sparse approximation use a common dictionary for all the signals. Most approaches assume that the data signals are simultaneously sparse (they share a common support), which is often referred to as the Multiple Measurement Vector case [START_REF] Tropp | Algorithms for simultaneous sparse approximation. part I: Greedy pursuit[END_REF].

However, simultaneous sparsity is a very restrictive assumption which does not hold when the peaks have different delays in the spectra. Dynamic approaches like recursive sparse approximation [START_REF] Vaswani | Recursive recovery of sparse signal sequences from compressive measurements: A review[END_REF] compute the sparse approximation of a spectrum from the knowledge of the sparse approximation of the previous spectra. Social sparsity was introduced in [START_REF] Kowalski | Social sparsity! Neighborhood systems enrich structured shrinkage operators[END_REF] by promoting a structure between the sparse representation of consecutive data signals with less restrictive support assumptions. However, these approaches require a precise definition of the neighborhood which is not always possible, and hence they are not well-suited to the problem of peak estimation and labeling. Sparse approximation methods have also been used for deep learning [START_REF] Papyan | Theoretical foundations of deep learning via sparse representations: A multilayer sparse model and its connection to convolutional neural networks[END_REF] in case of a convolutional dictionary. However, such approaches require a large dataset of already decomposed data which is not always available.

Finally, the problem can be seen as a delayed source separation problem [START_REF] Comon | Handbook of Blind Source Separation: Independent Component Analysis and Applications[END_REF] where peaks with the same label are associated with unknown sources, and observed spectra are seen as mixtures. The reader is referred to [START_REF] O'grady | Survey of sparse and non-sparse methods in source separation[END_REF] for a review of source mixing models, e.g, instantaneous, anechoic and echoic, and the related source separation approaches. Delayed source separation is a specific case of convolutive source separation [START_REF] Castella | Quadratic higher order criteria for iterative blind separation of a MIMO convolutive mixture of sources[END_REF][START_REF] Kawamoto | Blind deconvolution of MIMO-FIR systems with colored inputs using second-order statistics[END_REF] where the filter kernels are Dirac functions. The main advantage of this approach is that estimation and labeling of peaks are done simultaneously, thus avoiding a two-step approach. Note that in optical spectroscopy, the number of sources is often smaller than the number of mixtures, yielding overdetermined source separation. Many contributions of the literature were addressed in audio signal processing where delayed source separation is also referred to as anechoic source separation. The sources are not parameterized and have to be estimated either in a blind or semi-blind framework. Since the problem is ill-posed, regularization is necessary, e.g., by imposing non-negativity of the sources and mixing coefficients (see [START_REF] Moussaoui | Separation of non-negative mixture of non-negative sources using a Bayesian approach and MCMC sampling[END_REF] in the case of instantaneous mixtures). Besides, strong assumptions such as source independence and decorrelation are often made. When dealing with delays, most approaches rely on delay linearization strategies such as applying Taylor expansion on the mixtures in the temporal domain [START_REF] Chabriel | An instantaneous formulation of mixtures for blind separation of propagating waves[END_REF][START_REF] Jiang | Covariance and time-scale methods for blind separation of delayed sources[END_REF]. Another strategy is to analyze the mixtures in the frequency [START_REF] Mørup | Shifted independent component analysis[END_REF][START_REF] Ouedraogo | A frequency method for blind separation of an anechoic mixture[END_REF] or time-frequency [START_REF] Nion | A time-frequency technique for blind separation and localization of pure delayed sources[END_REF][START_REF] Omlor | Anechoic blind source separation using Wigner marginals[END_REF] domains, so that delays become phases which allows one to use instantaneous source separation approaches coupled with phase estimation methods. However, the independence and non-correlation assumptions are often not valid in many realworld applications [START_REF] Duan | Restoration and separation of piecewise polynomial signals[END_REF]. Efforts have been done though to relax these assumptions.

In [START_REF] Puigt | Time-frequency ratio-based blind separation methods for attenuated and time-delayed sources[END_REF], the strong W-disjoint orthogonality assumption needed in the DUET algorithm [START_REF] Yilmaz | Blind separation of speech mixtures via time-frequency masking[END_REF] and the assumption that the sources must be orthogonal in the time-frequency domain are relaxed: each source must be dominant at least in one time-frequency window. However, this assumption becomes invalid when the evolution of a source delay is fast from one mixture to another or when two delayed source signals significantly overlap in some mixtures. In optical spectroscopy, some delayed source separation methods have been proposed.

In [START_REF] Harshman | Shifted factor analysis -Part I: Models and properties[END_REF][START_REF] Hong | Shifted factor analysis -Part II: Algorithms[END_REF], a time warping strategy was used to cancel the effect of the delays on each mixture after finding the delays over a predefined discrete grid using an exhaustive search strategy. The sources and amplitudes were then estimated in the linear least squares sense. However, the exhaustive search is not feasible when the delays take many values. In [START_REF] Mørup | Shifted non-negative matrix factorization[END_REF], the Non-negative Matrix Factorization (NMF) method was extended to consider delays: non-negative amplitudes and sources were found by using multiplicative updates, and delays were estimated with a gradient descent algorithm coupled with a maximizing cross-correlation procedure to reduce the effect of local minima. However, our experiments (see Section 5) showed that this method performs poorly for highly correlated sources.

In addition, it is not straightforward to impose a slow evolution of the delays.

Hereafter, the terms mixture and source respectively represent an observed spectrum and a peak having the same label in the observed spectra. An original algorithm is proposed to address delayed source separation when the sources are defined as parameterized functions. The problem is stated as an optimization problem in Section 2 and a sparse-based Alternating Least Squares (ALS) strategy is proposed in Section 3. This method was sketched in [START_REF] Mortada | Separation of delayed parameterized sources[END_REF] but in this paper we provide more detailed problem formulation and we analyze the separation limits of the method. The latter limits motivate us to develop the brand new method of Section 4: the slow delay evolution prior is considered to better discriminate spectrally overlapping and similar sources. Finally, results on synthetic and real data are presented in Section 5.

The notations are as follows. Bold and lowercase variables correspond to vectors. Bold and uppercase variables correspond to matrices. The i-th row and j-th column of a matrix M are respectively denoted as M i: and M :j .

The notation u + refers to u + = max(u, 0). The arrow notations and respectively refer to the decrement and increment operators.

Parameterized Source Separation

Let us consider I mixtures x i (λ), each being the noisy sum of J parameterized sources s(λ; w j ) such that:

x i (λ) = J j=1 a ij s(λ -c ij ; w j ) + n i (λ) i = 1, . . . , I (1)
where a ij ∈ R + is the amplitude of source j in mixture i (the application considered in this paper shows emission lines but not absorption lines), c ij ∈ R is the delay of source j in mixture i, w j ∈ R is the shape parameter of source j, and n i (λ) represents observation and modeling errors. Throughout the paper, the number of sources J is supposed to be known. Since the sources are parameterized, the delayed source separation problem comes to the estimation of the amplitudes a ij , delays c ij and shape parameters w j . For the sake of clarity, we will consider that the shape parameter is scalar, but the extension to a multidimensional parameter is straightforward. Furthermore, the sources are supposed to be modeled with the same parameterized function, however, the proposed method works even when the functions are different. Of course, using the same parameterized function yields highly correlated sources; we will see that the proposed method remains accurate in this case.

It is known [START_REF] Nion | A time-frequency technique for blind separation and localization of pure delayed sources[END_REF][START_REF] Rivet | Solving the indeterminations of blind source separation of convolutive speech mixtures[END_REF] that even in the noiseless case, model (1) suffers from at least three indeterminacies, namely, scale, permutation and phase inderminacies, the latter being related to the fact that a source signal is known up to some arbitrary delay. The use of parameterized sources helps to overcome the scale and the phase ambiguities by imposing their energy to be equal to a specified value and to set their maximum amplitude for λ = 0. Besides, the permutation ambiguity cannot be alleviated but does not yield any aftereffect.

The delays are supposed to be discretized over a grid with step ∆, thus

c ij = ij ∆ with ij ∈ N.
Note that ∆ may be lower than the sampling step of the mixture signals, which is set to 1 without loss of generality. In addition, the sources are supposed to be normalized. As said before, normalization allows us to overcome the scale indeterminacy. We denote by

x i = x i (1) x i (2) . . . x i (N )
T the vector gathering the samples of the i-th mixture. Similarly, the samples of each delayed source are gathered in a vector s[ ij ; w j ] defined as the unit-norm vector satisfying:

s[ ij ; w j ] ∝ s(1 -ij ∆; w j ) s(2 -ij ∆; w j ) . . . s(N -ij ∆; w j ) T , (2)
where ∝ refers to proportionality. Equation (1) now reads:

x i = J j=1 a ij s[ ij ; w j ] + n i , i = 1, . . . , I, (3)
where n i refers to the noise in mixture i. Supposing the noise to be white and Gaussian, the maximum likelihood estimator is obtained by minimizing the criterion:

E(A, L, w) = I i=1 ε(A i: , L i: , w), (4)
where ε(A i: , L i: , w) is the quadratic error related to mixture i:

ε(A i: , L i: , w) = x i - J j=1 a ij s[ ij ; w j ] 2 2 , (5)
and

A ∈ R I×J + , L ∈ N I×J , w ∈ R J
+ respectively gather the amplitudes a ij , delays ij and shape parameters w j for mixtures i and sources j:

A =      a 11 . . . a 1J . . . . . . . . . a I1 . . . a IJ      =      A 1:
. . .

A I:      (6)
L =      11 . . . 1J . . . . . . . . . I1 . . . IJ      =      L 1:
. . .

L I:      (7) w = w 1 . . . w J T . (8)
Therefore, the source separation problem is formulated as the following constrained minimization problem:

min A≥0,L,w E(A, L, w). (9)

Sparse-based Alternating Least Squares

The optimization problem ( 9) is challenging because of the non-convexity of the criterion E, induced by the nonlinearity of model ( 3) with respect to L and w. As an example, Fig. 2 displays the variations of criterion E(A, L, w) with respect to L in the case of I = 1 mixture and J = 2 sources: one can see that it admits multiple local minimizers as well as flat regions, making its optimization difficult even for this simple example.

ALS is an iterative descent strategy in which criterion E is minimized with respect to a block of variables while fixing the others, and vice versa. The algorithm stops if the criterion decrease at one iteration becomes lower than a tolerance ρ. ALS is not guaranteed to converge towards the global minimizer of ( 9) since it is a block minimization of a non-convex criterion. However, this method is often used in delayed source separation, where e.g., the sources, delays and amplitudes are alternately estimated [START_REF] Hong | Shifted factor analysis -Part II: Algorithms[END_REF][START_REF] Mørup | Shifted independent component analysis[END_REF][START_REF] Nion | A time-frequency technique for blind separation and localization of pure delayed sources[END_REF]. However, the latter methods suffer from several limitations, as stated in section 1, making them unsuitable for the considered problem.

The proposed ALS scheme is given in Algorithm 1. The criterion is alternately Algorithm 1: ALS scheme for min A,L,w E(A, L, w).

Initialization: A = L = 0 I×J , w∼ U [wmin,wmax] J 1 do 2 (A 0 , L 0 , w 0 ) ← ( A, L, w) 3 for i = 1 I do 4 ( A i: , L i: ) ← argmin Ai:,Li: ε(A i: , L i: , w) (see Algorithm 2) 5 end 6 w ← argmin w E A, L, w 7 while E(A 0 ,L 0 ,w 0 )-E( A, L, w) E(A 0 ,L 0 ,w 0 ) ≥ ρ
optimized with respect to the source shape parameters w on the one hand and the delays L and amplitudes A on the other hand. The shape estimation subproblem is a continuous non-linear least-squares problem:

w ← argmin w E(A, L, w), (10)
and is solved using the Levenberg-Marquardt algorithm [START_REF] Marquardt | An algorithm for least-squares estimation of nonlinear parameters[END_REF]. In the first iteration, the shape parameters are initialized from the uniform distribution U [wmin,wmax] J

where w min and w max are the extreme values defined by the user. In the following iterations, w is initialized with the estimates obtained at the previous iteration.

The amplitude and delay estimation is the main challenge:

( A, L) ← argmin A≥0,L E(A, L, w). (11)
It is detailed hereafter.

Amplitude and Delay Estimation

It follows from (4) that problem ( 11) is separable to I independent subproblems: min Ai:≥0,Li: We propose to minimize each ε(A i: , L i: , w) using a sparse approximation algorithm aiming to sparsely represent mixture i in an overcomplete dictionary. This choice is justified because spectroscopic signals contain very few peaks.

ε(A i: , L i: , w) ∀i. ( 12
)
x i ≈ S 1 S 2 S 3 α i3 α i2 α i1

Dictionary Formulation

The dictionary is a block matrix [S 1 , . . . , S J ], with S j ∈ R N ×M the block gathering M delayed versions of source j (see Fig. 3):

S j = s[0; w j ] s[1; w j ] ... s[M -1; w j ] . (13)
Each mixture x i is approximated as:

x i ≈ J j=1 S j α ij ( 14
)
where α ij = [0, . . . , a ij , . . . , 0] T ∈ R M + is a 1-sparse vector, so that each source appears at most once in each mixture (see Fig. 3). The value and index of the non-zero element in α ij respectively indicate the amplitude a ij and delay ij . Thus, the optimization problem ( 12) can be rewritten as:

min ∀j,αij ≥0 x i - J j=1 S j α ij 2 2 s.t. ∀ j, α ij 0 ≤ 1, ( 15
)
where the 0 "norm" • 0 counts the number of non-zero coefficients in a vector.

NN-OMP-like Implementation for Delayed Source Separation

Greedy algorithms are effective and efficient when the sparsity level J is small and known [START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF]. In addition, their structure is simple and can be easily adapted to the recovery of structured sparse representations. The non-negative orthogonal matching pursuit (NN-OMP) algorithm [START_REF] Bruckstein | On the uniqueness of nonnegative sparse solutions to underdetermined systems of equations[END_REF] is an iterative algorithm composed of three steps:

1. the so-called forward selection step consists in choosing the column of the dictionary that is the most positively correlated with the residual;

2. the amplitudes corresponding to the chosen columns are updated by solving a non-negative least-squares estimation problem [START_REF] Lawson | Solving least squares problems[END_REF];

3. the residual is updated by removing the contributions of the chosen columns.

NN-OMP considers the constraint in (15) to be J j=1 α ij 0 ≤ J. However, this constraint does not enforce the sources to appear at most once in each mixture.

Therefore, the proposed implementation (Algorithm 2) consists in forcing the sparse vector to be structured in blocks, each α ij being 1-sparse:

1. the dictionary column that is the most positively correlated with the residual is selected, yielding the corresponding source j and its delay i j (lines 2-6). The selected source j is then added to the list J of selected sources (line 7) so that the dictionary columns embedded in the blocks S j for sources j ∈ J will not be tested in the next iterations; 2. the amplitudes of the selected sources gathered in J are estimated by solving the non-negative linear least-squares problem:

A i: ← argmin Ai: ε(A i: , L i: , w) s.t.    A iJ ≥ 0 A iJ = 0, (16)
where J denotes the complementary subset of J (line 8); Algorithm 2: Implementation of ( A i: , L i: ) ← argmin Ai:,Li: ε(A i: , L i: , w)

Initialization: A i: = L i: = 0 1×J , J = ∅, r i = x i 1 for k = 1 J do 2 for j ∈ {1, . . . , J}\J do 3 ˜ ij ← argmax r T i s[ ; w j ] + 4 end 5 j ← argmax j / ∈J r T i s[ ˜ ij ; w j ] + 6 i j ← ˜ i j 7 J ← J ∪ { j } 8
Update amplitudes A iJ according to ( 16)

9 r i ← x i - j∈J a ij s[ ij ; w j ] 10 end 11 if ε(A 0 i: , L 0 i: , w) < ε( A i: , L i: , w) then ( A i: , L i: ) ← (A 0 i: , L 0 i: ) end
3. lastly, the residual vector r i is updated (line 9).

The main difference between NN-OMP and our implementation lies in the first step. Also, line 11 ensures a decrease of the criterion by invalidating the estimates if they produce a criterion value greater than the value obtained at the previous iteration of Algorithm 1. Furthermore, Algorithm 2 offers the possibility to obtain a variable number of sources per mixture: an additional stopping criterion may be added such that the loop breaks if the residual norm r i 2 becomes lower than a threshold, e.g., related to the noise variance.

Finally, we also investigated recent methods to refine delay estimation, which is constrained to be on a grid. Indeed, the discretization induced by sampling the sources in the dictionary yields a bias in the estimation, as already been discussed

in [START_REF] Trigano | Sparse reconstruction algorithm for nonhomogeneous counting rate estimation[END_REF]. The easiest way to reduce this bias is to decrease the value of ∆, but this results in a bigger dictionary, and in turn, in an increase of the computational burden. So, we tested an interpolation extension proposed in [START_REF] Ekanadham | Recovery of sparse translation-invariant signals with continuous basis pursuit[END_REF][START_REF] Fyhn | Compressive parameter estimation for sparse translation-invariant signals using polar interpolation[END_REF]. It turned out that this extension is more time-consuming than a simple decrease in the sampling step ∆. Therefore, we rather choose to work with a fine grid.

Shape Discriminating Limit

Algorithm 1 is able to assign the estimated peaks to the right source because the sources can be discriminated by their shape parameter (see Fig. 3). We are interested in finding the resolution limit, that is the least difference between the shape parameters of two sources beyond which the sources can be discriminated.

For this purpose, we consider I = 40 mixtures each with N = 200 samples, and J = 2 Gaussians sources of widths w 1 and w 2 and with constant delays through the mixtures. We gradually vary the ratio w 2 /w 1 from 0.5 to 1.5. For each ratio, we measure the switch percentage defined as the percentage of wrongly assigned peaks over the total number of peaks (a peak is wrongly assigned if it belongs to source 1 while it is assigned to source 2 and vice versa). The experiment is repeated for three SNRs (Fig. 4). The SNR is defined as ten times the log-ratio of the mean energy of the noiseless mixtures with the noise variance σ 2 n :

SNR = 10 log 10 1 σ 2 n • I • N I i=1 J j=1 a * ij s[ * ij ; w * j ] 2 2 , ( 17
)
where the * symbol represents the value of the ground-truth parameters. The results show that for high SNR, the proposed method can separate the sources whatever the ratio w 2 /w 1 (except of course when both parameters are equal).

The ability to separate the sources decreases with the SNR. This motivates us to introduce a regularization term to overcome this limitation.

Slow Delay Evolution Enforcement

Regularized Criterion

This section aims at promoting slow evolution for each source delay L through the mixtures for two reasons. First, the slow evolution is, in practice, the consequence of a short acquisition time between measurements (such as in photoelectron spectroscopy [START_REF] Gloaguen | Experimental evidence for ultrafast electronic relaxation in molecules, mediated by diffuse states[END_REF]) or neighboring sensors (such as in galaxy kinematics [START_REF] Gilmore | Kinematics, chemistry, and structure of the galaxy[END_REF] or audio recorded mixtures [START_REF] Chabriel | An instantaneous formulation of mixtures for blind separation of propagating waves[END_REF]). Second, it would help to discriminate highly correlated sources. To this end, a regularization term R(L) is added to the data-fit term E(A, L, w) defined in (4):

F (A, L, w) = E(A, L, w) + τ R(L) ( 18
)
where τ is the regularization parameter to be set by the user. For the sake of clarity, we assume in the sequel that i is a one-dimensional index. The regularization term R(L) measures the sum of squared differences between consecutive delays:

R(L) = I i=2 J j=1 ij -(i-1)j 2 . ( 19
)
Criterion F is optimized using Algorithm 3 which is based on an ALS scheme.

Since R(L) does not depend on w, the shape estimation (line 4) is identical to [START_REF] Fyhn | Compressive parameter estimation for sparse translation-invariant signals using polar interpolation[END_REF]. Thus, it is computed using the Levenberg-Marquardt algorithm. On the contrary, the amplitude and delay estimation step (line 3) differs from ( 11)- [START_REF] Gloaguen | Experimental evidence for ultrafast electronic relaxation in molecules, mediated by diffuse states[END_REF] because F does not read as a separable sum with respect to L i: . Indeed, the Algorithm 3: ALS scheme for min A,L,w F (A, L, w).

Initialization: A = L = 0 I×J , w∼ U [wmin,wmax] J 1 do 2 (A 0 , L 0 , w 0 ) ← ( A, L, w) 3 ( A, L) ← argmin A,L F (A, L, w) (see Algorithm 4) 4 w ← argmin w F A, L, w 5 while F (A 0 ,L 0 ,w 0 )-F ( A, L, w) F (A 0 ,L 0 ,w 0 ) ≥ ρ terms ( ij -(i-1)j ) 2
appearing in the regularization term R(L) not only depend on ij but also on (i-1)j . Let us now detail the update rules for A and L.

Algorithm 4: Implementation of ( A, L) ← argmin A,L F (A, L, w)

Initialization: A = L = L temp = 0 I×J , J = ∅, r i = x i ∀i 1 for k = 1 J do 2
for j ∈ {1, . . . , J}\J do 3

Compute L :j defined in (22) using Algorithm 5 At each iteration, the source j inducing the largest decrease of criterion F is selected (lines 2-5) and added to the list of selected sources J (line 6). The corresponding delays L temp : j are computed. Then, the amplitudes of the selected sources in J in all the mixtures are estimated using a non-negative least squares solver (line 8). Finally, the residual vectors of all the mixtures are updated (line 11). Note that adding a new source results in a decrease of the data-fit term E(A, L, w) but makes the regularization term R(L) increasing. In other words, the criterion F (A, L, w) = E(A, L, w) + τ R(L) can either increase or decrease.

4 end 5 j ← argmax j / ∈J I i=1 r T i s[ ij ; w j ] 2 + -τ I i=2 ij -(i-1)j 2 6 J ← J ∪ { j } 7 L temp : j ← L : j 8 A temp ← argmin A F (A, L temp , w) s.t. {A :J ≥ 0, A :J = 0} 9 if F (A temp , L temp , w) > F ( A, L, w) then Break end A, L ← A temp , L temp for i = 1 I do r i ← x i - j∈J a ij s[ ij ; w j ] end end if F (A 0 , L 0 , w) < F ( A, L, w) then ( A, L) ← (A 0 , L 0 ) end

Amplitude and Slow Delay Estimation

Therefore, the condition in line 9 is set to break the loop in case of an increase of F , so that Algorithm 4 is indeed a descent algorithm. Similarly, the condition in line 13 ensures a decrease of the criterion by invalidating the estimates if they produce a criterion value that is larger than the value obtained at the previous iteration of Algorithm 3.

Delay Estimation with an ICM-like Algorithm

Let us now specify the rule for selecting the source j among all candidate sources j / ∈ J . Since j is defined as the source yielding the largest decrease of criterion F , lines 2-4 in Algorithm 4 aim at estimating, for each source j that has not already been selected, the value of the corresponding delays L :j (denoted by L:j ) as well as the value of F obtained while considering the set of sources J ∪ {j}. To do so, one needs to consider the minimization of F (A, L, w) with respect to A :j and L :j , while fixing the values of L :j and A :j for j ∈ J :

L :j ← argmin L:j min A:j ≥0 I i=1 r i -a ij s[ ij ; w j ] 2 2 + τ I i=2 ( ij -(i-1)j ) 2 . ( 20
)
The minimization of (20) with respect to A :j while fixing L :j has a closed form solution:

∀i, a ij = r T i s[ ij ; w j ] + . (21)
Plugging back ( 21) into ( 20), [START_REF] Malinowski | Factor analysis in chemistry[END_REF] simplifies to:

L :j ← argmax L:j I i=1 r T i s[ ij ; w j ] 2 + -τ I i=2 ( ij -(i-1)j ) 2 . ( 22
)
which is the cost function appearing at line 5 in Algorithm 4.

Because ( 22) is a combinatorial problem, we resort to the Iterated Conditional Modes (ICM) algorithm [START_REF] Besag | On the statistical analysis of dirty pictures[END_REF] which is a popular coordinate-wise optimization method in image processing. While it converges to a local optimizer, it generally gives good results. The proposed implementation to solve ( 22) is given in Algorithm 5. Specifically, at each iteration (lines 9-15), all mixtures i ∈ {1, . . . , I} are swept using the following way: a starting mixture i 0 is randomly chosen (line 10), then the delays in mixtures i 0 + 1 to I are sequentially estimated, and the same procedure is used for mixtures i 0 -1 to 1. Whenever a mixture i ∈ {1, . . . , I} is visited, the related delay ij is estimated (in lines 12-14) by maximizing the criterion in [START_REF] Masson | Coupled electronic and structural relaxation pathways in the postexcitation dynamics of Rydberg states of BaAr N Clusters[END_REF] with respect to ij whilst fixing the other delays i j with i = i. The convergence is reached when the relative distance between the estimates of two consecutive iterates is smaller than a tolerance threshold ξ (line 16).

The initialization is performed in the same way (lines 1-7), the major difference being that, at the first iteration, some delays are not yet estimated and, in consequence, cannot be considered in [START_REF] Masson | Coupled electronic and structural relaxation pathways in the postexcitation dynamics of Rydberg states of BaAr N Clusters[END_REF]. So, we simply discard the terms with unknown values from the equation. Note also that the first mixture to be considered is the one maximizing the data-fit term (line 1).

Remarks

If no regularization is considered (i.e., τ = 0), then F (A, L, w) = E(A, L, w).

We then recommend to use Algorithm 1 rather than Algorithm 3 since the former exploits the separability of the criterion. Conversely, when τ tends to infinity, the Algorithm 5: Implementation of

L :j ← argmax L:j I i=1 r T i s[ ij ; w j ] 2 + -τ I i=2 ( ij -(i-1)j ) 2
Initialization:

L :j = 0 I×1 1 (i 0 , i0j ) ← argmax i, r T i s[ ; w j ] 2 + 2 for i = (i 0 + 1) I do 3 ij ← argmax r T i s[ ; w j ] 2 + -τ ( -(i-1)j ) 2
4 end

5 for i = (i 0 -1) 1 do 6 ij ← argmax r T i s[ ; w j ] 2 + -τ ( -(i+1)j ) 2 7 end 8 do 9 L old :j ← L :j 10 i 0 ∼ U {1,...,I} 11
for i = i 0 I and i = (i 0 -1) 1 do

12 if i = 1 then ij ← argmax r T i s[ ; w j ] 2 + -τ ( -(i+1)j ) 2 end 13 if i = I then ij ← argmax r T i s[ ; w j ] 2 + -τ ( -(i-1)j ) 2 end 14 if 1 < i < I then ij ← argmax r T i s[ ; w j ] 2 + -τ ( -(i-1)j ) 2 -τ ( -(i+1)j ) 2 end 15 end 16 while L :j -L old :j 2 2 / L old :j 2 2 < ξ
delays related to each source are necessarily constant: ij = i j , ∀i, i . Therefore, the sparse vectors α ij ∀i share a common support (see Fig. 5(c)). The estimation of L :j in ( 22) becomes:

L :j = [ j , . . . , j ] with j ← argmax I i=1 r T i s[ ; w j ] 2 + , (23)
and can be obtained using an exhaustive search by testing the M possibilities for . This is very similar to the S-OMP algorithm [START_REF] Tropp | Algorithms for simultaneous sparse approximation. part I: Greedy pursuit[END_REF] (which is a greedy algorithm

α i1 α i2 α i3 (a) i α i1 α i2 α i3 (b) i α i1 α i2 α i3 (c)
i for sparse recovery of vectors having a common support) with the difference that the vectors α ij are 1-sparse.

Numerical Results

Comparison on Synthetic Mixtures

Algorithms 1 and 3 are compared with state-of-the-art methods that are able to separate highly correlated sources.

Evaluation of Non-Parameterized Methods

First of all, we evaluate two competing methods dealing with delayed source separation [START_REF] Mørup | Shifted non-negative matrix factorization[END_REF] and convolutive source separation [START_REF] Castella | Quadratic higher order criteria for iterative blind separation of a MIMO convolutive mixture of sources[END_REF]. The latter estimates the i.e., s(λ; w j ) = exp(-λ 2 /2w 2 j ). The amplitudes and delays of each Gaussian waveform are continuously generated from polynomials of degree 2, 3 or 4.

Besides, the shape parameters, i.e., the standard deviations of the Gaussian sources are set to w * = 4 4 1.5 6 T . The data are displayed in Fig. 6 and show noticeable behavior: two sources (around λ = 50) are very close and highly correlated (Gaussians with same width, equal to 4); there are several overlaps of the sources; and the number of sources per mixture is not constant.

The estimated sources with the non-parameterized methods [START_REF] Castella | Quadratic higher order criteria for iterative blind separation of a MIMO convolutive mixture of sources[END_REF][START_REF] Mørup | Shifted non-negative matrix factorization[END_REF] are shown in Fig. 7, and the mixture reconstruction and estimated amplitudes and delays are displayed in Fig. 8. The mixture reconstruction obtained with the method of [START_REF] Mørup | Shifted non-negative matrix factorization[END_REF] is half satisfactory. On the one hand, the sources s 1 and s 2 are Gaussian-shaped and their amplitudes and delays roughly correspond to the ground truth. On the Figure 7: Ground truth and estimated sources with the methods of [START_REF] Mørup | Shifted non-negative matrix factorization[END_REF] and [START_REF] Castella | Quadratic higher order criteria for iterative blind separation of a MIMO convolutive mixture of sources[END_REF].

other hand, the sources s 3 and s 4 are bimodal (so they do not match with the ground truth) and actually fit the two ground-truth sources with similar shapes around λ = 50. Regarding the method of [START_REF] Castella | Quadratic higher order criteria for iterative blind separation of a MIMO convolutive mixture of sources[END_REF], the estimated sources present a dominant peak with small oscillations, so they are not close to the ground truth.

Besides, the mixture reconstruction is very noisy and the parameter estimates do not completely match with the ground-truth. In conclusion, this example shows that the methods of [START_REF] Castella | Quadratic higher order criteria for iterative blind separation of a MIMO convolutive mixture of sources[END_REF][START_REF] Mørup | Shifted non-negative matrix factorization[END_REF] are not able to deal with highly correlated sources, partly because they do not consider the highly informative knowledge of parametric sources.

Comparison with a Parameterized Method

We now compare the non-regularized and regularized methods (which respectively refer to Algorithms 1 and 3) with a Bayesian method [START_REF] Mazet | Unsupervised joint decomposition of a spectroscopic signal sequence[END_REF] that explicitely exploits the knowledge of parametrized source shapes for decomposing a sequence of spectra. The methods are tested with the data displayed in Fig. 6. The three methods model the sources by Gaussians so the shape parameter is the Gaussian width. This is a positive scalar, therefore a trust-region reflective algorithm is used instead of the Levenberg-Marquardt algorithm because it is able to consider the positivity constraint. The stopping constant of the proposed methods is ρ = 10 -4 ; the delay sampling step is set to ∆ = 0.2; the regularization parameter is empirically tuned to τ = 1.5 • 10 -2 . The ICM stopping constant is set to Reconstruction Delays Amplitudes

Method of [START_REF] Mørup | Shifted non-negative matrix factorization[END_REF] The results are displayed in Fig. 9. The reconstruction of the mixtures is equally good for the proposed methods and the method of [START_REF] Mazet | Unsupervised joint decomposition of a spectroscopic signal sequence[END_REF]. This is confirmed by Tab. 1 which compares the methods in terms of (i) computation time; (ii)

Mean Squared Error (MSE) defined as

I i r i 2
2 /(N • I); (iii) amplitude, delay and shape error, defined as:

C * -C 2 F C * 2 F , A * -A 2 F A * 2 F and w * -w 2 2 w * 2 2 ( 24
)
where • F denotes the Frobenius norm and C = ∆L is the delay matrix for the proposed methods.

The estimated shape parameters are very close to the ground-truth, respectively w = 4.0 3.9 1.4 5.6

T and w = 3.9 4.0 1.4 5.9

T for the non-regularized and regularized methods. With the method of [START_REF] Mazet | Unsupervised joint decomposition of a spectroscopic signal sequence[END_REF], the shape parameters can vary within a source, and a source does not necessarily appear in all the mixtures. But the variations remain small for the considered data so that we can consider the means of the estimated shape parameters: w = 3.9 3.9 1.5 6.0 T which is also close to the ground-truth.

Reconstruction Delays Amplitudes

Non-regularized method Considering the non-regularized method, the delay and amplitude estimates are not satisfactory despite the good reconstruction and source estimation. This can be explained by the fact that when some sources have the same shape parameter, the identification becomes ambiguous (see Section 3.2). On the contrary, the regularized method yields improved estimates as well as the method of [START_REF] Mazet | Unsupervised joint decomposition of a spectroscopic signal sequence[END_REF]. These results are numerically validated in Tab. 1.

Conclusion

These simulations show that a parametrized method is of interest when dealing with highly correlated sources. The proposed regularization ensures to non-regularized regularized Table 1: Numerical performance of the compared methods on the synthetic mixtures plotted in Figure 6. The delay and amplitude error are defined in [START_REF] Mazet | Unsupervised joint decomposition of a spectroscopic signal sequence[END_REF].

separate correctly the sources with very similar shape parameters. The proposed methods give accurate results within a very low computation time. On the contrary, the method of [START_REF] Mazet | Unsupervised joint decomposition of a spectroscopic signal sequence[END_REF], while very effective, is 15 to 30 times slower than the proposed methods. Note however, that the latter model is more versatile since the shape parameter are allowed to vary within a source and additional priors on the shapes and amplitudes are considered.

Influence of the SNR

We now investigate the influence of the SNR on the performance of the proposed methods. A first example of synthetic mixtures with I = 15 and J = 2 is displayed in Fig. 10 to show that the proposed methods are robust to the SNR: the conclusion of section 5.1 remains, that is, both proposed methods yield good reconstruction and source quality whatever the SNR. Also, Fig. 10 shows the benefit of using a regularization on the delays. Recall that the permutation indetermination, which can be clearly noticed in this example, has no influence on the parameter estimation.

Furthermore, we perform statistical simulations to compare both proposed methods with the method of [START_REF] Mazet | Unsupervised joint decomposition of a spectroscopic signal sequence[END_REF]. The methods of [START_REF] Castella | Quadratic higher order criteria for iterative blind separation of a MIMO convolutive mixture of sources[END_REF] and [START_REF] Mørup | Shifted non-negative matrix factorization[END_REF] are not tested amplitudes are generated by polynomial functions with random coefficients and degrees varying between 2 and 4; the delay sampling step is equal to ∆ = 0.2; the stopping constant is set to ρ = 10 -4 for both proposed methods; the ICM stopping constant is set to ξ = 10 -4 and the regularization parameter is set to τ = 4 • 10 -3 . The results are shown in Fig. 11 where, for each SNR, the average of the results obtained for the 100 generated datasets is plotted.

As expected, the MSE decreases as the SNR increases for the three methods (Fig. 11a), and it should be noticed that the variations are very similar. The delay and amplitude errors of the non-regularized method improve when the SNR increases but it is shown that the non-regularized method has lower performance than its competitors. On the contrary, the regularized method gives improved estimation of the amplitudes and delays. In addition, introducing a regularization on the delays yields estimates that are robust to the noise variation; this is not the case for the amplitudes for which no assumption of slow evolution is considered.

The method of [START_REF] Mazet | Unsupervised joint decomposition of a spectroscopic signal sequence[END_REF] behaves similar to the regularized method in terms of delay error, whilst its amplitude error is generally worse than the proposed methods (except for low SNR). Again, it is shown that the method of [START_REF] Mazet | Unsupervised joint decomposition of a spectroscopic signal sequence[END_REF] is much slower than the proposed methods (approximately 20 times slower). Therefore, the regularized method is competitive with the method in [START_REF] Mazet | Unsupervised joint decomposition of a spectroscopic signal sequence[END_REF] with much lower computation time.

Results on Real Photoelectron Data

Time-resolved photoelectron spectroscopy [START_REF] Stolow | Femtosecond time-resolved photoelectron spectroscopy[END_REF] is an experimental tool which allows to study the energy relaxation occurring after absorption of a photon by an isolated molecule, atom or a blend of both. The energy relaxation is probed by ionization of the excited system with another delayed photon, thereby ejecting a so-called photoelectron. The distribution of the photoelectrons according to their energy is measured at different times to get a temporal sequence of photoelectron spectra and is given by the estimated source characteristics (amplitudes, delays, shapes). The delays C indicate the energy from which photoemission is occurring and the area under each source, which depends on both the amplitudes A and shapes w, corresponds to the relative number of emitted photoelectrons at the corresponding energy. The experiment presented in this paper studies the relaxation of an atom of barium [START_REF] Masson | Dynamics of highly excited barium atoms deposited on large argon clusters. I. General trends[END_REF][START_REF] Mazet | Unsupervised joint decomposition of a spectroscopic signal sequence[END_REF]. The regularized method was applied to the sequence of photoelectron spectra shown in Fig. 12: the goal is to determine how the energy, intensity, and width of the peaks (i.e., delays, amplitudes and shapes of the sources) evolve through the sequence, indicating the temporal changes undergone by the studied system. Figure 13: The reconstruction and the estimated continuum, delays and amplitudes using the regularized method (first column) and the method of [START_REF] Mazet | Unsupervised joint decomposition of a spectroscopic signal sequence[END_REF] (second column). The estimated shape parameters of regularized method are w = [0.057, 0.049, 0.025, 0.09, 0.06] T and for the method of [START_REF] Mazet | Unsupervised joint decomposition of a spectroscopic signal sequence[END_REF] after averaging the estimates for each source are w = [0.066, 0.064, 0.039, 0.037, 0.071, 0.05] T .

functions, which is the usual model for the peaks in photoelectron spectroscopy1 .

The grid sampling step is set to ∆ = 5 • 10 -4 and the regularization term is set to τ = 2 • 10 -5 . The stopping constant ρ and the ICM stopping constant are respectively set to ρ = 10 -4 and ξ = 10 -4 . Besides, the photoelectron spectra include a background that must be estimated and removed from the spectra.

The background in each mixture i is modeled as an exponential of the form α i exp(-λ/β). An additional step is added to the ALS scheme in Algorithm 3

to estimate the new unknowns (the I weights α i and a single β value) by using the Levenberg-Marquardt algorithm.

The results obtained with the regularized method and the method of [START_REF] Mazet | Unsupervised joint decomposition of a spectroscopic signal sequence[END_REF] are presented in Fig. 13: both approaches appear roughly similar concerning delay and shape error. However, differences can be observed: (i) the estimated continuum provided is much smoother for the regularized method, which was a source of noise in the method [START_REF] Mazet | Unsupervised joint decomposition of a spectroscopic signal sequence[END_REF]; (ii) the regularized method rather improves the shape of the main band by addition of nearby new sources than focusing on the low intensity energy ones as seen for method [START_REF] Mazet | Unsupervised joint decomposition of a spectroscopic signal sequence[END_REF]. Such behavior looks reasonable according to the shape of the decay simulated [START_REF] Masson | Coupled electronic and structural relaxation pathways in the postexcitation dynamics of Rydberg states of BaAr N Clusters[END_REF]. (iii) The method of [START_REF] Mazet | Unsupervised joint decomposition of a spectroscopic signal sequence[END_REF] is able to deal with a varying number of peaks through the mixtures, whereas in the proposed method the J sources are present in all the mixtures.

Nevertheless, small amplitudes can be cancelled by applying a threshold. (iv)

Last but not least, the proposed method is about 560 times faster (7.5 seconds for the regularized method versus 70 minutes for the method of [START_REF] Mazet | Unsupervised joint decomposition of a spectroscopic signal sequence[END_REF]). In conclusion, the proposed method appears to be both effective and efficient.

Software

The Matlab codes associated with this article will be freely available in case of publication.

Figure 1 :

 1 Figure 1: (a) Synthetic sequence of spectra, each containing 3 Gaussian waveforms (peaks).

12 ∆Figure 2 :

 122 Figure 2: (a) A noisy mixture of two Gaussian sources (I = 1, J = 2) with a * 11 = 6.0, a * 12 = 2.7, * 11 ∆ = 50, * 12 ∆ = 65, w * 1 = 5, w * 2 = 4, N = 100. (b) The criterion L 1: → ε(A * 1: , L 1: , w * ) admits local minimizers and flat surfaces. The red × indicates the global minimizer.

Figure 3 :

 3 Figure 3: Sparse representation model of a mixture x i with J = 3 sources. Each block S j of the dictionary gathers the delayed versions of s[ ij ; w j ]. The sparse representation is structured so that each block vector α ij is 1-sparse.

Figure 4 :

 4 Figure 4: Switch percentage with respect to the ratio of Gaussian widths, for three SNRs.

Algorithm 4

 4 is proposed to optimize F (A, L, w) with respect to A and L. It is a greedy algorithm in the spirit of Algorithm 2, and, in this way, reproduces the three steps of the NN-OMP framework already stated in Section 3.1.2. However, Algorithm 4 takes all mixtures as inputs while Algorithm 2 considers a single mixture. The estimation of the source delays in line 3 of Algorithm 4 is deferred to Algorithm 5.

Figure 5 :

 5 Figure 5: Examples of sparse model for J = 3 sources and I = 6 mixtures where the white and the colored squares respectively represent the zero and non-zero amplitudes. Each column corresponds to a mixture i, and is divided into J blocks representing the sources (each source is indicated with a unique color). (a) Case of no delay regularization. (b) Case of a slow-moving regularization. (c) Case of a very strong regularization, resulting in a simultaneous sparse approximation with constant support for all mixtures.

Figure 6 :

 6 Figure 6: Synthetic data used for comparison: 40 mixtures with 4 sources and SNR = 15 dB.

  Castella et. al[START_REF] Castella | Quadratic higher order criteria for iterative blind separation of a MIMO convolutive mixture of sources[END_REF]

Figure 8 :

 8 Figure 8: Results on the synthetic data of Fig. 6. The columns respectively the reconstructed mixtures, the estimated delays and the estimated amplitudes. Each source is represented with unique color and marker. The ground truth delays and amplitudes are plotted in gray lines.

Figure 9 :

 9 Figure 9: Results on the synthetic data of Fig. 6. The first column displays the reconstructed mixtures as well as the estimated sources. The second and third columns display the estimated delays and amplitudes. Each source is represented with unique color and marker. The ground truth delays and amplitudes are plotted in gray lines.

Figure 10 :

 10 Figure 10: Synthetic data with 2 sources at two different SNR, and the related reconstruction and estimated sources with both proposed methods.

  Delay error, defined in[START_REF] Mazet | Unsupervised joint decomposition of a spectroscopic signal sequence[END_REF].

Figure 11 :

 11 Figure11: Comparison of the non-regularized (green ) and the regularized (red •) proposed methods with the method of[START_REF] Mazet | Unsupervised joint decomposition of a spectroscopic signal sequence[END_REF] (blue +) with respect to the SNR.

Figure 12 :

 12 Figure 12: Photoelectron spectra (I = 44 mixtures).

  The sequence gathers I = 44 spectra (covering a duration of 3.47 ps), each of N = 181 samples (from 0.02 eV to 2.52 eV). The sources are modeled by

Note that choosing a model not "consistent" with the shape of real peaks may yield an increase of the number of estimated sources.

Acknowledgment

This work was supported by the DSIM project (grant ANR-14-CE27-0005).

Conclusion

This paper studies the delayed source separation when source signals are supposed to be parameterized, which is a valid assumption in many spectroscopy applications. An ALS scheme is proposed to alternatively estimate the shape parameters and the amplitudes and delays. For the challenging delay parameter estimation step, we propose a greedy strategy using parametric dictionaries whereas the sparse representation is structured to respect the mixing model. We present two algorithms that follow the same ALS scheme: in the first, only the data-fit criterion is considered whilst in the second an additional regularization term is added to promote slow delay evolution within each source. The proposed methods outperform the state-of-the-art delayed source separation methods when sources are highly correlated. Furthermore, it is as effective as the best competitors with much lower computation time. Results on real data confirm the effectiveness and the efficiency of the proposed methods. Future works will be dedicated to studying the benefit of a warm-start initialization of A and L in Algorithms 2 and 4, so as to benefit from the estimated parameters of previous ALS iterations. Although this idea seems natural, the update of delay parameters would then require to design discrete search algorithms having a more complex structure than greedy search algorithms. Another perspective of this work would be to estimate the number of sources: strategies such as Akaike information criterion or Bayesian information criterion [START_REF] Stoica | Model-order selection: a review of information criterion rules[END_REF] could be a first step to explore. Finally, a third perspective would be to deal with sources whose shape parameters evolve through the mixtures.