Fabien Lauer
email: fabien.lauer@esstin.uhp-nancy.fr

René Vidal

Gérard Bloch
email: gerard.bloch@esstin.uhp-nancy.fr

A product-of-errors framework for linear hybrid system identification

Keywords: switched systems, piecewise affine systems, linear hybrid systems, system identification, global optimization, large-scale problems, robustness to outliers

 and support vector regression (Lauer and Bloch, 2008a,b) methods as particular cases. Inspired by these approaches, we then propose an optimization framework that relies on the minimization of a product of loss functions. Here, the identification problem is recast as a nonlinear and non-convex, though continuous, optimization program that involves only the model parameters as variables. As a result, its complexity scales linearly with the number of data and it can easily be solved using standard global optimization methods. Moreover, we show that by choosing a saturated loss function, such as Hampel's loss function, the algorithm can efficiently deal with noise and outliers in the data. The final result is a general framework for linear hybrid system identification that can deal efficiently with noise, outliers, and large data sets. Numerical experiments demonstrate the efficiency and robustness of the proposed approach.

INTRODUCTION

In this paper, we are concerned with the identification of a class of discrete-time hybrid systems of the ARX form y i = f λi (x i ) + e i , (1) where x i = [y i-1 . . . y i-na , u i-n k . . . u i-n k -nc+1 ] T is the continuous state (or regression vector) of dimension p containing the lagged n c inputs u i-k and n a outputs y i-k , λ i ∈ {1, . . . , n} is the discrete state (or mode) determining which one of the n subsystems {f j } n j=1 is active at time step i, and e i is an additive noise term. In particular, we concentrate on the problem of finding a hybrid model f = {f j } n j=1 of the form (1) from input-output data {(x i , y i )} N i=1 . In the following, we assume the number of modes n to be known.

Related work.

Recently, six main approaches to hybrid system identification have been developed: the clustering approach [START_REF] Ferrari-Trecate | A clustering technique for the identification of piecewise affine systems[END_REF][START_REF] Nakada | Identification of piecewise affine systems based on statistical clustering technique[END_REF], the Bayesian approach (Juloski et al., 2005b), the mixed integer programming (MIP) approach [START_REF] Roll | Identification of piecewise affine systems via mixed-integer programming[END_REF], the bounded-error approach [START_REF] Bemporad | A boundederror approach to piecewise affine system identification[END_REF], the algebraic approach [START_REF] Vidal | An algebraic geometric approach to the identification of a class of linear hybrid systems[END_REF][START_REF] Ma | Identification of deterministic switched ARX systems via identification of algebraic varieties[END_REF], and the Support Vector Regression (SVR) approach (Lauer and Bloch, 2008a;[START_REF] Lauer | From Support Vector Machines to Hybrid System Identification[END_REF]. The first four focus on the problem of PieceWise Affine (PWA) system identification, where the discrete state λ i depends on the continuous state x i . The first two approaches rely on alternating minimization and are sensitive to initialization, while the next two involve solving a combinatorial optimization problem, which can be computationally intensive. Both the bounded-error and Bayesian approaches can also be used to identify a broader class of systems, known as switched linear systems, where the discrete state evolves independently of the continuous state. The algebraic approach [START_REF] Vidal | An algebraic geometric approach to the identification of a class of linear hybrid systems[END_REF] gives a closed form solution to this latter problem. However, it is sensitive to noise compared to the clustering-based or bounded-error methods, as shown in Juloski et al. (2005a). The SVR approach provides a convenient way of dealing with noisy data by incorporating regularization into the optimization framework. However, it optimizes over a number of variables that grows with the number of data points, hence it is only applicable to small data sets. The bounded error approach suffers from the same drawback. To the best of our knowledge, other than the bounded error approach, none of the existing approaches deals explicitly with outliers in the data. We refer the reader to [START_REF] Paoletti | Identification of hybrid systems: a tutorial[END_REF] for a review and comparison of some of these methods.

Paper contributions. Our main contribution is a hybrid system identification algorithm that is computationally efficient with respect to the number of data points, and robust with respect to outliers in the data. The proposed method is based on the minimization of a product of loss functions plus a regularization term. This provides a general framework for hybrid system identification, in which any suitable loss function can be used. In particular, data corrupted by non-Gaussian noise can be treated by a suitable choice of the loss function. The interpretation of the method in a maximum likelihood framework also allows one to choose the loss functions with respect to the noise probability density function, as in standard estimation theory.

We also show that unbounded loss functions that are robust to outliers in classical linear estimation problems cannot guarantee robustness for the product-of-errors approach. Thus, we propose to use saturated loss functions, such as the Hampel's loss function [START_REF] Cichocki | Neural Networks for Optimization and Signal Processing[END_REF], which are shown to be robust in this case. Experiments show that this loss function outperforms the classical loss functions in the presence of outliers and additive noise.

The paper also proposes a reformulation of the hybrid system identification problem as an unconstrained optimization program. Though non-convex, this nonlinear program involves a low number of variables, equal to the number of parameters of the hybrid model. As we will see, its complexity scales linearly with the number of data and thus it is efficiently solvable by standard global optimization algorithms. Experiments show that the global optimum is always found in few seconds for data sets of up to hundreds of thousands of data points.

Paper organization. The paper starts by presenting a general optimization framework in §2, including its minmin and maximum likelihood interpretations ( §2.1-2.2). We also show how two existing hybrid system identification algorithms can be interpreted as particular cases of this framework ( §2.3-2.4). §3 presents the proposed product-oferrors approach to hybrid system identification, including an analysis of its robustness to outliers ( §3.1), and the global optimization method ( §3.2). The paper ends with numerical experiments in §4 and conclusions in §5.

GENERAL FRAMEWORK

The general principle behind all hybrid system identification methods is to find a collection of models {f i } n i=1 that best fit the given collection of data points {(x i , y i )} N i=1 . As a consequence, one can pose the hybrid system identification problem as an optimization problem of the form min f1,...,fn

R(f 1 , . . . , f n ) + J(f 1 , . . . , f n , {(x i , y i )} N i=1 ). (2)
The first term is called the regularizer, and measures the model smoothness or the model complexity. The second term is called the fitting error and measures the fidelity of the model with respect to the data. Data fidelity is often measured using a loss function of the error e ij = y i -f j (x i ) incurred by assigning the ith data point to the jth model.

The methods included in the framework differ in the choice of the regularizer, in the choice of the loss function, or in how these loss functions are combined across different models to form the fitting error. Before delving into the details, we first review some of the well known loss functions used in estimation theory and shown in Fig. 1.

Squared loss function. The squared loss function, l(e) = e 2 , (3) is perhaps the most well known. As will be seen in §2.3, the algebraic method [START_REF] Vidal | An algebraic geometric approach to the identification of a class of linear hybrid systems[END_REF] can be interpreted in the proposed framework on the basis of this loss function.

Absolute loss function. Another well known loss function is the absolute loss, which given by l(e) = |e|. In classical linear estimation theory, this loss function offers a certain robustness to outliers compared to the squared loss function. As will be seen in §3.1, this property does not hold for the proposed product-of-errors approach to hybrid system identification.

ε-insensitive loss function. The ε-insensitive loss function, defined by [START_REF] Vapnik | The Nature of Statistical Learning Theory[END_REF] for Support Vector Regression (SVR), l(e) = max(0, |e|ε),

(5) builds a tube of insensitivity of radius ε, inside which the errors are meaningless. Errors larger than ε are penalized linearly, which ensures a certain robustness to outliers in linear estimation theory. The ε-insensitive loss function has been used for hybrid system identification in Lauer and Bloch (2008a,b) as a relaxation to the bounded-error approach [START_REF] Bemporad | A boundederror approach to piecewise affine system identification[END_REF].

Min-min error estimator

The min-min estimator (MME) assigns sample (x i , y i ) to the submodel f j that gives the best estimate

f j (x i ) of y i : λi = arg min j=1,...,n l(y i -f j (x i )), i = 1, . . . , N, (6)
where l is, for example, one of the loss functions above. Thus, the fitting error to be minimized is the sum of the errors made after assigning point i to submodel λi , i.e., min f1,...,fn

J MME = N i=1 min j=1,...,n l(y i -f j (x i )) . (7)
One of the main difficulties with this MME framework is that it involves optimization over both discrete and continuous variables. To see this, notice that one can rewrite the optimization problem in (7) as

min {fj },{βij} N i=1 n j=1 β ij l(y i -f j (x i )), (8)
where β ij ∈ {0, 1} and j β ij = 1. The discrete variables β ij encode the assignment of points to submodels, while the continuous variables encode the parameters of each submodel. One way to solve this mixed (discrete and continuous) program is to use alternating minimization: given the submodels, compute the assignment of points to submodels using (6); and given the assignments, compute one submodel for each group of points. This approach is indeed effective when the submodels are linear and the squared loss (3) is used, because the estimation of each submodel is a linear system identification problem. However, this approach is sensitive to initialization. Note that problem (8) can also be solved by using mixed integer programming techniques, as proposed in [START_REF] Roll | Identification of piecewise affine systems via mixed-integer programming[END_REF] for hinging hyperplane models. These latter optimization techniques can guarantee to find the global minimum, but, due to their high complexity, they can only be used in practice for small data sets.

Maximum likelihood of the most likely estimator

The Maximum Likelihood (ML) approach consists in finding the model f that most likely generated the data. The Maximum Likelihood of the Most Likely (MLML) estimator assigns each training sample (x i , y i ) to the most likely submodel f j , i.e. the one with maximal likelihood of the sample w.r.t. f j . This leads to

f (x i ) = f λi (x i ), for λi = arg max j=1,...,n p(y i |x i , f j ). ( 9
)
As p(y i |x i , f ) = p(y i |x i , f λi ), the likelihood of a sample w.r.t. the hybrid model f can be written as

p(y i |x i , f ) ∝ max j=1,...,n p(y i |x i , f j ). ( 10
)
After replacing the maximization in ( 10) by the minimization of the negative log-likelihood, we obtain the following MLML estimator over all the samples min f1,...,fn

J MLML = N i=1 min j=1,...,n -ln p(y i |x i , f j ) . ( 11
)
Notice that the MLML problem amounts to a min-min optimization program, hence there is an equivalence between MLML and MME estimators. Indeed, choosing a particular loss function l for an MME estimator (7) corresponds to a particular choice of the noise probability density function, p(y i |x i , f j ), for an MLML estimator (11). Therefore, MLML suffers from the same drawbacks as MME. That is, the solution to problem (11) involves alternating minimization, which is sensitive to initialization, and possibly computationally costly for non-Gaussian errors.

Algebraic approach

One way of reducing the computational burden of mixed optimization is to use soft assignments β ij ∈ R + rather than hard assignments β ij ∈ {0, 1}. The algebraic approach [START_REF] Vidal | An algebraic geometric approach to the identification of a class of linear hybrid systems[END_REF] uses a clever choice of the assignments that leads to a closed form solution for linear models. Specifically, the algebraic approach is based on the observation that, in the absence of noise, y i = w T λi x i . This leads to the following hybrid decoupling constraints n j=1

(y i -w T j x i ) = 0, i = 1, . . . , N, (12)
from which one can algebraically solve for the model parameters as shown in [START_REF] Vidal | An algebraic geometric approach to the identification of a class of linear hybrid systems[END_REF]; [START_REF] Ma | Identification of deterministic switched ARX systems via identification of algebraic varieties[END_REF]. In case of noisy data, the authors solve the equations in (12) in a least squares sense by minimizing min w1,...,wn

J A = N i=1 n j=1 (y i -w T j x i ) 2 . ( 13
)
This cost function is the same as that in (8) with l(e) = e 2 , but with soft assignments 13) is a nonlinear optimization problem, the authors find a linear solution by minimizing over the coefficients of the product polynomial, rather than over the model parameters w j . This approximation results in a convex (quadratic) optimization problem, which can be solved very efficiently using linear techniques. However, this approximation comes at the cost of sensitivity to noise and outliers. Whenever pertinent, we will refer to these two approaches as linear or nonlinear algebraic approach.

β ij = k =j l(y i -w T k x i ). Though in principle (

Support Vector Regression approach

None of the approaches described so far includes regularization over the submodel parameters. A natural way of including regularization is to use SVR with ℓ 2 regularization, which in the case of linear models, f j (x) = w T j x, yields j w T j w j . Inspired by the hybrid decoupling constraint (12), the work of Lauer and Bloch (2008a) combines this regularizer with a fitting error which is essentially the product of upper bounds on the errors ξ ij ≥ |y iw T j x i |δ j , for thresholds δ j . This leads to the following optimization program min

wj ,ξij ≥0 1 n n j=1 w T j w j + C N N i=1 n j=1 ξ ij (14) -ξ ij -δ j ≤ y i -w T j x i ≤ δ j +ξ ij , i = 1, . . . , N, j = 1, . . . , n.
where C is the parameter that tunes the trade-off between the regularizer and the fitting error. A similar formulation of this method with ℓ 1 -norm regularization can be found in [START_REF] Lauer | From Support Vector Machines to Hybrid System Identification[END_REF].

However, the problem with this formulation is that it involves a number of variables equal to n × (p + N ), which grows with the number of samples N for a fixed number of modes n and a fixed number of parameters per mode p. To resolve this issue, notice that the SVR method can be reformulated in the proposed general framework by using the ε-insensitive loss function in (5). This leads to the optimization problem min w1,...,wn

1 n n j=1 w T j w j R SV R + C N N i=1 n j=1 max(0, |y i -w T j x i |-δ j ) J SV R
, which involves only n × p variables. This reformulation of the SVR approach is the basis for the product-of-errors approach, which we propose next.

PRODUCT-OF-ERRORS ESTIMATOR

In the previous section, we presented a general optimization framework for solving the hybrid system identification problem, and showed how different methods can be obtained via different choices of the regularizer and fitting error, as summarized in Table 1. Some methods, e.g., the algebraic approach, lead to simple optimization problems, but suffer from robustness issues. Other methods, e.g., the SVR approach, can be made robust to noise, but not to outliers and involve optimization over a large number of variables.

To circumvent these issues, in this section we propose a product-of-errors (PE) estimator of the form min f1,...,fn

1 n n j=1 R(f j ) R P E + C N N i=1 n j=1 l(y i -f j (x i )) J P E . ( 15
)
where R(f j ) is a regularizer for submodel f j and l is a robust loss function, to be defined below. The PE estimator has several important properties. First of all, notice that for linear models, f j (x) = w T j x, the PE estimator is as efficient as the nonlinear algebraic approach, because it Table 1. Particular choice of regularizer and loss function leading to the existing methods.

Method

Regularizer R(f 1 , . . . , fn) Fitting error J(f 1 , . . . , fn, {(x i , y i )} N i=1 )

SVR approach (Lauer and Bloch, 2008a) [START_REF] Vidal | An algebraic geometric approach to the identification of a class of linear hybrid systems[END_REF] N i=1 n j=1 (y iw T j x i ) 2 Product-of-errors estimator (this paper)

1 n n j=1 w j 2 2 N i=1 n j=1 max(0, |y i -w T j x i | -δ j ) Algebraic approach
1 n n j=1 R(w j ), for any R(w j ) N i=1 n j=1 l(y i -f j (x i )), for any l(y i -w T j x i )
optimizes over a number of variables that does not depend on the number of data points. In fact, J P E coincides with the algebraic error in (13) when l(e) = e 2 . However, the PE estimator can be made robust to outliers by a suitable choice of the loss function, as we will show in §3.1.

Notice also that the PE estimator coincides with the SVR approach when R(w j ) = w T j w j and l is chosen as the εinsensitive loss function. However, we will see in §3.1 that the product of ε-insensitive loss functions is not robust to outliers, and propose an alternative choice. Finally, although the PE estimator involves solving a nonlinear and non-convex program, we will show in §3.2 that global optimization techniques can be used to find the global minimum, under the mild requirement of continuity of the loss function.

Robustness to outliers

In this paper, we consider only outliers with respect to the output value, i.e., points with arbitrary value y i , but exact regression vector x i . For an estimator to be robust to such outliers, the effect of a single point on the estimation must be bounded. For instance, in classical linear estimation problems, the influence function of the squared loss (3), i.e. its derivative with respect to y i , is unbounded. On the other hand, the absolute loss (4) has an influence function bounded by 1 and thus leads to robust estimators. However, a simple calculation shows that this property does not hold when using a PE estimator for hybrid systems. In this case, the influence of a single point on the objective function (15) becomes

∂J P E ∂y i = n j=1 ∂l(y i -f j (x i )) ∂y i k∈{1,...,n}\j l(y i -f k (x i )), ( 16
)
where ∂l(

y i -f j (x i ))/∂y i = ∂l(y i -f j (x i ))/∂(y i -f j (x i ))× ∂(y i -f j (x i ))/∂y i = l ′ (e ij ).
Though the derivative l ′ (e) of the absolute loss is bounded by 1, the values l(y if k (x i )) cannot be bounded for an unbounded y i . Thus the influence of a single point cannot be bounded.

In order to obtain a robust PE estimator, one needs a loss function l leading to l ′ (e) = 0, for large values of the error |e|. For this purpose, we propose to use a saturated absolute loss function defined as l(e) = min(δ, |e|), (17) with derivative l ′ (e) = -1, if e ∈ [-δ, 0[, 1, if e ∈ [0, δ[, and 0 otherwise.

A similar methodology can be applied to define other robust loss functions, such as the saturated squared loss l(e) = min(δ 2 , e 2 ), also known as Talvar's loss function, or the smoother Hampel's loss function

l(e) = δ 2 /π (1 -cos(πe/δ)) , if |e| ≤ δ, 2δ 2 /π, otherwise. ( 18
)
0 0 0 0 0 0 δ 2δ 2 /π δ 2 -δ -δ δ δ δ -δ
Fig. 2. Saturated absolute, saturated squared (or Talvar's), and Hampel's loss functions l(e) versus the error e.

These loss functions are plotted in Fig. 2.

Optimization of the product-of-errors function

For linear submodels of the form f j (x) = w T j x, j = 1, . . . , n,, the product-of-error estimator is given by min w1,...,wn

1 n n j=1 R(w j ) + C N N i=1 n j=1 l(y i -w T j x i ), ( 19
)
where the parameter vectors w j to estimate are of dimension p = n a + n c (or p = n a + n c + 1 for affine submodels). [START_REF] Huyer | Global optimization by multilevel coordinate search[END_REF], that is guaranteed to converge if the objective is continuous in the neighborhood of the global minimizer. This optimizer uses only function values (when required, derivatives are estimated from these) and alternates between global and local search. The local search, done via sequential quadratic programming (SQP), speeds up the convergence once the global part has found a point in the bassin of attraction of the global minimizer.

NUMERICAL EXPERIMENTS

In this section, we present some numerical results supporting the approach. In particular, we consider the identification of a switched linear system from data with noise and outliers in §4.1, and large-scale experiments in §4.2.

All the experiments are performed on a 2.4GHz Core 2 duo laptop with 2GB of memory using Matlab. Though the MCS algorithm can deal with unbounded variables, box constraints are used to limit the search space and restrain the variables to the interval [-10, 10] (in practice, these bounds can be larger with little influence). Beside this, the default parameters of MCS are used. For all the problems, N samples are generated by

y i = θ T λi x i + v i , i = 1, . . . , N, (20)
where the θ j ∈ R p are the true parameters to recover and v i ∼ N (0, σ 2 v ) is a Gaussian noise. The methods are compared on the basis of the normalized Mean Squared Error (MSE) on the parameters, NMSE = n j=1 θ jw j 2 2 / θ j 2 2 , where the w j are the estimated parameters.

Switched linear system identification

Consider the example taken from [START_REF] Vidal | Recursive identification of switched ARX systems[END_REF]. The aim is to recover, from N = 1000 samples, the parameters θ 1 = [0.9, 1] T and θ 2 = [1, -1] T of a dynamical system, arbitrarily switching between n = 2 modes, with continuous state

x i = [y i-1 , u i-1 ] T .
Linear vs. nonlinear algebraic approach. As seen in §2.3, the algebraic method can be implemented either as the direct minimization over the model parameters, as in ( 13), or as a linear problem solved w.r.t. the product of parameters. These two algorithms are compared in Table 2, which shows the average and standard deviation of the NMSE over 100 trials with different noise sequences. These results highlight the gain in solving the problem directly for the model parameters w j versus the optimization over the product of parameters. The larger error obtained with the linear method for large noise levels is due to the fact that this method neglects one constraint: the polynomial must be factorable into the original form (12). However, the gain in NMSE obtained by solving directly for the w j comes at the cost of solving a nonlinear optimization program instead of a linear problem. This leads to computing times about 30 times larger. The data are corrupted with 20% of outliers by forcing the target outputs y i , at random time steps i, to take uniformly distributed random values in the interval [-10, 10].

Table 3 shows the resulting NMSE for the squared, absolute, saturated absolute and Hampel's loss functions with parameter δ set to δ = 2.

These results emphasize that the squared loss function cannot be used to accurately estimate the parameters in the presence of outliers. On the other hand, the absolute loss function provides a certain level of robustness even for a PE estimator and still leads to accurate estimates. However, the loss functions satisfying the robustness condition for PE estimators, i.e. l ′ (e) = 0 for |e| > δ, lead to a lower NMSE. In these experiments, the minimum NMSE is always obtained with Hampel's loss function. The benefits of using saturated loss functions compared to the absolute loss are also emphasized by the plots in Fig. 3, which show the NMSE for an increasing percentage of outliers and additive Gaussian noise (σ v = 0.3).

Large-scale experiments

Large data sets. The method is now evaluated over 100 large-scale problems, generated by (20) with uniformly distributed random parameters θ j ∈ [-2, 2] p , j = 1, 2, and regression vectors x i ∈ [-5, 5] p , a random switching sequence {λ i } and σ v = 0.2. The PE estimator (19) with squared loss is considered without regularization. For n = 2 and p = 3, Table 4 shows the resulting NMSE and the computing times, averaged over the 100 problems. Note that the algorithm, including the MCS optimization, is entirely implemented in non-compiled Matlab code. The times in Table 4 show that the proposed algorithm has a complexity scaling linearly with the number of data N . As a result, the method can be applied to large data sets with hundreds of thousands of data.

Large model structures. The computing time of the method heavily relies on the number of model parameters n × p. Thus the method may not be suitable for hybrid models with numerous modes. However, as shown by Fig. 4, the average computing time remains below 4 minutes for models with up to 40 parameters. The average is evaluated over 10 runs of the method using Hampel's loss function for N = 10 000 data points, generated by random sets of parameters.

Application to real data

The method is now applied to real data from a pickand-place machine, described in the comparison paper by Juloski et al. (2005a). In this paper, most methods (except the algebraic approach) had to subsample the data set of 60 000 samples to 750 samples in order to be applied in reasonable time. This subsampled data set is then divided in two overlapping subsets of 500 samples. The PE estimator (19) with squared loss and without regularization is applied to the first subset (with normalized output and a 1 appended to the regression vector) to build a hybrid model with 2 modes and orders n a = n c = 2. This model, obtained in 2 seconds, is then tested on the second subset, leading to a one-step-ahead prediction MSE = N i=1 (y if λi (x i )) 2 = 0.0590. When considering the entire data set, divided in two subsets of 30 000 samples, the model, obtained in 18 seconds, leads to MSE = 5.6248 × 10 -6 . The results reported in [START_REF] Ma | Identification of deterministic switched ARX systems via identification of algebraic varieties[END_REF] on this application are 0.1195 and 5.3426 × 10 -6 for the two settings, respectively. In comparison, the PE estimator leads to a lower MSE for the subsampled data set. Figure 5 shows the simulation of the model on the entire data set, where the data y i and the model output f (x i ) cannot be distinguished.

CONCLUSION

We have proposed a new framework for hybrid system identification, which includes some of the most recent approaches as special cases. We also proposed a specific optimization framework, which relies on the minimization of a product of loss functions to circumvent some computational issues. The loss function can be chosen on the basis of the noise model thanks to a maximum likelihood interpretation. In addition, we gave a condition for a loss function to be robust to outliers in the product-of-errors context. Numerical experiments showed that the resulting algorithm can identify linear hybrid systems from both large data sets and data sets corrupted with outliers.

Fig. 1 .

 1 Fig. 1. Squared, absolute, and ε-insensitive loss functions l(e) versus the error e.

Fig. 3 .

 3 Fig. 3. Average (left) and standard deviation (right) of the NMSE over 100 trials vs. the percentage of outliers.

Fig. 4 .

 4 Fig. 4. Average computing time in seconds (left) and NMSE (right) vs. the number of parameters per submodel for n modes.

Fig. 5 .

 5 Fig. 5. Model simulation on the real data.

Table 2 .

 2 Average NMSE and computing time over 100 trials.

		σv	Nonlinear algebraic	Linear algebraic
		0.00	0.00000	0.00000
		0.02	0.00000 ± 0.00000	0.00003 ± 0.00012
	NMSE	0.10	0.00007 ± 0.00007	0.00447 ± 0.02962
		0.20	0.00036 ± 0.00028	0.03968 ± 0.21762
		0.30	0.00069 ± 0.00055	0.04449 ± 0.19270
	Time (sec.)		0.32 ± 0.02	0.01 ± 0.00
	Robustness to outliers. We now study the robustness
	to outliers of the PE estimator (19) without regularization.

Table 3 .

 3 Average and standard deviation of the NMSE for a data set with 20% of outliers.

	Loss function		σv		NMSE				Time (sec.)
	Squared				0.0		0.22898 ± 0.10684		0.36 ± 0.06
	Absolute				0.0		0.00708 ± 0.01618		0.54 ± 0.02
	Saturated absolute	0.0		0.00165 ± 0.00442		0.65 ± 0.04
	Hampel's			0.0		0.00009 ± 0.00007 0.76 ± 0.21
	Squared				0.1		0.21267 ± 0.10700		0.37 ± 0.05
	Absolute				0.1		0.00276 ± 0.00801		0.61 ± 0.03
	Saturated absolute	0.1		0.00152 ± 0.00418		0.65 ± 0.03
	Hampel's			0.1		0.00013 ± 0.00010 0.68 ± 0.10
	Squared				0.2		0.18548 ± 0.10161		0.41 ± 0.07
	Absolute				0.2		0.00099 ± 0.00128		0.63 ± 0.03
	Saturated absolute	0.2		0.00033 ± 0.00035		0.69 ± 0.05
	Hampel's			0.2		0.00024 ± 0.00021 0.71 ± 0.13
	Squared				0.3		0.22177 ± 0.12003		0.40 ± 0.06
	Absolute				0.3		0.00170 ± 0.00141		0.64 ± 0.04
	Saturated absolute	0.3		0.00075 ± 0.00059		0.69 ± 0.04
	Hampel's			0.3		0.00050 ± 0.00042 0.74 ± 0.12
	0.012								0.012
		Absolute loss							Absolute loss
	0.01	Saturated absolute loss Hampel's loss					0.01	Saturated absolute loss Hampel's loss
	0.008								0.008
	0.006								0.006
	0.004								0.004
	0.002								0.002
	-10 0	0	10	20	30	40	50	60	-10 0	0	10	20	30	40	50	60

Table 4 .

 4 Average error and computing time over 100 randomly generated problems.

	Number of data N	NMSE (×10 -3 ) Time (seconds)
	10 000	0.0046 ± 0.0213	1.27 ± 0.17
	50 000	0.0010 ± 0.0012	4.77 ± 0.56
	100 000	0.0008 ± 0.0010	11.62 ± 1.92
	500 000	0.0007 ± 0.0010	63.70 ± 8.23

Available at http://www.mat.univie.ac.at/˜neum/software/mcs/.

Future work will focus on the estimation of the number of modes and of the subsystem orders in the presence of outliers, as well as on the extension of the proposed framework to nonlinear hybrid system identification, e.g. with nonlinear submodels in kernel form as studied by Lauer and Bloch (2008b). In addition, robustness to outliers in the regressors also requires further investigation.


Marion Dohen
email: dohen@icp.inpg.fr

Hélène Loevenbruck

Audiovisual Production and Perception of Contrastive Focus in French: a Multispeaker Study

This study examines the visual cues to prosodic contrastive focus in Hexagonal French and their role in visual speech perception. Two audiovisual corpora were recorded (from two male native speakers of French) consisting of sentences with a subject-verb-object (SVO) syntactic structure. Four conditions were studied: focus on each phrase (S,V,O) and broad focus. The corpora were first acoustically validated. Then lip area and jaw opening were extracted from the video. For each speaker, we identified a set of visible correlates of contrastive focus. The combined results showed that there were consistent visible articulatory correlates of contrastive focus across speakers: a) an increase in lip area and its first derivative on the focused item b) a lengthening of the focal syllables. There were also speaker-specific strategies in the amount of a) prefocal anticipation or b) post-focal hypo-articulation.

Visual only perception tests were then conducted to see if the identified correlates were valid cues in perception. They showed that contrastive focus was well perceived visually for both speakers. The scores were better for the first speaker who displayed greater focal hyper-articulation. We also found that presence and salience of the visual cues enhances perception.

'S1 said: Denis revived the good-looking mother?'

Introduction

Studies of French prosody have mainly focused on laryngeal and pulmonic correlates. A few supralaryngeal analyses exist, mostly considering tongue movements, e.g. [START_REF] Loevenbruck | An Investigation of Articulatory Correlates of the Accentual Phrase in French[END_REF], or spectral consequences of differences in articulation [START_REF] Tabain | Effects of prosodic boundary on /aC/ sequences: articulatory results[END_REF]. The few studies that have examined visual cues to prosody have focused on facial cues [START_REF] Burnham | Visual discrimination of Cantonese tones by tonal but non-Cantonese speakers and by nontonal language speakers[END_REF] such as eyebrow movements [START_REF] Granström | Prosodic Cues in Multimodal Speech Perception[END_REF] or on head movements [START_REF] Munhall | Visual prosody and speech intelligibility: Head movement improves auditory speech perception[END_REF]. Only few studies have examined visible mouth correlates [START_REF] Bernstein | Single-channel vibrotactile supplements to visual perception of intonation and stress[END_REF][START_REF] Jong | The supraglottal articulation of prominance in English: linguistic stress as localized hyper-articulation[END_REF][START_REF] Keating | Optical Phonetics and Visual Perception of Lexical and Phrasal Stress in English[END_REF] and none have done so for French. "Visible" mouth correlates include articulatory correlates such as mouth opening and durational ones, such as lengthening. The purpose of this study is to relate tonal and visual characteristics of contrastive focus in French and to tell whether the visible correlates are used in perception.

Contrastive focus is used to emphasize a word or group of words in an utterance as opposed to another. In French, it can be either syntactic ("C'est xxx qui court." It is xxx who runs.) or prosodic ("XXX F court." XXX F runs.). This study deals with audiovisual prosodic contrastive focus in French.

Experimental material

The audiovisual data

Two audiovisual corpora were recorded, which consisted of sentences with a Subject-Verb-Object structure (SVO) and with CV syllables. Each sentence was likely to be produced as a single Intonational Phrase (IP) consisting of 3 Accentual Phrases (APs). In the broad focus condition, following Jun & Fougeron's findings [START_REF] Fougeron | A Phonological Model of French Intonation[END_REF][START_REF] Fougeron | Realizations of Accentual Phrases in French Intonation[END_REF], the expected default tonal pattern is {[LHiLH*] S [LHiLH*] V [LHiLL%] O }. Sonorants were favoured in order to facilitate F0 tracking.

The audio-visual recording

Corpus 1 was recorded for speaker S1 (male) with front and profile cameras (see Figure 1) and was entirely analyzed. This led to an optimisation of the corpus and corpus 2 was thus recorded for speaker S2. For each corpus, four conditions were elicited: subject-, verb-and object-focus (narrow focus) and broad focus (neutral version). In order to trigger focus, the speakers had to perform a correction task by focusing a phrase which had been mispronounced in the prompt. The recording went as follows (where capital letters signal focus):

Audio prompt: S1 : Romain ranima la jolie maman. S2 : S1 a dit : Denis ranima la jolie maman ?

Speaker uttered: ROMAIN ranima la jolie maman. The speakers were given no indication on how to produce focus (e.g. which syllables should be accented). For speaker S1, two speaking modes were recorded: real and reiterant speech. For speaker S2, only real speech was recorded. Reiterant speech was produced by replacing all the syllables with [ma]. The purpose of reiterant speech is to compare the acoustic and articulatory features across all the syllables.

The first step was to acoustically validate the corpora. It was checked, for both speakers, that the focused utterances displayed a typical focused intonation as described in [START_REF] Dohen | Pre-focal Rephrasing, Focal Enhancement and Post-focal Deaccentuation in French[END_REF].

Measurement techniques

Figure 1 shows an example of the images that were recorded. A program designed at Institut de la Communication Parlée (ICP) [START_REF] Lallouache | Un poste Visage-Parole couleur. Acquisition et traitement automatique des contours de lèvres[END_REF][START_REF] Audouy | Traitement d'images vidéo pour la capture des mouvements labiaux[END_REF] enabled us to extract parameters describing lip shape and protrusion and jaw position from a sequence of digitalized frames. The mouth opening gesture was studied through a blue marker on the jaw (see Figure 1). The lip contour was automatically detected from the video signal and lip -height, -spreading, -area and -protrusion were derived.

Preliminary study: reiterant speech

Before studying real speech, a preliminary study was conducted on reiterant speech for speaker S1 [START_REF] Dohen | Visual perception of contrastive focus in reiterant French speech[END_REF]. The purpose was to determine a set of possible visible correlates to contrastive focus. These results showed that the large jaw opening gestures associated with high opening velocities on all the focused syllables and the long lip closure for the first segment of the focused group could be interpreted as a set of visual cues to the perception of focused reiterated [ma] sequences. Additional cues may be prefocal lengthening and post-focal hypo-articulation.

A visual only perception experiment showed that the visual cues described above are used for the perception of contrastive focus in French for reiterant speech.

Production studies 4.1. Preliminary analysis of the problem

Possible articulatory correlates

There are many possible visible articulatory correlates: jaw opening, lip -height, -area, -spreading, -protrusion, etc. The problem is to identify the one(s) which will vary the most significantly across conditions and the most invariantly across syllables. In our preliminary study [START_REF] Dohen | Visual perception of contrastive focus in reiterant French speech[END_REF] we found that the main articulatory consequence of contrastive focus is hyperarticulation. Hyper-articulation can be achieved in various ways, including increase in the amplitude of lip and/or jaw opening and closing movements, increase in lip spreading or narrowing. The parameter affected by hyper-articulation varies, depending on which syllable is uttered: for a hyperarticulated /a/ the mouth will be more opened thus the lip opening and the jaw opening will therefore be larger, for a hyper-articulated /i/, lip spreading, but not lip height will increase, and for a hyper-articulated /u/, lip protrusion will increase but not lip height nor spreading. The parameters which are most likely to be affected by hyper-articulation are thus lip height (LH), lip spreading (LS) and lip protrusion (LP). The lip area parameter (LA) takes into account the variations of both LH and LS. The articulatory parameters studied were thus LP, LA and LA's first derivative.

Possible durational correlates

The major durational correlates of focus identified in the study of reiterant speech were: focal lengthening, prefocal lengthening and what was called "lenghtening of initial lip closure" for the first [ma] in the focalized phrase. Similarly, in this study, focal and prefocal duration were measured, as well as the duration of the first phoneme of the focused sequence. This last parameter will thereafter be referred to as "first segment duration". In so doing , we wanted to find out if the lengthening of the initial lip closure measured for reiterant speech (see 3) was only an artifact of the syllable used or a general correlate of contrastive focus in French.

Measurements

All the maxima of LA and LP were detected. The duration of all the syllables were also computed. As explained before, the fact that real speech is studied here induces a great deal of variability. Even if the LA and LP parameters will account for most hyper-articulation strategies, there still remains a comparison issue. Lip area is indeed not comparable from an /a/ to an /i/ and the same can be put forward concerning duration. In order to compare data across the corpus, a "normalization" had to be performed. Therefore, the values of all the maxima (resp. of the corresponding syllable durations) for each parameter were divided by the mean value of the maxima of that parameter (resp. the mean value of the durations) for both broad focused utterances. All the broad focus maxima (resp. durations) therefore correspond to the value 1 and for the other focus types a value above 1 implies an increase of the considered parameter and a value below 1 a decrease of the considered parameter.

Results for corpus 1 (speaker S1)

Articulatory measurements

Figure 2 shows the grand mean of the "normalized" values of LA over each syntactic phrase and over all the identical syntactic phrases of the corpus. For example, the 1st column was computed in the following way: all the peaks of lip area were detected and "normalized" for the subject of the broad focused utterances, the means over each subject were then calculated and the means of these means were plotted. The mean increase of LA from a broad focus condition to a focused condition is of 48.7% (significant: p<0.05). The mean increase of LA's first derivative is of 49.8% (significant: p<0.05). Concerning the pre-focal sequence, we found a mean increase of LA (resp. LA's 1 st derivative) of 23.3% (resp. 13.3%). As for the post-focal sequence, we found a mean decrease of LA (resp. LA's 1 st derivative) of 2.2% (resp. 7.9%), this was not statistically significant. LP was not studied for this speaker because there were not enough protruded syllables in this corpus.

Duration measurements

Focal lengthening: The mean duration of the focal syllables is significantly higher (p<0.05) than the duration of the same syllables in the broad focused condition. The mean lengthening from broad to narrow focus is of 33.6%.

Prefocal lengthening:

The duration of the last syllable of a phrase was measured as significantly higher (p<0.05) when the following phrase was focused: +19.7%.

First segment lengthening:

The first segment of a phrase is significantly lengthened by 53.2% (p<0.05) when the phrase it belongs to is focused. The first segment is therefore more lengthened than the rest of the focused phrase (only 33.6%).

Conclusion: S1's focus strategy

Considering the results presented above, we can summarize S1's visible correlates to contrastive focus: -pre-focal anticipation: as had already been found and explained in the preliminary study [START_REF] Dohen | Visual perception of contrastive focus in reiterant French speech[END_REF], S1 displays an anticipation strategy i.e. he increases both duration and lip area just before focus. -focal hyper-articulation: The focal syllables display a significantly larger lip area (and lip area's 1 st derivative). These syllables are also significantly lengthened. We noted that both cues were used simultaneously: S1 does not either increase lip area or duration but both. -post-focal hypo-articulation: Unlike what had been found in the preliminary study [START_REF] Dohen | Visual perception of contrastive focus in reiterant French speech[END_REF], S1 does not seem to hypoarticulate the post-focal sequence. The duration of the post focal syllables does not significantly change. The mean increase of LA from a broad focus condition to a focused condition is of 8.5% (significant: p<0.05). The mean increase of LA's first derivative is of 14.9% (significant: p<0.05).

Concerning the pre focal sequence, we found a mean decrease of LA of 3.3% (not statistically significant) but an increase of LA's 1 st derivative of 31.5%. As for the post focal sequence, we found a mean decrease of LA (resp. LA's 1 st derivative) of 22.2% (resp. 22.3%). S2 also increases lip protrusion by 11.2% for the focused syllables.

Duration measurements

Focal lengthening: The mean lengthening of the syllables from the broad focus case to the focus case is of 20.4% (significant: p<0.05).

Prefocal lengthening:

The duration of the syllable preceding a focused phrase did not significantly change when the following phrase was focused (decrease of 2%; not significant).

First segment lengthening:

The mean lengthening of the first segment of a focused phrase is of 20.4%. This exactly corresponds to the focal lengthening (20.4%). For speaker S2, the first segment is thus not more lengthened than the rest of the focused phrase.

Conclusion: S2's focus strategy

Considering the results presented above, we can summarize S2's visible correlates to contrastive focus: -pre-focal anticipation: the results show that S2 does not develop an anticipation strategy (no rise in lip area or duration just before focus). -focal hyper-articulation: The focal syllables display a significantly larger lip area (and lip area's 1 st derivative) and lip protrusion. These syllables are also significantly lengthened. We noted that S2 increases both LA and duration in only 40% of the cases and increases only one of the two parameters in 40% of the cases. In 20% of the cases he increases neither. -post-focal hypo-articulation: the post-focal sequence displays an important decrease in lip area and its first derivative. The duration of the post-focal sequence however does not significantly change.

Comparison between S1 and S2

The rise for lip area and duration in S2 is not as important as for S1 (LA: S1: 48.7% S2: 8.5%; LA's 1 st derivative: S1: 49.8% S2: 14.9%; duration: S1: 33.6% S2: 20.4%).

Perception studies

Visual only perception tests were conducted to check if the visible correlates identified above are used for perception.

Description of the experiments

The participants were told that they would be witnessing a conversation between two speakers. The first speaker would pronounce an utterance which they would first hear (audio prompt). They were told that one element (Subject, Verb or Object) in this sentence was misunderstood by the second speaker, who would therefore repeat the sentence as a question. This question would neither be heard nor seen by the participants. The first speaker would then repeat the sentence and put focus on the misunderstood phrase. The participants saw a video recording of that speaker but heard no sound.

Below is an example of how the test went: Speaker 1 (audio only): Romain ranima la jolie maman. Speaker 2 (nothing): Denis ranima la jolie maman ? Speaker 1 (video only): ROMAIN ranima la jolie maman. The participants were told that, in some cases, there was no misunderstanding (corresponding to a broad focus case). They were asked to determine which phrase (S, V, O or broad) had been misunderstood and thus focused. The participants used a highlighter pen to mark the constituent they perceived as focused on an answer sheet presented as below and highlighted the empty cell when they perceived broad focus.

Romain

ranima la jolie maman.

Test 1: speaker S1

We used four sentences from the corpus for their nearly balanced structures (almost the same number of syllables in S, V and O): s2, s4, s6 and s7. A total of 32 sentence pairs (1 pair: audio only unfocused utterance and visual only focused utterance) were available (4 sentences, 4 focus conditions, 2 repetitions). Five tests consisting of five random combinations of the 32 pairs were presented to each participant. These five tests were the same for all participants but the presentation order was different. Therefore, each person was presented with a total of 160 pairs of sentences. Both front and profile views were shown at the same time. A total of 33 native speakers of French participated in the experiment.

Test 2: speaker S2

We used nine sentences from the corpus for their nearly balanced structures: from s9 to s17. A total of 72 sentence pairs were available (9 sentences, 4 focus conditions, 2 repetitions). Two tests consisting of two random combinations of the 72 pairs were presented to each participant. The participants were tested on both views (front and profile) separately: some were presented with the first test front and the others with the first test profile and vice-versa for the second test. Therefore, each person was presented with a total of 144 pairs of sentences.

A total of 27 native speakers of French participated in the experiment.

Results

The results showed that the participants successfully perceived the focus through the visual modality alone. For test 1, the percentage of correct answers was of 71.45% and for test 2 it was of 43% (chance level for both tests: 25%). The scores for the test using speaker S1 were better, a finding which we expected since hyper-articulation was more salient for that speaker. However both scores are well above chance. A detailed analysis of the results showed that poorly perceived stimuli corresponded to unsalient visible correlates (articulatory and durational). This supports the hypothesis that the correlates perceived are those identified in the production studies. However, it was also found that the stimuli with the highest scores displayed all the correlates but none of them was either very unsalient nor very salient. This could mean that all the correlates are necessary to best identify focus even if they are not highly significant.

Discussion & Conclusions

The measurements and experiments suggest that there are lower face visual correlates of contrastive focus in French which can intervene in audiovisual speech perception. It is highly possible that more subtle facial correlates are also used in the visual perception of focus. Those could be head and/or eyebrow movements as suggested in [START_REF] Granström | Prosodic Cues in Multimodal Speech Perception[END_REF][START_REF] Munhall | Visual prosody and speech intelligibility: Head movement improves auditory speech perception[END_REF][START_REF] Keating | Optical Phonetics and Visual Perception of Lexical and Phrasal Stress in English[END_REF]. We are currently analyzing Optotrak data to assess this issue.

corpus 1 :

 1 it consists of eight sentences: s1.[Jean]S1 [veut ménager]V3 [nos jolis nouveaux navets]O7. 'Jean wants to spare our fine new turnips.' s2. [Romain]S2 [ranima]V3 [la jolie maman]O5. 'Romain revived the good-looking mother.' s3. [Mélanie]S3 [vit]V1 [les mauvais loups malheureux]O7. 'Melanie saw the unhappy bad wolves.' s4. [Véroniqua]S3 [mangeait]V2 [les mauvais melons]O5. 'Veroniqua was eating the bad melons.' s5. [Les mauvais loups]S4 [mangeront]V3 [Jean]O1. 'The bad wolves will eat John.' s6.[Mon mari]S3 [veut ranimer]V4 [Romain]O2. 'My husband wants to revive Romain.' s7.[Les loups]S2 [suivaient]V2 [Marilou]O3. 'The wolves were following Marilou.' s8.[Le beau marin]S4 [vit]V1 [Véroniqua]O4. 'The good-looking sailor saw Veronica.' corpus 2: it consists of thirteen sentences. The first four sentences correspond to s2, s4, s6 and s7 from corpus1.s9. [La nounou]S3 [mariera]V3 [Li.]O1 'The nurse will marry Li.' s10. [Le lama lent]S4 [lu]V1 [Marinella.]O4 'The slow lama read Marinella.' s11. [Marinella]S4 [va laminer]V4 [Numu.]O2 'Marinella will laminate Numu.' s12. [Lou]S1 [mima]V2 [le lama.]O3 'Lou mimed the lama.' s13. [Le nominé]S4 [lu]V1 [les longs mots.]O3 'The nominee read the long words.' s14. [La nounou]S3 [vit]V1 [Lou.]O1 'The nurse saw Lou.' s15. [Les loups]S2 [mimaient]V2 [Marilou.]O3 'The wolves mimed Marilou.' s16. [Lou]S1 [ramena]V3 [Manu.]O2 'Lou gave a lift back to Manu.' s17. [Li]S1 [ralluma]V3 [les moulinets.]O4 'Li lighted the wheels again.'

Figure 1 :

 1 Figure 1: (left) Video signal recorded (profile and front views) and (right) symbolic representation of the lip area parameter.

Figure 2 :

 2 Figure 2: Grand mean of (left) the max of LA (cm²) and (right) the max of LA's 1 st derivative over S, V & O.

4. 3 . Results for corpus 2 (Figure 3 :

 323 Figure 3: Grand mean of (left) the max of LA (cm²) and (right) the max of LA's 1 st derivative over S, V & O.

Figure 3

 3 Figure 3 corresponds to Figure2for speaker S2. The mean increase of LA from a broad focus condition to a focused condition is of 8.5% (significant: p<0.05). The mean increase of LA's first derivative is of 14.9% (significant: p<0.05). Concerning the pre focal sequence, we found a mean decrease of LA of 3.3% (not statistically significant) but an increase of LA's 1 st derivative of 31.5%. As for the post focal sequence, we found a mean decrease of LA (resp. LA's 1 st derivative) of 22.2% (resp. 22.3%). S2 also increases lip protrusion by 11.2% for the focused syllables.

Acknowledgements

We thank G. Rolland for designing and recording corpus 1 and C. Savariaux and A. Arnal for their technical help. We also thank M.-A. Cathiard, J.-L. Schwartz and P. Welby for their comments on our work.