Krylov-Subspace-Based Order Reduction Methods Applied to Generate Compact-Electro-Thermal Models for MEMS

Keywords: compact thermo-electric model, order reduction, Krylov-subspace, Arnoldi process, Lanczos algorithm

The high power dissipation density in today's miniature electronic/mechanical systems makes on-chip thermal management crucial. In order to achieve quick-to-evaluate, yet accurate electro-thermal models, needed for the thermal management of microsystems, model order reduction is necessary. In this paper, we use Krylov-subspace methods for the order reduction of a electro-thermal MEMS model, illustrated by a novel type of micropropulsion device. Comparison between different moment-matching algorithms including a new two-sided Arnoldi algorithm, is performed.

INTRODUCTION

Modeling of thermo-electric processes becomes increasingly important for a variety of applications, including power transistors, thick-film circuits, prediction of electrostatic discharge, hotplate sensors, oxide-confined verticalcavity lasers, and so on. Therefore, it is necessary to develop a electro-thermal model which computes the dependence between power dissipation and temperature distribution over the device. Moreover, such a heat transfer analysis needs to be done quickly in response to every design alteration. The model must also provide good accuracy in order to return precise temperature values.

The main problem of electro-thermal modeling is that ordinary differential equation (ODE) systems resulting from finite element, finite difference or some other spatial discretization method can easily reach an order of 100 000 or more. Even with the increasing speed of modern computers it is not possible to perform effective system-level simulation without simplification or model reduction.

Conventionally, the reduction of thermo-electric models for micro-electronic and MEMS devices is performed through a lumped-element decomposition of the model followed by parameter fitting [START_REF] Alessandro | A Critical Review of Thermal Models for Electro-thermal Simulation[END_REF]. Such a non-automatic approach requires the designer to choose the corect reduced model structure without strict guidelines, and to perform a time-consuming parametrization including one or more simulations of the full-scale model. Moreover, sufficient accuracy is provided only for a limited parametric domain.

In order to achieve both efficiency and accuracy in thermal management of microsystems, we propose a different, automatic order reduction approach, suitable for linear electro-thermal models, and based on Krylov-subspace methods. In adittion to the Arnoldi [START_REF] Ramaswami | Automatic Generation of Small-Signal Dynamic Macromodels from 3D Simulation[END_REF] and Lanczos [START_REF] Freund | Krylov-subspace Methods for Reduced Order Modeling in Circuit Simulation[END_REF] algorithm reported previously, we introduce a new method called twosided Arnoldi to find two basis necessary for projection and calculating the reduced order model. All three algorithms were tested and compared with regard to their computational complexity, accuracy of approximation, numerical stability, preservation of the stability and passivity of the original system, invariance properties and approximation of the complete output.

As a MEMS test case a new class of high energy actuator, which integrates solid fuel with three silicon micromachined wafers [START_REF] Rudnyi | Solid Propellant Microthruster: Theory of Operation and Modelling Strategy[END_REF], was used. This microthruster is ignited by passing an electric current through a polysilicon resistor embedded in the membrane, as shown in Fig. 1. The present work considers the initial heating phase of the fuel, right up to the onset of ignition, described through the following equations: [START_REF] Alessandro | A Critical Review of Thermal Models for Electro-thermal Simulation[END_REF] where is the thermal conductivity, is the specific heat capacity, is the mass density, is the temperature distribution, is the heat generation, is the spartially varying electric current density vector and is the specific electric conductivity.

We use a two dimensional axi-symmetric model, which after the finite element (FE) based spatial discretization of the governing equations (1) (assuming that the heat generation is uniformly distributed within the heating area) results in a linear system of about 1000 ordinary differential equations (ODEs) of the form: [START_REF] Ramaswami | Automatic Generation of Small-Signal Dynamic Macromodels from 3D Simulation[END_REF] where are the global heat conductivity and heat capacity matrix, are the temperature (state), the load and the output vector respectively and is the dimension of the system. The electric current through the heater with electric resistivity is the input to the system.The equation ( 2) represent a Single-Input-Single-Output (SISO) system. The present work also considers a special case when is an identity matrix, that is, , which we call a Single-Input-Complete Output (SICO) system.

MODEL ORDER REDUCTION

Most of the practical work in model reduction of large linear dynamic systems has been tied with moment matching of the transfer function via Krylov subspaces by means of either the Arnoldi or Lanczos process. They define a projection from the high dimensional space of the original model (2) to a lower dimensional space and thereby define the reduced order model. By applying the transformation where and r < n , to system (2) and then multiplying the state equation by transpose of some matrix , a model with reduced order r can be found:

(3)

The reduced system matrices and load vector are computed by: (4)

The key question is: how to find V and W ?

Let us rewrite the system (2) as: [START_REF] Salimbahrami | Invariance Properties of Krylov Subspace Methods in Linear Model Order Reduction[END_REF] with and . As already mentioned, the basic idea behind the Krylov-subspace-based algorithms is to write down the transfer function of (5) in the frequency domain using a Taylor series in the Laplace variable s around :

where is called the i th moment, and then to find a much lower order system (of the same form as ( 5)) whose transfer function has the same moments as up to some degree. Due to the numerical instability of ( 6), the moments are not computed explicitly. Instead, a Krylov subspace (of the dimension r ) defined as: [START_REF] Silveira | A Coordinate-transformed Arnoldi Algorithm for Generating Guaranteed Stable Reduced-Order Models of RL C Circuits[END_REF] is used. The vectors that span the subspace are called the basic vectors.

Arnoldi Algorithm

In case of Arnoldi algorithm . This algorithm generates a set of orthonormal vectors (with length one and orthogonal to each other) which simultaneously represent a basis for the given Krylov-subspace [START_REF] Silveira | A Coordinate-transformed Arnoldi Algorithm for Generating Guaranteed Stable Reduced-Order Models of RL C Circuits[END_REF]. They are saved as columns of the matrix V . This means that . The algorithm further generates an upper Hessenberg matrix ( ) , related to the system matrix A as follows:

(8)

The matrix H A can be considered as an orthogonal projection of the matrix A onto the Krylov-subspace [START_REF] Silveira | A Coordinate-transformed Arnoldi Algorithm for Generating Guaranteed Stable Reduced-Order Models of RL C Circuits[END_REF], and it is equal to the system matrix of the reduced system A r .. It can be proved that the first moments of and match [START_REF] Salimbahrami | Invariance Properties of Krylov Subspace Methods in Linear Model Order Reduction[END_REF].

Lanczos algorithm

For the Lanczos algorithm . The columns of the matrix V form a basis (but not an orthonormal one) of the Krylov-subspace [START_REF] Silveira | A Coordinate-transformed Arnoldi Algorithm for Generating Guaranteed Stable Reduced-Order Models of RL C Circuits[END_REF], which is also called an input Krylovsubspace. The columns of the matrix W form a nonorthonormal basis of the output Krylov-subspace defined as:

(9) κ T ∇ ( ) ∇• Q ρC p t ∂ ∂T - + 0 Q j 2 σ ---- = , = κ C p ρ T Q j σ Q C [ ]T ˙K [ ]T + FI t ( ) 2 R = y E T T ⋅ = K [ ] C [ ] R n n × ∈ , T t ( ) F E , R n ∈ , n I t ( ) R E y T = T V [ ] T r ⋅ = V R n r × ∈ W R n r × ∈ W T C [ ]VT r ˙WT K [ ]VT r + W T FI t ( ) 2 R = y r E T VT r ⋅ = C [ ] r W T C [ ]V K [ ] ; r W T K [ ]V F r ; = W T F = = A [ ]T ˙T bu t ( ) + = y c T T ⋅ = A K [ ] - 1 -C [ ] b K [ ] - 1 -F = , = c E = s 0 0 = G s ( ) { } m { } i s i i 0 = ∞ ∑ - = m { } i c T A ( ) i b = G r s ( ) { } G s ( ) { } K r A b , { } span b Ab … A r 1 -b , , , { } = V W = V T V ⋅ I = a ij 0 i j , ( ) j i - 1 > , ∀ , = H A R r r × ∈ V T A V ⋅ ⋅ H A = r G r s ( ) { } G s ( ) { } V W ≠ K r A * c , { } span c A * c … A * ( ) r 1 - c , , ,       =
where A * is the conjugate transpose of the matrix A. The basis V and W of the subspaces ( 7) and (9) are biorthogonal, which means that . The algorithm generates further a tridiagonal matrix T L related to the original system matrix as:

(10)

The matrix T L can be considered as an oblique projection of the matrix A onto the input Krylov-subspace [START_REF] Silveira | A Coordinate-transformed Arnoldi Algorithm for Generating Guaranteed Stable Reduced-Order Models of RL C Circuits[END_REF], and it is equal the system matrix of the reduced system A r . In this case the first moments of and match [START_REF] Salimbahrami | Invariance Properties of Krylov Subspace Methods in Linear Model Order Reduction[END_REF].

Two-Sided Arnoldi Algorithm

The two-sided Arnoldi algorithm computes matrix V as an orthonormal basis for the input Krylov-subspace [START_REF] Silveira | A Coordinate-transformed Arnoldi Algorithm for Generating Guaranteed Stable Reduced-Order Models of RL C Circuits[END_REF], and matrix W as an orthonormal basis for the output Krylov-subspace (9), by twice using the one-sided Arnoldi algorithm (section 2.1). This means that and . The reduced system is computed as in (3), and the first moments of and match as in the Lanczos algorithm. Moreover, the transfer function of the reduced model gained from the two-sided Arnoldi method equals the one gained from the Lanczos algorithm [START_REF] Salimbahrami | Invariance Properties of Krylov Subspace Methods in Linear Model Order Reduction[END_REF]. Further details can be found in [START_REF] Salimbahrami | Two-sided Arnoldi Algorithm and Its Application in Order Reduction of MEMS[END_REF].

RESULTS

In Fig. 2 a relative error between the full and the 5th order reduced model, using both the one-sided Arnoldi and two-sided Arnoldi algorithm, is shown. The results for Lanczos algorithm are identical as for two-sided Arnoldi.

For the microthruster model, the simple SISO setup for both the one-sided and the two-sided Arnoldi algorithm was sufficient to approximate not only a single output response but also the transient thermal response in all the finite element nodes of the microthruster (SICO setup). Fig. 3 shows the mean relative difference for all the nodes between the full-scale and the reduced different order models. Hence it was possible, after the simulation of the reduced model, to recover the full solution (for all the 1071 nodes) by applying the projection . In Table 1 the maximal relative error for a single output node (node 1 in Fig 1) e 1 and a maximal mean square relative differences for approximating the complete output e 2 for different orders of the reduced model are given. The results for the Lanczos algorithm are identical as for the two-sided Arnoldi.

DISCUSSION AND CONCLUSION

Accuracy of approximation: As mentioned in chapter 2, the one-sided Arnoldi algorithm matches only r moments of and , whereas the two-sided Arnoldi and Lanczos algorithm match 2r moments. Hence, the twosided Arnoldi and Lanczos algorithm are optimal in a sense that they match as many moments as there are free coefficients in the reduced order transfer function (r poles and r zeroes). This provides a unique reduced order model, i. e., invariance regarding the representation and the realization of the original system, but does not guarantee the preservation Fig. 2 Relative error between the full and the 5th order reduced model for a single node (node 1 in Fig 1).

W T V ⋅ I = W T A V ⋅ ⋅ T L = 2r G r s ( ) { } G s ( ) { } V T V ⋅ I = W T W ⋅ I = 2r G r s ( ) { } G s ( ) { } 0 0.
MSRD in % T V [ ] T r ⋅ = G r s ( ) { } G s ( ) { }
of stability and passivity. In many applications it is better to trade some of the "optimality" of two-sided methods to gain guaranteed stable and passive reduced order models, for example by using the one-sided Arnoldi reduction algorithm.

Preservation of the stability and passivity of the original system: In circuit simulation, reduced-order modeling is mostly applied to large passive linear subcircuits (such as RLC networks), and preservation of passivity is crucial for the stability of the simulation of the whole circuit. Unfortunately, the reduced models gained by two-sided Arnoldi and Lanczos algorithm are not passive and stable in general [START_REF] Freund | Krylov-subspace Methods for Reduced Order Modeling in Circuit Simulation[END_REF]. In the microthruster example, the reduced model of order 9 is unstable. It has been shown that the coordinate transformed one-sided Arnoldi algorithm on the other side, can generate guaranteed stable reduced-order models [START_REF] Silveira | A Coordinate-transformed Arnoldi Algorithm for Generating Guaranteed Stable Reduced-Order Models of RL C Circuits[END_REF].

Computational complexity: One disadvantage of the Arnoldi method is that each new Arnoldi vector (new column of matrix V) should be orthogonal to all previously generated vectors. This means that the computational cost for orthogonalisation over the r steps of algorithm grows as O(2r 2 n) with the dimensions n of the full space and r of the Krylov-subspace. Additionally r steps of the Arnoldi procedure require r matrix-vector products at the cost of 2rNz(A), where Nz(A) is a number of nonzero elements of A (Nz(A)=n 2 for a dense matrix). Thus, on average the computational costs for (one-sided) Arnoldi algorithm grows as O(2r 2 n + 2rNz(A)) . For the two-sided Arnoldi algorithm the costs are double. The Lanczos algorithm has smallest computational costs for orthogonalisation. In each step it is necessary to deal with just two previously generated vectors (matrix T L is tridiagonal), which makes the orthogonalisation costs over the r steps of the algorithm grow only as O(16rn). Together with matrix-vector products, the Lanczos algorithm´s costs are of O(16rn+4rNz(A)). Hence, for large r ´the Lanczos algorithm is the fastest of three. However, it is numerically less stable than both the Arnoldi and the twosided Arnoldi process.

Numerical stability: Because of rounding errors in all three algorithms, the column vectors of basis V and W for the Krylov-subspaces ( 7) and (9) may become non-orthogonal. How quickly this happens depends on the chosen dimension r of the Krylov-subspaces. For the one-and twosided Arnoldi algorithm each new vector should be orthogonal to all previously generated vectors so that the rounding errors accumulate slower than by Lanczos algorithm, where each new vector is orthogonalised only with respect to the last two generated vectors.

Invariance properties: Changing the representation (multiplying the state equation with some matrix) or the realization (change of state vector) of the original system does not change the input-output behavior of the reduced order models generated by two-sided methods [START_REF] Salimbahrami | Invariance Properties of Krylov Subspace Methods in Linear Model Order Reduction[END_REF]. In case of one-sided Arnoldi algorithm these invariance properties are not given, due to matching only r moments which is half of the number of unknowns in . Approximation of the complete output: In general the two-sided Arnoldi and Lanczos algorithms produce reduced-order models which are "optimized" for particular output(s). This is due to their taking into account the output vector c by using a basis W of the output Krylov-subspace (9). It means that, even if in the case of the microthruster model two-sided Arnoldi properly describe the complete state vector (see e 2 in Table 1), in general this may not work. The one-sided Arnoldi algorithm on the other hand, doesn´t take into account the output vector at all, so we can expect the approximation of the complete output also in the general case.

Conclusion: Our results show that, for moderate values of r the two-sided Arnoldi is definitely better than Lanczos algorithm. Stability and passivity properties of two-sided methods however, need to be further researched. Meanwhile, we consider the one-sided Arnoldi algorithm as an acceptable tool for compacting linear electro-thermal models.

Fig. 1

 1 Fig. 1 Microthruster Structure.

Table 1 :

 1 Results for one-sided and two-sided Arnoldi

		0
	error in %	-4 -2	One-sided Arnoldi Two-sided Arnoldi
		-6
		-8
			01 0.02 0.03 0.04 0.05 0.06 0.07
			time in s

ACKNOWLEDGMENTS

This work is partially funded by the EU through the project MICROPYROS (IST-1999-29047), partially by the DFG project MST-Compact (KO-1883/6) and partially by an operating grant of the University of Freiburg.

G r s

( ) { }


Ioannis Z Emiris
email: emiris@di.uoa.gr

Guillaume Moroz
email: guillaume.moroz@inria.fr

The assembly modes of rigid 11-bar linkages

Keywords: 11-bar linkage, assembly modes, polynomial system, mixed volume, distance matrix, cross entropy, simulated annealing

Designing an m-bar linkage with a maximal number of assembly modes is important in robot kinematics, and has further applications in structural biology and computational geometry. A related question concerns the number of assembly modes of rigid mechanisms as a function of their nodes n, which is uniquely defined given m. Rigid 11-bar linkages, where n = 7, are the simplest planar linkages for which these questions were still open. It will be proven that the maximal number of assembly modes of such linkages is exactly 56.

The rigidity of a linkage is captured by a polynomial system derived from distance, or Cayley-Menger, matrices. The upper bound on the number of assembly modes is obtained as the mixed volume of a 5 × 5 system. An 11-bar linkage admitting 56 configurations is constructed using stochastic optimisation methods. This yields a general lower bound of Ω(2.3 n ) on the number of assembly modes, slightly improving the current record of Ω(2.289 n ), while the best known upper bound is roughly 4 n . Our methods are straightforward and have been implemented in Maple. They are described in general terms illustrating the fact that they can be readily extended to other planar or spatial linkages.

This version (2017) typesets correctly the last figure 5 so as to include all 28 configurations modulo reflection.

Introduction

Rigid mechanisms (or linkages) constitute an old but still very active area of research in mechanism and linkage theory, e.g. [START_REF] Collins | Forward kinematics of planar parallel manipulators in the Clifford algebra of P 2[END_REF][START_REF] Faugère | The combinatorial classes of parallel manipulators combinatorial classes of parallel manipulators[END_REF][START_REF] Walter | On a 9-bar linkage, its possible configurations and conditions for paradoxical mobility[END_REF][START_REF] Walter | A spatial 9-bar linkage, possible configurations and conditions for paradoxical mobility[END_REF][START_REF] Wunderlich | Gefärlice annahmen der trilateration und bewegliche afchwerke i[END_REF] as well as computational geometry and structural bioinformatics, e.g. [BS04, EM99, JRKT01, Hav98, TD99].

A given linkage may be represented by a graph with edge set E of lengths l ij ∈ R + , for (i, j) ∈ E. An assembly mode, or Euclidean embedding, in R d is a mapping of its vertices to a set of points in R d , such that l ij equals the Euclidean distance between the images of the i-th and j-th vertices, for (i, j) ∈ E. Euclidean embeddings impose no requirements on whether the edges cross each other or not. A linkage is (generically) rigid in R d if and only if, for generic edge lengths it can be embedded in R d in a finite number of ways, modulo rigid motions. A graph is minimally rigid if and only if it is no longer rigid once any edge is removed.

V [ 7 ] V [ 6 ] V [ 5 ] V [ 4 ] V [ 3 ] V [ 2 ] V [ 1 ]
Figure 1: An 11-bar linkage.

Let us focus on planar linkages and the associated graphs. A graph is called Laman if and only if the number of edges is 2n -3, where n is the number of vertices and, additionally, all of its vertex-induced subgraphs on 3 ≤ k < n vertices have ≤ 2k -3 edges. This is essentially the Grübler-Kutzbach-Chebychev formula on the degrees of freedom for mechanical linkages, e.g. [START_REF] Angeles | Rational Kinematics[END_REF]. It is a fundamental theorem that the class of Laman graphs coincides with the generically minimally rigid planar linkages [START_REF] Laman | On graphs and rigidity of plane skeletal structures[END_REF]. The problem studied in this paper is to compute the number of distinct planar assembly modes of rigid mechanisms, up to rigid motions. In particular, the maximal number of assembly modes of an 11-bar linkage (Figure 1) will be presented for the first time.

Existing work

The algebraic approach is to define a well-constrained polynomial system expressing the length constraints, such that the system's real solutions correspond precisely to the different assembly modes. When defining a straightforward system such as (1), where the unknowns are the nodes' coordinates, all nontrivial equations are quadratic. For planar linkages, there are 2n -4 equations hence, by applying the classical Bézout bound on the number of common roots, we obtain 4 n-2 . It is indicative of the hardness of the problem that efforts to substantially improve these bounds have failed.

Today, the best general upper bound is roughly 4 n-2 / π(n -2). This was obtained using determinantal varieties defined by distance matrices [START_REF] Borcea | The number of embeddings of minimally rigid graphs[END_REF]. Straightforward application of mixed volumes (discussed in Section 3) yields an upper bound of 4 n-2 [ST10].

The best general lower bounds are 24 (n-2)/4 2.21 n and 2 • 12 (n-3)/3 2.289 n /6, obtained by a caterpillar and a fan1 construction, respectively [START_REF] Borcea | The number of embeddings of minimally rigid graphs[END_REF]. Both bounds are based on the Desargues graph (Figure 2), which admits 24 assembly modes.

In applications, it is crucial to know the number of assembly modes for specific (small) values of n. The most important result in this direction was to show that the Desargues linkage, also known as the planar parallel robot, admits precisely 24 assembly modes in the plane [START_REF] Borcea | The number of embeddings of minimally rigid graphs[END_REF][START_REF] Hunt | Structural kinematics of in parallel actuated robot arms[END_REF][START_REF] Gosselin | Solutions polynomiales au problème de la cinématique directe des manipulateurs parallèles plans à trois degrés de liberté[END_REF][START_REF] Lazard | The (true) Stewart platform has 12 configurations[END_REF]. Moreover, the K 3,3 , or 9-bar, linkage admits exactly 16 assembly modes in the plane [START_REF] Walter | On a 9-bar linkage, its possible configurations and conditions for paradoxical mobility[END_REF]. This paper demonstrated such a linkage through an elaborate construction which was not required in our case.

All existing bounds for n ≤ 10 are found in

Our contribution

To upper bound the number of planar rigid mechanisms, we explore adequate polynomial systems leading to tight root bounds. Modeling physical systems by appropriate polynomial systems is a deep and hard question, with a wide range of applications in different fields.

We employ two powerful algebraic tools for defining polynomial systems and for counting the system's common roots. For the former, we use distance matrices, also known as Cayley-Menger matrices, which contain all known and unknown distances between the graph's nodes. The signs of the matrix minors capture rigidity and embeddability in Euclidean spaces, as described in Section 3. Our results indicate that such matrices are advantageous to using the coordinates' formulation, such as system (1), when constructing the polynomial system.

Our second tool is the mixed volume of a well-constrained polynomial system, which exploits the sparseness of the equations. The mixed volume bounds the number of common roots, by Bernstein's Theorem 1, as described in Section 3. This bound is never larger than Bézout's, and is typically much tighter for systems whose equations do not contain all possible terms for a given total degree. It turns out that such are the systems encountered here and could be of wider interest in similar enumeration problems.

Moreover, mixed volume is a bound so it only needs to consider which are the nonzero coefficients, without considering specific values. This is in sharp contrast to solving methods often used for root counting, where one assigns random values to the coefficients.

We thus obtain a tight bound for linkages with n = 7 nodes by appropriately formulating the polynomial system based on distance matrices. We compute the mixed volume of a 5 × 5 polynomial system, which equals 56. This is tight since we demonstrate a construction with as many assembly modes, see Figure 5. We also show how our stochastic optimisation methods are generalizable for studying further planar or spatial linkages. Our construction also yields a general lower bound of Ω(2.3 n ) on the number of configurations, thus improving the available bound of Ω(2.289 n ), whereas the best available upper bound is O(4 n / √ n). The rest of the paper is structured as follows. Section 2 presents rigid mechanisms and how our study focuses on a specific 11-bar linkage. Section 3 presents our algebraic tools and obtains the upper bound of 56 for 7 nodes and 11 bars. Section 4 actually constructs a linkage with 56 assembly modes.

We conclude with open questions.

Rigid linkages

This section studies the number of assembly modes of planar rigid linkages or, equivalently, rigid graphs and their number of embeddings in R 2 .

Such linkages admit inductive constructions that begin with a triangle, followed by a sequence of so-called Henneberg steps. There are two types of such steps, each adding one new vertex and increasing the total number of edges by 2. A graph is Laman, or the associated mechanism is rigid, if and only if it can be constructed by a sequence of the corresponding Henneberg steps. By exploiting this fact, all rigid graphs in R 2 , but also in R 3 , were constructed using the Henneberg steps in [START_REF] Emiris | Algebraic methods for counting Euclidean embeddings of rigid graphs[END_REF], and classified up to graph isomorhism.

Let us consider the Henneberg steps defining Laman graphs, each adding a new vertex and a total number of two edges. A Henneberg-1 (or H 1 ) step connects the new vertex to two existing vertices. A Hennenerg-2 (or H 2 ) step connects the new vertex to 3 existing vertices having at least one edge among them, and this edge is removed. Both steps are illustrated in Figure 3. We represent each Laman graph by s 4 . . . , s n , where s i ∈ {1, 2}; this is known as its Henneberg sequence. Note that this sequence is by no means unique. A Laman graph is called H 1 if and only if it can be constructed using only H 1 steps; it is called H 2 otherwise.

Since two circles intersect generically in two points, a H 1 step at most doubles the number of assembly modes and this is tight, generically. In Table 1 most lower bounds follow from this fact. The lower bound for n = 9 follows from the Desargues fan [START_REF] Borcea | The number of embeddings of minimally rigid graphs[END_REF].

For a Laman graph on 6 vertices, a tight uppper bound of 24 follows by examining the only three possibilities: the graph is either H 1 , it is K 3,3 , or is the Desargues graph. Now, H 1 graphs on 6 vertices have exactly 16 assembly modes. The K 3,3 graph has precisely 16 assembly modes [START_REF] Walter | On a 9-bar linkage, its possible configurations and conditions for paradoxical mobility[END_REF], a fact first conjectured in [START_REF] Wunderlich | Gefärlice annahmen der trilateration und bewegliche afchwerke i[END_REF]. The Desargues graph has precisely 24 assembly modes: the upper bound was first shown in [START_REF] Hunt | Structural kinematics of in parallel actuated robot arms[END_REF] and proven more explicitly, along with the lower bound, in [START_REF] Gosselin | Solutions polynomiales au problème de la cinématique directe des manipulateurs parallèles plans à trois degrés de liberté[END_REF]. This fact was independently rediscovered in [START_REF] Borcea | The number of embeddings of minimally rigid graphs[END_REF] and proven via coupler curves.

In [START_REF] Emiris | Algebraic methods for counting Euclidean embeddings of rigid graphs[END_REF] further upper bounds were obtained from the following family of simple systems, employed in order to express graph embeddability in R 2 . Let x i , y i denote the coordinates of the i-th vertex, where the l ij are the given lengths:

x i = a i , y i = b i , i = 1, 2, a i , b i ∈ R, (x i -x j ) 2 + (y i -y j ) 2 = l 2 ij , (i, j) ∈ E -{(1, 2)} (1)
Vertex (x 1 , y 1 ) was fixed to (0, 0), to discard translations, and (x 2 , y 2 ) to (1, 0) to remove rotations and scaling, assuming without loss of generality that the corresponding edge exists. Then the mixed volume of the system gave the upper bounds in the table which, unfortunately, are loose, even after removing roots at infinity [START_REF] Emiris | Algebraic methods for counting Euclidean embeddings of rigid graphs[END_REF]. Below we overcome this limitation by considering polynomial systems derived from distance matrices. We now study the case n = 7. If a H 1 step is applied to any linkage with n = 6, the resulting linkage with n = 7 admits exactly 48 assembly modes, which is reflected to the corresponding lower bound in Table 1. To maximize the number of assembly modes, we shall apply a H 2 step to any graph with n = 6. By checking graph isomorphisms and taking into account the various symmetries, it was shown [START_REF] Emiris | Algebraic methods for counting Euclidean embeddings of rigid graphs[END_REF] that there are only 3 relevant graphs to be considered. These are obtained by a H 2 step applied to the Desargues graph as follows, where notation refers to Figure 2: We remove edge (4, 5), and add edges (4, 7), (5, 7), and

• (1, 7), or • (3, 7), or • (6, 7).
The first case corresponds to the topology that shall be studied extensively in the sequel, since it leads to the maximum number of 56 assembly modes. The linkage is shown in Figure ref11barLinkage and in the next section we will proove that 56 is the maximum number of assembly modes possible for such a linkage. The other two linkages admit at most 44 and 48 assembly modes, respectively, hence they are not studied any further. These upper bounds are obtained as mixed volumes of polynomial systems, as explained in Section 3.

Algebraic methods for the upper bound

This section discusses polynomial systems, describes distance matrices, and derives an upper bound on the number of assembly modes.

Mixed volume

We first discuss multivariate polynomial systems and introduce sparse elimination theory in order to exploit their structure and sparseness; for details, see [START_REF] Cox | Using Algebraic Geometry[END_REF][START_REF] Emiris | Sparse Elimination and Applications in Kinematics[END_REF]. Classical elimination theory characterizes every polynomial by its total degree. For a well-constrained system of polynomial equations, the classical Bézout bound on the number of isolated roots equals the product of the polynomials' total degrees. One disadvantage of this bound is that it counts projective roots and hence increases when there are roots at projective infinity.

In sparse (or toric) elimination theory, a polynomial is characterized by its support. Given a polynomial f in n variables, its support is the set of exponents in N n corresponding to nonzero terms (or monomials). The Newton polytope of f is the convex hull of its support and lies in R n . Consider (Newton) polytopes P i ⊂ R n and parameters λ i ∈ R, λ i ≥ 0, for i = 1, . . . , n. We denote by λ i P i the corresponding scalar multiple of P i . Consider the Minkowski sum of the scaled polytopes

λ 1 P 1 + • • • + λ n P n ∈ R n ; its Euclidean volume is a homogeneous polynomial of degree n in the λ i . The coefficient of λ 1 • • • λ n is the mixed volume of P 1 , . . . , P n . If P 1 = • • • = P n ,
then the mixed volume is n! times the Euclidean volume of P 1 .

We now focus on the topological torus C * = C -{0} in order to state Bernstein's root bound in terms of mixed volume.

Theorem 1 [Ber75] Let f 1 = • • • = f n =
0 be a polynomial system in n variables with real coefficients, where the f i have fixed supports.

The number of isolated common solutions in (C * ) n is bounded above by the mixed volume of the Newton polytopes of the f i . This bound is tight for a generic choice of coefficients of the f i 's.

One alternative to using general bounds such as mixed volume is to manipulate a system, with coefficients chosen randomly, so as to bound the number of common real roots, e.g. by means of Gröbner bases or homotopy continuation. We actually used Gröbner bases for all systems studied in this paper, where we have used random coefficients. In particular, for the linkage analyzed below, we computed the total-degree Gröbner basis of system (1) with random distances in Maple, and obtained the Hilbert polynomial of the corresponding ideal. This is a univariate polynomial where, by setting its variable to 1, one deduces an upper bound on the number of complex common roots. This was indeed 56, which coincides with mixed volume. However, this only offers an indication on the number of roots, not an actual bound, since it depends on the choice of coefficients.

The advantage of mixed volume is that it treats entire classes of systems defined by their nonzero terms, without considering specific coefficient values. Lastly, a tight mixed volume implies that one can solve the system efficiently either by sparse resultants, e.g. [START_REF] Canny | A Subdivision-Based Algorithm for the Sparse Resultant[END_REF], or by sparse homotopies, e.g. [START_REF] Verschelde | Algorithm 795: PHCpack: a general-purpose solver for polynomial systems by homotopy continuation[END_REF].

Distance geometry

The theory of distance geometry has been well developed, e.g. [START_REF] Blumenthal | Theory and Applications of Distance Geometry[END_REF][START_REF] Schönberg | Remarks to M. Frechet's article "Sur la définition axiomatique d'une classe d'espaces vectoriels distanciés applicables vectoriellement sur l'espace de Hilbert[END_REF], with several applications especially in structural bioinformatics, e.g. [START_REF] Emiris | Computer algebra methods for studying and computing molecular conformations[END_REF][START_REF] Havel | Distance geometry: Theory, algorithms, and chemical applications[END_REF].

Definition 1 Given an n-vertex graph, the corresponding distance or Cayley-Menger matrix is a symmetric (n + 1) × (n + 1) matrix B, indexed by i = 0, . . . , n, such that B(0, i) = 1 for i ≥ 1, B(i, i) = 0 for i ≥ 0, and B(i, j) equals the squared distance between vertices i, j, for j > i ≥ 1.

To be more precise, to obtain the Cayley-Menger matrix of the given n vertices, we must multiply B by -1/2. But for testing embeddability, the two formulations are equivalent, as seen in the following theorem, due to the work of Cayley and Menger, see [START_REF] Blumenthal | Theory and Applications of Distance Geometry[END_REF][START_REF] Schönberg | Remarks to M. Frechet's article "Sur la définition axiomatique d'une classe d'espaces vectoriels distanciés applicables vectoriellement sur l'espace de Hilbert[END_REF]. Let D(i 1 , . . . , i k ) denote the (k + 1) × (k + 1) diagonal minor of B indexed by rows and columns 0, i 1 , . . . , i k , where i

1 < • • • < i k ∈ {1, . . . , n}.
Theorem 2 Suppose we are given a matrix B of the form specified in Definition 1. Then B corresponds to a (complete) graph embeddable in R d , if and only if 1. for k = 2, . . . , d + 1 and any {i 1 , . . . , i k } ⊂ {1, . . . , n}, it holds (-1) k D(i 1 , . . . , i k ) ≥ 0, and 2. rank(B) = d + 2.

The second condition, due to Cayley, yields a (large) number of equalities, which are typically not independent.

The first condition, due to Menger, yields inequalities: For k = 2, it expresses the fact that all entries must be non-negative. For k = 3, it captures the triangular inequality. If we apply it, without loss of generality, to indices 1, 2, 3, using c ij = l 2 ij for the given entries, then the condition states that D(1, 2, 3) ≤ 0, where this quantity equals -(l 12 + l 13 + l 23 )(l 12 + l 13 -l 23 )(l 12 + l 23 -l 13 )(l 13 + l 23 -l 12 ). This vanishes precisely when the corresponding points are collinear, whereas D(1, 2, 3) < 0 when the points define a triangle. Equivalently, D(i, j, k) ≤ 0 can be written l ik + l jk ≥ l ij for all triplets i, j, k ∈ {1, . . . , n}. For k = 4 the condition captures the tetrangular inequality.

Upper bound

We shall employ the rank condition in order to derive equality constraints on the unspecified distances of our linkage. Here is the matrix, where the c ij = l 2 ij correspond to the fixed distances, and x ij are the unspecified distances.

v 1 v 2 v 3 v 4 v 5 v 6 v 7 v 1 v 2 v 3 v 4 v 5 v 6 v 7             0 1 1 1 1 1 1 1 1 0 c 12 c
           
By Theorem 2, any 5 × 5 minor of this matrix must vanish, which yields polynomial equations on the variables x ij . There are 7 3 = 35 such minors, each in 2 to 4 variables. Among these polynomials, no 4 × 4 subsystem exists that corresponds to a rigid mechanism, as can be verified by checking Laman's condition. This is indispensable for the system to have a finite number of solutions.

However, it is possible to find certain 5 × 5 subsystems whose subgraph is Laman and, moreover, uniquely define the configuration of the overall linkage (Figure 4). One of these systems has 4 bivariate equations and one trivariate equation, and is defined by taking the following diagonal minors: (2) defines a system of 3 quadratic and two cubic equations in x 15 , x 16 , x 35 , x 45 , x 67 ; the cubics are the first and last polynomials. Now, this system's mixed volume turns out to be 56, computed using the software from [START_REF] Emiris | Efficient incremental algorithms for the sparse resultant and the mixed volume[END_REF] 3 . This bounds the number of (complex) common system's solutions, hence the number of (real) assembly modes.

               D(4,
Notice that this bound does not take into account solutions with zero coordinates, in other words some zero length. However, a linkage has assembly modes with some zero length only when the input bar lengths form a singular set, in the sense that they would satisfy a non-generic algebraic dependency. For example, by letting some input distance be exactly 0, some 11-bar linkage may theoretically have infinitely many configurations. However, generically, it is impossible to have such a linkage.

Thus, the mixed volume bound guarantees that a true 11-bar linkage can have at most 56 assembly modes.

Lower bound

To prove that 56 is a lower bound on the number of configurations, it is sufficient to exhibit a linkage with 56 configurations. Unfortunately, we have 11 design parameters. Thus our search space is homeomorphic to R 11 + . Even if we reduce the search space to the integers between 1 and 100 for each parameter, an exhaustive search would lead us to consider 10 22 linkages.

Another approach was used in [START_REF] Lazard | The (true) Stewart platform has 12 configurations[END_REF] to find a planar parallel robot with 12 real solutions. The authors found a singular positions with many real solutions, and use deformations to get the desired result. The advantage of this approach is that the space of singular configurations has a lower dimension than the full design space. However, in our case, the dimension was still too high and this approach could not allow us to conclude.

Stochastic methods

Random sampling is a common approach to search in high dimensional space is random sampling. The system at hand being homogeneous, we can specify the last coordinate l 10 to 100 without restriction of generality. Moreover, in the following, we restrict our domain to integer coordinates.

Our goal is to maximize the function:

N : N 10 → N
(l 0 , ..., l 9 ) → 0, if system (2) has infinitely many solutions. Number of solutions of system (2), otherwise.

A first idea is to sample random points using a Gaussian law centered on an arbitrary initial point. This naive approach did not allow us to find a linkage with 56 configurations, but we could finally find such a point by using more sophisticated methods.

Evaluation of the objective function

The cost of evaluating the objective function is the bottleneck of the stochastic methods. The function N needs to be evaluated several hundreds of time before converging toward a point where it is maximal.

In our case, evaluating N means solving a system of polynomial equations with finitely many solutions, and counting the number of its real positive roots. Experimentally, we use the Maple function RootFinding[Isolate]. This function computes isolating boxes around the real solutions of the system to solve. The underlying algorithm uses Grbner bases and Rational Univariate Representation.

The choice of the system modeling the 11-bar manipulator is critical. A naive modelling yields a system with 14 variables, one for each of the coordinates of the 7 vertices, and 11 equations, one for each link in the manipulator. In Maple, solving this system is roughly 10 times slower than solving System (2). That can be explained partly by the fact that, generically, the time complexity to solve a 0-dimensional system is roughly exponential in the number of variables, e.g. [START_REF] Cox | Using Algebraic Geometry[END_REF][START_REF] Lazard | Solving systems of algebraic equations[END_REF].

Simulated annealing

The Monte-Carlo simulation and its simulated annealing variant [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF] are the most well spread stochastic optimisation methods. These methods have already been used in robotics for path planning in [START_REF] Kavraki | Probabilistic roadmaps for path planning in high-dimensional configuration spaces[END_REF]. The convergence of these methods has been well studied and the simulated annealing simulations has been successfully used in different field such as biology and chemistry. We implemented this method to maximize the function N .

The Monte-Carlo simulation depends on a parameter T called temperature. In the simulated annealing variant, the temperature T is decreased according to a specific schedule at each step. The optimal way to decrease the temperature depends on each problem. In our case, we chose arbitrarily a linear cooling schedule ( [START_REF] Szu | Fast simulated annealing[END_REF][START_REF] Nourani | A comparison of simulated annealing cooling strategies[END_REF]). The corresponding pseudo-code is summarized in Algorithm 1.

Cross entropy method

A new method was introduced recently by Rubinstein in [START_REF] Rubinstein | Optimization of computer simulation models with rare events[END_REF] for the simulation of rare events. This method is especially well-suited for combinatorial and continuous optimization ( [START_REF] Rubinstein | The cross-entropy method for combinatorial and continuous optimization[END_REF]). The idea of this approach consists in minimizing the distance between specific probability laws appearing during the computation (see [START_REF] De Boer | A tutorial on the crossentropy method[END_REF] for more details). For our problem, we use the scheme developed for continuous multi-extremal function optimization presented in [START_REF] Kroese | The cross-entropy method for continuous multi-extremal optimization[END_REF]. The function N is not continuous, but behaves smoothly: generically, if we modify the lengths of a given linkage by small enough values, the number of configurations is increased or decreased only by 2. For our problem it yields Algorithm 2.

Results

In Table 2, we compare the results of three stochastic methods used to optimize the function N with respect to the lengths of the 11-bar linkage.

The first column shows the result of the direct random sampling. Each line correspond to 600 evaluation of N on points chosen according to a Gaussian law of center 100 and standard deviation 100. This operation has been run 10 times, yielding a total of 6000 evaluations of N on random points. However, this did not allow us to find a linkage with 56 configurations.

The second approach is a Simulated Annealing simulation. We ran it 10 times, limiting the number of evaluations of N to 600. One of our simulation returned a set of link lengths yielding a linkage with 56.

Finally, the last column reports the results of the Cross Entropy Method. Each line correspond to a simulation stopped after 600 evaluations of N . This approach was the most efficient. Four simulations out of ten returned a linkage with 56 configurations.

The results of our stochastic methods are summarized in Table 2. In particular, we found that the manipulator with the following design parameters has exactly 56 assembly modes. Some of its assembly modes are shown in Figure 5.

Our simulations show that the set of 11-bar linkages with 56 assembly modes has a non-zero volume. With the direct random sampling, we did not find any linkage with more than 48 assembly modes. This indicates that linkages with bar lengths between 0 and 300 have a very small probability to have 56 assembly modes. Our experiments show that the Cross Entropy method is a general approach well suited to generate such configurations. Moreover, it can be used easily for other mechanisms, such as larger planar or spatial linkages for example. A variant of this method can also be used to compute an estimation of the probability to find a linkage with a maximal number of assembly mode. The constant under the exponent is 4 √ 28 2.3003 while in the previous lower bound the constant was 3 √ 12 2.2894.

Further work

Undoubtedly, the most important and oldest problem in rigidity theory is the full combinatorial characterization of rigid graphs in R 3 . We believe that one can extend our methods to spatial linkages. We expect our approach leads to an algorithmic process for obtaining good algebraic representations of the enumeration problem at hand, namely low mixed volumes for the systems derived from distance matrices. One issue is that the number of equations produced is quite large, with algebraic dependancies among them. The question becomes then to choose the best well-constrained system. The structures studied in this paper are point-and-bar structures; they generalize to body-and-bar, where edges can be connected to different points of a rigid body. It is known that a body-and-bar structure in R d is rigid if and only if the associated graph is the edge-disjoint union of d+1 2 spanning trees [START_REF] Whiteley | Generating isostatic frameworks[END_REF]. Body-and-bar structures are our next object of study.

V[6] V[5] V[4] V[3] V[2] V[1]

[ 7 ] V [ 6 ] V [ 5 ] V [ 4 ] V [ 3 ] V [ 2 ] V [ 1 ]

  13 c 14 x 15 x 16 c 17 1 c 12 0 c 23 x 24 c 25 x 26 x 27 1 c 13 c 23 0 x 34 x 35 c 36 x 37 1 c 14 x 24 x 34 0 x 45 c 46 c 47 1 x 15 c 25 x 35 x 45 0 c 56 c 57 1 x 16 x 26 c 36 c 46 c 56 0 x 67 1 c 17 x 27 x 37 c 47 c 57 x 67 0

  5, 6, 7)(c 46 , c 47 , c 56 , c 57 , x 45 , x 67 ) = 0 D(1, 4, 6, 7)(c 14 , c 17 , c 46 , c 47 , x 16 , x 67 ) = 0 D(1, 4, 5, 7)(c 14 , c 17 , c 47 , c 57 , x 15 , x 45 ) = 0 D(1, 2, 3, 5)(c 12 , c 13 , c 25 , c 23 , x 15 , x 35 ) = 0 D(1, 3, 5, 6)(c 13 , c 36 , c 56 , x 15 , x 16 , x 35 ) = 0

l 0 =

 0 180 l 1 = 70 l 2 = 200 l 3 = 205 l 4 = 210 l 5 = 205 l 6 = 80 l 7 = 200 l 8 = 70 l 9 = 200 l 10 = 100

Algorithm 1

 1 Simulated annealing with linear cooling to maximize the number of solutions of System (2) T 0 ← 4 maxStep ← 1000 σ ← [10, ..., 10] point ← {l 0 = 100, ..., l 9 = 100} value ← N (l 0 , ..., l 9 ) for n from 1 to maxStep while value < 56 do newP oint ← GaussianN eighbour(point, σ) newV alue ← N (l 0 , ..., l 9 ) treshold ← Uniform random value in [0,1] point ← newP oint value ← newV alue end if T ← T 0 (1 -n/maxStep)

Figure 2 :

 2 Figure 2: The Desargues linkage or, equivalently, a planar parallel robot.

Figure 3 :

 3 Figure 3: The planar Henneberg steps; the bottom vertex is new.

V

Figure 4 :

 4 Figure 4: When the edge lengths [15], [16], [35], [45] and [67] are specified as a solution of System (2), the 11-bar linkage has only one assembly mode.

Figure 5 :

 5 Figure5: 28 of the 56 configurations of our linkage. The other 28 can be deduced by symmetry with respect to the horizontal axis.

Table 1 :

 1 Bounds for the number of embeddings of rigid graphs with n ≤ 10 [ETV09].

	n =	3 4 5 6	7	8	9	10
	upper 2 4 8 24 64 128 512 2048
	lower 2 4 8 24 48 96 288 576

Table 1 from [ETV09] 2 .

This slightly corrects the exponent of the original statement.

see[START_REF] Emiris | Mixed volume and distance geometry techniques for counting Euclidean embeddings of rigid graphs[END_REF] for an update.

http://www.di.uoa.gr/∼emiris/

Acknowledgements. I.Z. Emiris is partially supported by FP7 contract PITN-GA-2008-214584 SAGA: Shapes, Algebra, and Geometry. He also acknowledges fruitful discussions with Elias Tsigaridas. Both authors were inspired to work on this problem at the Workshop on Discrete and Algebraic Geometry in September 2010, at Val d'Ajol, France.

Improvement on the general lower bound

This specific linkage allows us to improve slightly the lower bound on the maximal number of assembly modes of a planar linkage with n pivot joints. We use the same fan construction as in [START_REF] Borcea | The number of embeddings of minimally rigid graphs[END_REF].

We construct a linkage with n nodes composed of k sub-mechanisms as follows. Each of the k sub-mechanisms is a 11-bar linkage with the topology we have considered so it has 56 assembly modes, and 28 when one triangle is fixed. All sub-mechanisms share this triangle so there are no degrees of freedom other than those within each sub-mechanism, and the overall mechanism is rigid. The total number of assembly modes equals the product of remaining assembly modes per sub-mechanism, since these sub-mechanisms can be configured independently. Hence the overall number is 28 k when one triangle is specified, and 2 • 28 k in total. For n pivot joints, we let k be the greatest integer ≤ n-3 4 .

Algorithm 2 Cross Entropy Method to maximize the number of solutions of System (2)


Amel Fraisse
email: fraisse@limsi.fr

Patrick Paroubek

Twitter as a Comparable Corpus to build Multilingual Affective Lexicons

Keywords: Affective Lexicon, Comparable Corpus, Sentiment Analysis

The main issue of any lexicon-based sentiment analysis system is the lack of affective lexicons. Such lexicons contain lists of words annotated with their affective classes. There exist some number of such resources but only for few languages and often for a small number of affective classes, generally restricted to two classes (positive and negative). In this paper we propose to use Twitter as a comparable corpus to generate a fine-grained and multilingual affective lexicons. Our approach is based in the co-occurence between English and target affective words in the same emotional corpus. And it can be applied to any number of target languages. In this paper we describe the building of affective lexicons for seven languages (en, fr, de, it, es, pt, ru).

Introduction

Research in Sentiment Analysis and Opinion Mining, has flourished in the past years. The growing interest in processing emotions and opinions expressed in written text is motivated by the birth and rapid expansion of the Social Web that made it possible for people all over the world to share, comment or consult content on any given topic. In this context, opinions, sentiments and emotions expressed in Social Media texts have been shown to have a high influence on the social and worldwide economic behavior. In spite of the growing body of research in the area in the past years, dealing with affective phenomena in text has proven to be a complex and interdisciplinary problem that remains far from being solved. As any emergent field, its challenges include the need to develop linguistic resources to perform computational tasks. In our case, we are interested in the sentiment classification task which is performed either with statistical approaches or with lexicon-based approaches. In the two cases, the lack and the scarcity of affective lexicons present a real issue for sentiment analysis system. Multilingual affective lexicons are central components for cross-lingual sentiment analysis systems. Their manual construction is a hard, long and costly process. While often it is impossible to consider for most under-resourced languages because of the scarcity or even lack of experts. Existing affective lexicons are always monolingual and often developed for English. Furthermore, many of these lexicons are very simple, i.e. they consist of a list of words divided into only two classes : positive and negative. To our knowledge, there is no fine grained affective and multilingual lexicons. Most previous work addressing the problem of bilingual lexicon extraction are based on parallel corpora. However, despite serious efforts in the compilation of corpora [START_REF] Armstrong | A presentation of mlcc : Multilingual corpora for cooperation[END_REF], [START_REF] Church | Introduction to the special issue on computational linguistics using large corpora[END_REF], to our knowledge, there is no available affective parallel corpus for the field of sentiment analysis. On the other hand, with the rapidly growing volume of resources on the Web, the acquisition of non-parallel texts is usually much easier. Thus, as mentioned by [START_REF] Rapp | Identifying word translations in nonparallel texts[END_REF] and [START_REF] Rapp | Automatic identification of word translations from unrelated english and german corpora[END_REF] it would be desirable to have an approach that can extract lexicons from comparable or even unrelated texts. In this paper, we propose to use Twitter as a comparable corpus to extract multilingual affective corpus. Our approach is motivated by the fact that, nowadays, social media user's and in particular twitter users' express and share their sentiments, opinions and emotions on a variety of topics and discuss current issues over the world. In fact, many people can talk about the same event and describe their emotional state triggered by this event in different languages. Hence, Twitter could be considered as a comparable corpus as we could group tweets (messages written by users) by emotion/opinion/sentiment expressed in different languages. We have tested our approach to build seven affective lexicons for English, French, German, Spanish, Italian, Portuguese and Russian.

Related Work

There are two ways to cover the lack of sentiment analysis resources. The first way is to create manually a lexicon in a source language as [START_REF] Bradley | Affective norms for english words (anew)[END_REF] who developed the Affective Norms of English Words (ANEW) which is a set of normative emotional ratings for 1034 English words. And then localize the source lexicon into target languages. [START_REF] Redondo | The spanish adaptation of anew (affective norms for english words)[END_REF] have adapted the ANEW into Spanish, [START_REF] Vo | The berlin affective word list reloaded (bawl-r)[END_REF] localized it into German. This approach requires human translators to ensure the quality of the localized resource and therefore is cost expensive and not scalable. [START_REF] Strapparava | Wordnet-affect : an affective extension of wordnet[END_REF] developed the WordNet Affect which is a manually created extension of the Word-Net, including a subset of synsets suitable to represent affective concepts correlated with affective words. The second approach is automatic construction of a lexicon. The most common method is bootstrapping. This method starts with seed words with a known polarity (e.g. good, happy, wonderful for a positive class, bad, sad, terrible for a negative class). Next, the seed words are used to find related words and assign them the same class or estimate their po-larity. [START_REF] Qadir | Bootstrapped learning of emotion hashtags #hashtags4you[END_REF] present a bootstrapping algorithm to automatically learn English twitter hashtags that convey emotion. [START_REF] Mohammad | Emotional tweets[END_REF] use the pointwise mutual information to measure the association between a word and a given emotion. So he builds a word emotion association lexicons which are lists of words and associated emotions. For example, the word victory may be associated with the emotions of joy and relief. [START_REF] Pak | Construction d'un lexique affectif pour le franc ¸ais à partir de twitter[END_REF] proposed to use Twitter to collect a dataset of emotional texts in French. Using the collected dataset, they estimated the affective norms of words present in the corpus and built a polarity classifier. Both for manual and automatic approaches, existing affective lexicons are always monolingual.

Word-Opinion/Sentiment/Emotion association lexicon

In a previous work [START_REF] Fraisse | Toward a unifying model for opinion, sentiment and emotion information extraction[END_REF], we have presented 20 semantic categories including all types of emotions, sentiments and opinions. Each semantic class correspond to one type of emotion/sentiment/opinion and is referred to by means of a multi-word label that regroups various subjective words generally associated to one of the various sentiments contained in the considered class ( For each label of the Table 1 and for each language, we whish to extract the associated lexicon. Table 2 illustrate an example of comparable tweets in four languages ; the four tweeter talked about the same topic violence in ukraine expressing the same emotion Sadness in different languages. So, based on such data our approach aims to extract, accross different languages, and for each affective label the 2 associated lexicon (Table 3).

Our approach for multilingual affective lexicon construction

Hashtags are a distinctive characteristic of tweets (Jackiewicz and Vidak, 2014). They are a community created convention for providing meta-information about a tweet.

Hashtags are made by adding the hash symbol # as a prefix to a word. Thus, a hashtag is simply a way for people to search for tweets that have a common topic. In general, the tweeter (one who tweets) use emotion-word hashtag, to notify others of the emotions associated with the message he or she is tweeting. Consider the tweet bellow :

Oh okay all the people I fancy are taken ...that's cool watch them be happy as I sit in a corner and cry #sad

The tweeter has used the emotion word hashtag #sad, to convey that he or she is sad. And as English is considered the reference language on the Web, the tweeter use generally the emotion word hashtag in their native languages and give the corresponding English one as shown in the following French tweet :

Je suis vraiment #triste aujourd'hui #sad.

So, our approach is based on the co-occurence between the English and the target emotion word hashtags in the tweets. To achieve this, we proceed in two steps ; firstly we construct emotional corpora in the following seven languages : English, French, Spanish, German, from the collected corpora based on the co-occurence between English and target emotional hashtags in the same affenctive class.

Corpora collection

Data collection from the Web usually involves crawling and parsing of HTML pages which is a solvable but at the same time a consuming task. In our case, collecting data from Twitter is much easier since it provides an easy and welldocumented API1 to access its content. In this work, we selected from the Table 1 the 8 prominent affective classes that are frequent in tweets : Negative surprise, Anger, Sadness, Fear, Displeasure, Boredom, Positive surprise, Pleasure and Love. For each selected class we have defined a list of English seed emotional words that are commonly used by English speakers to express their affective state on Twitter. Table 4 presents an extract of English seed emotional words that are used for the three affective classes Anger, Fear and Love. Then, we supplied the Twitter Search API with the English emotional hashtags queries and collected tweets written in their native languages and containing at least one hashtag of the English list. In fact, we noticed that when a user writes an affective tweet, he or she uses an emotional word hashtag in his or her native language and he or she, also, gives the corresponding English word. The characteristics of the gathered corpus are presented in the Table 5.

Lexicon construction

In the preprocessing of the collected corpus, we discarded tweets with the prefix Rt, RT, and rt, which indicate that the tweet that follow are re-tweets (re-postings of tweets sent earlier by somebody else). Second, we grouped the gathered tweets by language and by emotion (Figure 1). Then, for each emotion e i.e. SADNESS, PLEASURE, LOVE, etc., we extract all cooccurent hashtags and compute their correlation to e. In order to compute how much an hashtag h is correlated to an emotion e, we compute the Strength of Association the (SoA) between an hashtag h and an emotion e (Equation 1). We discarded short (less than 2 characters) and numerical hashtags.

SoA(h, e) = log ( f req(h, e) f req(h) • f req(e)

)

Where the f req(h, e) is the is the number of times h occurs in tweets belonging to the emotion e. And f req(h), f req(e) are the frequencies of h and e in the corpus.

If an hashtag appear in more than one emotion class, we associate it to the most correlated class. The size of the constructed lexicons is about 17.000 entries for the seven languages.

Conclusion

In this research we have presented a novel approach based on Twitter as a comparable corpus to extract automatically affective lexicons in seven langages (English, French, German, Italian, Spanish, Portuguese and Russian). Our approach was motivated by the fact, that non english speaker's, usually, use bilingual terms in their messages. So, we are based in the co-occurence between the English and the target affective terms to generate multilingual affective lexicons. The presented approach is generic as it could be applied for any language. Since the number of returned tweets is limited by the Twitter Search API, in a future work, we plan to use the Twitter Streaming API 2 , in order to collect a larger corpus and then obtain larger lexicons. Obtained lexicons, contains not only purely emotio- nal words but also some common-sence words that are associated to an affective class ; such as the german word Tod which is associated to the Fear affective class or the french term coeur which is assiciated to the Love class. So, for each langiage, we plan to divide the obtained lexicon into two sub-categories : purely emotional words and commonsence words.

#

  Ukraine #death toll rises as clashes continue #sad #grief. #Ukraine 60 morts aujourd'hui ! ! ! c'est vraiment #triste #chagrin Stop the #violence in #Ukraine 60 #tod heute #traurig Impresionantes imagenes de #Kiev que pasarian por fotogramas de una pelicula de guerra ! ! muchos #muertos ! ! estoy #triste

FIGURE 1 -

 1 FIGURE 1 -Extraction of Hashtags from the French corpus

TABLE 1

 1

	-uComp semantic categories of opi-
	nion/sentiment/emotion,	e=emotion,	s=sentiment,
	o=opinion, i=information, +=positive valence, -=negative
	valence.

TABLE 2 -

 2 Example of comparable tweets

	Affect. Label		Associated Words
		English French	German Spanish
		Sad	Triste	Traurig Triste
	SADNESS	Death	Mort	Tod	Muertos
		Grief	Chagrin

TABLE 3 -

 3 Multilingual affective lexicon associated to the tweets described in the Table

TABLE 4 -

 4 Seed Affective word for Anger, Fear and Pleasure affective classes

	Anger		Fear		Love
	#anger		#fear			#love
	#rage		#terror		#affection
	#irritation		#shyness		#care
	#nervousness #worry		#tenderness
	#impatience	#apprehension #fondness
	#annoyance	#terrified		#kindness
	#angry		#alarm		#attachment
	#edgy		# scare		#devotion
	#exasperated #scared		#passion
	#irritated					#envy
	#annoyed					#desire
	Affect.Cl.	En	Fr	De	It	Es	Ru	Pt
	DISCOMFORT	551	232	65	33	157	10	63
	FEAR	1677	156	123	15	488	35	124
	DISPLEASURE	1617	645	13	7	74	6	15
	SADNESS	283	211	204	209	459	110	272
	ANGER	1690	73	9	16	198	102	43
	CONTEMPT	506	606	53	43	310	69	68
	PLEASURE	2414	1952	1639	1099	2082	664	1198
		2452	434	595	632	2251	1369	933

Italian, Portuguese and Russian. Secondly, We extract affective lexicon

TABLE 5 -

 5 Number of document per affective class and per language.

TABLE 6 -

 6 The Top-10 entries of the French and German affective lexicons for the Anger, Fear, Love and Pleasure emotion classes.

	Affective Class French	German
	Anger	en colère	wütend
		fâcher	angepisst
		rage	Wut
		irriter	zerstoren
		rougir	Unterbrechung
		nervosité	Tollwut
		massacre	massaker
		énerver	Erregung
		exciter	schötteln
		furax	verärgert
	Fear	peur	angst
		terreur	terror
		violence	befürchten
		trombler	gestrandet
		mort	Tod
		terrifié	erschrocken
		appréhender achtgeben
		inquiétude	sorge
		timidité	eingeschüchtert
		anxiété	ängstlich
	Love	amour	Liebe
		Valentin	verheiratet
		coeur	verpassen
		mariage	schön
		manquer	verpassen
		aimer	lieben
		adorer	leidenschaft
		envie	Neid
		gentillesse	freundlichkeit
		affection	zuneigung
	Pleasure	heureux	Vergnügen
		content	glücklich
		génial	spielend
		bonheur	Musik
		plaisir	schön
		jouer	underschön
		vacances	Ferien
		podium	erstaunlich
		agréable	reizend
		amusant	lustig

2. https ://dev.twitter.com/docs/streaming-apis

Twitter API : https ://dev.twitter.com/docs


Céline De Looze1
email: deloozec@tcd.ie

Alain Ghio2
email: alain.ghio@lpl-aix.fr

Stefan Scherer1

Gilles Pouchoulin2

François Viallet2

AUTOMATIC ANALYSIS OF THE PROSODIC VARIATIONS IN PARKINSONIAN READ AND SEMI-SPONTANEOUS SPEECH

Keywords: prosodic variations, automatic measurements and tools, interpausal run, Parkinson disease, dysprosody

This work presents an automatic analysis of Parkinsonian speech prosody in read and spontaneous tasks. We first evaluate the reliability of existing algorithms for the automatic analysis of speech prosody and show that these tools are efficient for the analysis of Parkinsonian speech if a large amount of data is processed. We then investigate the prosodic characteristics of 25 control and 30 Parkinsonian subjects (in off (non-medicated) state) in both tasks and show that Parkinsonian read and spontaneous speech are characterized by longer pauses as well as smaller variations in pitch range, voice intensity level and articulation rate; these characteristics being all the more marked in semi-spontaneous speech. These results suggest that some linguistic functions, normally expressed via prosodic variations, may not be conveyed in Parkinsonian speech, particularly in a semi-spontaneous task.

Introduction

Parkinsonian speech associated with hypokinetic dysarthria has often been described as monotonous, monoloudness, stress-reduced, variable speech rate and marked by inappropriate silences and paroxysmal acceleration [START_REF] Canter | Speech characteristics of patients with Parkinson's disease: Intensity, pitch and duration[END_REF][START_REF] Monrad-Krohn | Dysprosody or altered melody of language[END_REF][START_REF] Darley | Motor speech disorders[END_REF][START_REF] Weismer | Acoustic description of dysarthric speech: Perception correlates and physiological inferences[END_REF][START_REF] Caekebeke | The interpretation of dysprosody in patients with Parkinson's disease[END_REF][START_REF] Ho Ak | Regulation of parkinsonian speech volume: the effect of interlocutor distance[END_REF][START_REF] Teston | L'évaluation objective de la prosodie in Les dysarthries[END_REF][START_REF] Duez | Syllable structure, syllable duration and final lengthening in parkinsonian French speech[END_REF][START_REF] Cheang | An acoustic investigation of Parkinsonian speech in linguistic and emotional contexts[END_REF]. While these descriptions were first based on perceptual assessments [START_REF] Canter | Speech characteristics of patients with Parkinson's disease: Intensity, pitch and duration[END_REF] [START_REF] Monrad-Krohn | Dysprosody or altered melody of language[END_REF], they were later augmented by instrumental analyses [START_REF] Teston | L'évaluation objective de la prosodie in Les dysarthries[END_REF][START_REF] Duez | Syllable structure, syllable duration and final lengthening in parkinsonian French speech[END_REF][START_REF] Cheang | An acoustic investigation of Parkinsonian speech in linguistic and emotional contexts[END_REF]. Acoustic analyses, which provide more objective measurements, have been generally considered as "augmenting the impressionistic nature of perceptual description" [START_REF] Kent | Voice quality measurement[END_REF]. This approach is, however, not an easy task. It requires transcriptions and annotations at different levels of analysis (e.g. phonemes, syllables, prosodic units), which can be both time-consuming and error-prone when manually carried out and when processing large amounts of data. An alternative practical approach is to use automatic tools. This requires however that they are reliable for the analysis of pathological speech, with remaining inaccuracies being compensated by the large amount of data processed. In terms of prosodic analyses, the definition of a temporal span or prosodic domain for the extraction of prosodic information is all the more crucial. While global measurements taken from entire recordings may provide a preliminary description on the prosodic characteristics of pathological speech, they cannot however capture the finer details of prosodic variations in speech. Yet, prosodic variations assume many communicative functions, e.g. give information about the hierarchical dimension and relational organization of discourse units and indicate a speaker's social states and attitudes [START_REF] De Looze | Analyse et interprétation de l'empan temporel des variations prosodiques en français et en anglais contemporain[END_REF]. In this study, we propose to use interpausal run (IR) as a unit for prosodic analysis as it can provide information on prosodic variations in speech and allows for an automatic analysis of large amounts of data. Acoustic analysis based on these interpausal run annotations specifically provides measurements of variations in pitch range, voice intensity and speech rate in order to further refine the current descriptions of Parkinsonian speech prosody. In this paper we first present an evaluation of three algorithms that we used for the automatic extraction of pitch extrema [START_REF] De Looze | Analyse et interprétation de l'empan temporel des variations prosodiques en français et en anglais contemporain[END_REF], syllable nuclei [START_REF] De | Praat script to detect nuclei and measure speech rate automatically[END_REF] and speech/pause intervals [START_REF] Boersma | Praat: doing phonetics by computer[END_REF] in order to investigate their reliability for the study of pathological speech. We then report our automatic analysis on the prosodic variation characteristics of Parkinsonian speech, which compares the read and semispontaneous speech of 25 control and 30 Parkinson subjects (off state).

Corpora

2.1

Selected subjects Thirty patients with Parkinson's disease (PD) (19 male and 11 female) were selected for this study. Speech recordings were carried out by expert phoneticians during a patient's visit to the Neurology Department of Aix-en-Provence Hospital, France. All patients met the UK Parkinson's disease Brain Bank Criteria for diagnosis of idiopathic PD. The average age of the patients at the time of recording was 68.1 ± 7.8 (standard deviation) years; mean PD duration was 8.8 ± 4.7 years. Patients were studied without (off), and then with (on) anti-Parkinsonian medication. For this pilot study, we selected the off medication condition only. For this condition, patients were recorded after a full overnight fast, i.e. at least 12 hours of PD treatment withdrawal. The patients' global motor disability was assessed using the motor section of the Unified Parkinson's Disease Rating Scale (UPDRS, part III, [START_REF] Fahn | Members of the UPDRS Development Committee. Unified Parkinson's Disease Rating Scale[END_REF]). Mean UPDRS in the off medication state was 30.1 ± 10.8. The specific item on speech disability (i.e. UPDRS item 18) was rated on average at 1.5. This item is a 5-point EAI scale (Equal Appearing Interval) where 0 = Normal speech; 1 = Slight loss of expression, diction and/or volume; 2 = Monotone, slurred but understandable; moderately impaired; 3 = Marked impairment, difficult to understand; 4 = Unintelligible. Our results (1.5 off) show that the PD speakers are globally slightly dysarthric. The speech of thirty-three control subjects (18 male and 15 female speakers) was recorded in the same manner. Mean age of the control group (C) at the recording time was 67.6 ± 10.8 years.

2.2

Recordings and task Within the Neurology department, a quiet room was set up with EVA2 computerised speech acquisition equipment (SQLab & LPL, Aix-en-Provence, France; www.sqlab.fr/). Speech was recorded with a Headset Cardioid Condenser Microphone with an adjustable headband (AKG C420). Speakers were first instructed to read a standardised text (i.e. "La chèvre de M. Seguin" by A. Daudet) at their most comfortable pitch and loudness, with natural effort. They were then asked to perform a semi-spontaneous task -describing a picture ("Vol de gateau"). All data are part of a larger dysarthric speech corpus collected within the Neurology department [START_REF] Ghio | How to manage sound, physiological and clinical data of 2500 dysphonic and dysarthric speakers[END_REF].

Experiment

3.1

Automatic prosodic annotations In this experiment, we first investigated whether the three Praat scripts that we used for the automatic analysis of speech prosody were suitable and stable enough for Parkinsonian (dysprosodic) speech.

-Pitch range adaptation script (De Looze [START_REF] De Looze | Analyse et interprétation de l'empan temporel des variations prosodiques en français et en anglais contemporain[END_REF]) for the automatic extraction of f0 extrema -Syllable nuclei detection script (De Jong et al. [START_REF] De | Praat script to detect nuclei and measure speech rate automatically[END_REF]) for the automatic detection of syllable nuclei and measurement of speech rate -To TextGrid (silences) Praat command [START_REF] Boersma | Praat: doing phonetics by computer[END_REF] for the automatic annotation of speech/pause intervals.

3.1.1Pitch range adaptation

Computation of f0 is generally based on default values, in particular in the case of extreme ones. For instance, when creating a Pitch Object in Praat, pitch floor and ceiling parameters are set to the interval [60 -600 Hz] for general context to [100 -500 Hz] for female or [75 -300 Hz] for male voices. Such large pitch ranges, however, generate pitch tracking errors with inappropriate extreme values. In order to avoid these errors, De Looze [START_REF] De Looze | Analyse et interprétation de l'empan temporel des variations prosodiques en français et en anglais contemporain[END_REF] proposed a multipass process to detect f0. In the first step, the extreme values are set to 60 and 600 Hz. In the second step, the extreme values for a speaker are adapted according to the first pass results: Pitch floor and pitch ceiling are respectively set to the values q15*0.83 (where 'q' stands for percentile) and q65*1.92. These formulae (as well as the formulae q25*0.75 -q75*1.5 and q35*0.72 -q65*1.90) have been shown to provide a better estimation of pitch extrema in 'normal' speech -as they help to exclude more octave errors at the extremes of the f0 distribution -than setting pitch floor and ceiling parameters to larger pitch ranges (e.g. [60 -600]; [100 -500] and [75 -300]). This pitch range adaptation script was therefore used to automatically extract pitch values in our experiment.

Syllable nuclei detection

Measuring articulation rate necessitates annotation in terms of syllables or phonemes. As mentioned in the introduction, this task can however be both time-consuming and error-prone in manual annotations. As a result, several automatic algorithms have been developed for phoneme and syllable detection. Most of them however (e.g. EasyAlign [START_REF] Goldman | EasyAligner: a semi-automatic phonetic alignment tool under Praat[END_REF]; LIA Aligner [START_REF] Audibert | Evaluation d'un alignement automatique sur la parole dysarthrique[END_REF]) are based on orthographic transcription of speech. In order to offer a fully automatic process, De Jong et al [START_REF] De | Praat script to detect nuclei and measure speech rate automatically[END_REF] developed an algorithm that detects syllable nuclei and measures speech rate automatically without requiring any transcription of the speech material as input. In this script, syllable nuclei correspond to peaks in intensity (dB) that are preceded and followed by dips in intensity, with unvoiced peaks being discarded. De Jong's algorithm has been shown to be suitable for the study of non-pathological Dutch speech. In this experiment, the script was used to automatically annotate syllable nuclei. The script's original output was modified to obtain a TextGrid with a tier annotated in syllables.

Silent pause detection

In order to annotate recordings into interpausal runs, the Praat command To TextGrid (silences) [START_REF] Boersma | Praat: doing phonetics by computer[END_REF] was used in this experiment. It provides an annotation of silent and sounding intervals, based on the sound's intensity contour. Several silence thresholds ranging from -20dB to -40dB were tested candidate borders between silence and sounding.

3.2

Automatic prosodic measurements The prosodic characteristics of Parkinson patients' and of control subjects' speech were investigated in terms of pitch range (level and span), voice intensity (relative level and span) and speech rate (articulation rate and pauses). Pitch level was measured by calculating the f0-median and given on a linear scale (i.e. Hertz); pitch span was measured by calculating the f0-standard deviation (sd-f0). The intensity of the voice was expressed in terms of the intensity curve median (Int-median) and standard deviation (Int-sd). Articulation rate was measured by calculating the number of syllable nuclei per second (syllsec). Silent pauses were measured in terms of number and mean pause duration. Number of interpausal runs (IR) and the amount of speech and pause for the recordings were also calculated to further flesh out the description of global speech/pause patterns. Prosodic features were extracted from two different temporal spans: (1) on the whole recording and (2) for each IR. Differences in f0-median (d_f0-median), sd-f0 (d_sd-f0), Intmedian (d_Int-median), Int-sd (d_Int-sd) and syllsec (d_syllsec) between consecutive IR were also calculated.

Results

4.1

Evaluation of automatic tools Pitch range adaptation, Syllable nuclei detection and Silence detection algorithms were first evaluated by comparing manual annotations with automatic annotations. Evaluations were carried out for the 25 control subjects (C) and the 30 Parkinson's Disease (PD) patients in the reading task.

4.1.1Evaluation of pitch range adaptation algorithm

Automatic and manual annotations of pitch extrema were compared using paired t-tests. Results show significant differences between manual and automatic detection of extrema, the mean of the difference being smaller for the pitch range adaptation method than others, i.e. either adjusted to the speaker's gender ([100-500]; [75-300]) or to default values ([60-600]) (cf. Table 1). In addition, ANOVA analyses reveal that differences are larger for PD patients than for C subjects (p=0.02 Results show that the proposed pitch range adaptation script [START_REF] Cheang | An acoustic investigation of Parkinsonian speech in linguistic and emotional contexts[END_REF] provides a better account of f0 measurements than other pitch range settings. However, even with this multi-pass process, differences in Hz are still large (especially for max).

Several f0 extractions are therefore needed to define a speaker's pitch range limits.

Evaluation of syllable nuclei detection

Automatic and manual annotations of syllable nuclei were compared by computing F-measures (F1), calculated as follows:

(1) F1 = 2* Pr*Re Pr+Re where Pr denotes the precision (ratio of hits to all hits and false alarms) and Re the recall (ratio of hits to all hits in the set). Hits correspond to syllable nuclei detection within 20ms of the manually annotated syllables while missed values when outside 20ms. F-measures were significantly smaller for PD patients than for C subjects (p= 0.00776), which means that syllable detection is less accurate for the PD population than it is for C subjects. However, as for both populations F-measures are relatively high (C=0.90; PD=0.86), the annotation in terms of syllable nuclei can be considered reliable.

Evaluation of silent detection

Automatic and manual annotations of speech/pause intervals were compared by computing F-measures and Hamming distances (Hd). F-measures were computed using equation (1), while Hd were found by calculating the distance between two binary vectors (respectively corresponding to manual and automatic annotations) with pauses coded as zeros and speech as ones. Among the different tested thresholds, results show that for our corpus a threshold set at -30dB provides higher Fmeasures and smaller hamming distances. In addition, results show that for this threshold, F-measures (F1) are high enough and Hd small enough to provide a reliable automatic annotation in terms of speech/pause intervals (F1=0.76; HD=0.10). In addition, no significant difference for Hd and Fmeasure were found between PD and C subjects.

Evaluation of extracted prosodic features

Automatic annotations were also evaluated by comparing the values of the prosodic features extracted from manual and automatic annotations using paired t-tests. Significance was set at p<0.01.

With respect to global measurements: T-tests based on global

measurements show that articulation rate and mean pause durations obtained with automatic annotation (AA) are not significantly different from those obtained with manual annotation (MA) for both population. Moreover, number of pauses and number of runs are not significantly different in AA and MA for PD patients. Results show however that mean run durations measured from AA are significantly different from those obtained from MA for both populations.

With respect to interpausal runs (IR) measurements: T-tests based on IR

show that prosodic values extracted from automatic annotations (AA) are not significantly different from those extracted from manual annotations (MA), except for Int-median and syllsec. Moreover, no significant differences are found between PD and C subjects for all parameters extracted. Prosodic values extracted from AA are thus as reliable for PD patients as they are for control subjects.

Prosodic characteristics of Parkinsonian speech

In order to investigate the prosodic characteristics of Parkinsonian read and semi-spontaneous speech, the data were analysed using the 'R' software using linear mixed-models, with Type (PD vs. Control) and Task (reading vs. semispontaneous) as fixed factors and Speaker as random factor. P-values were estimated using MCMC sampling and significance was considered for p<0.01. Tables 2 and3 provide the mean values obtained for each prosodic parameter according to Type and Task as well as the absolute differences of these mean values obtained between Parkinson (PD) and Control (C) subjects in read (R) and semi-spontaneous (S) task.

4.2.1Global measurements

Analyses based on global measurements revealed a highly significant main effect of Type on mean pause duration (mpauseD) and speech/pause ratio (sp_ratio). Results show that mpauseD and pause_ratio are higher for PD patients than they are for C subjects. No significant differences between C and PD subjects were found for f0-median, sd-f0, mean run duration, number of pauses and number of runs. Analyses also revealed a highly significant effect of Task on median_f0, mpauseD, pause_ratio and number of pauses. Results show that f0-median and npauses are higher in the R than in the S task; In addition, pause_ratio is smaller for R than for S. An interaction between Type and Task reached high significance for speakers' mean pause duration and speech/pause ratio. No interaction was found between Type and Task for the other prosodic parameters.

C/R C/S PD/R PD/S C-PD/R C-

Interpausal runs

Analyses based on interpausal runs (IR) revealed a highly significant effect of Type on differences between consecutive IR, in terms of f0-median, sd-f0 and syllsec. Results show that d_f0-median, d_sd-f0 and d_syllsec between runs are smaller for PD than they are for C. No significant differences were found between C and PD subjects in terms of IR f0-median, sd-f0 and syllsec, neither in terms of difference between consecutive IR in terms of Int-median and Int-sd.

Analyses also revealed a highly significant effect of Task on all prosodic parameters. IR f0-median and syllsec are shown to be higher in reading task (R) than in semi-spontaneous task (S); Differences in f0-median, sd-f0, syllsec, Int-median and Int-sd between IR are lower in R than in S. An interaction between Type and Task reached high significance for IR f0-median, sd-f0 and syllsec and for consecutive IR differences in f0-median, sd-f0, syllsec, Intmedian and Int-sd: f0-median is lower, sd-f0 is narrower, d_f0-median, d_sd-f0 and d_syllsec are smaller in S for PD patients than for C subjects; syllsec, d_Int-median and d_Intsd for S task are higher for PD patients than for C subjects.

Discussion and conclusion

In this paper, we have presented an automatic analysis of the prosodic variations in the speech of 25 control and 30 Parkinsonian subjects (in off (non-medicated) state) in reading and semi-spontaneous tasks. We first evaluated the reliability of existing algorithms for the automatic analysis of speech prosody and have shown that these tools are efficient for the analysis of Parkinsonian speech if a large amount of data is processed and if measurements rely on the extractions of prosodic features at different time intervals. Our evaluation of the Pitch range adaptation script [START_REF] De Looze | Analyse et interprétation de l'empan temporel des variations prosodiques en français et en anglais contemporain[END_REF] reveals it is all the more important for pitch extraction extrema, for which several candidates at the extreme of the f0-curve must be considered in order to define a speaker's pitch range limits.

In addition, our study reveals that Parkinsonian speech is characterised by longer pauses in both reading and spontaneous tasks, and by lower pitch level, narrower pitch range and faster articulation rate in spontaneous task, for which the difference in terms of pause length between control and PD subjects is all the more important. These findings corroborate those of the literature [START_REF] Canter | Speech characteristics of patients with Parkinson's disease: Intensity, pitch and duration[END_REF][START_REF] Monrad-Krohn | Dysprosody or altered melody of language[END_REF][START_REF] Darley | Motor speech disorders[END_REF][START_REF] Weismer | Acoustic description of dysarthric speech: Perception correlates and physiological inferences[END_REF][START_REF] Caekebeke | The interpretation of dysprosody in patients with Parkinson's disease[END_REF][START_REF] Ho Ak | Regulation of parkinsonian speech volume: the effect of interlocutor distance[END_REF][START_REF] Teston | L'évaluation objective de la prosodie in Les dysarthries[END_REF][START_REF] Duez | Syllable structure, syllable duration and final lengthening in parkinsonian French speech[END_REF][START_REF] Cheang | An acoustic investigation of Parkinsonian speech in linguistic and emotional contexts[END_REF] with the exception that we found no difference in terms of pitch range and articulation rate between control and PD speakers in reading task. This could be explained by the fact that PD patients in our experiment are slightly dysarthric, their dysprosodic characteristics being only revealed in semi-spontaneous speech. These findings highlight the importance of evaluating PD patients in different tasks.

Our study also reveals that differences in terms of pitch range and articulation rate between consecutive interpausal runs are smaller in PD patients than in control subjects in both reading and spontaneous speech, these differences between IR being all the more smaller in semi-spontaneous speech. Our results further show that differences in voice intensity level between IR are smaller and in voice intensity range larger in PD speech than in control subjects' in semi-spontaneous task. These findings come to augment the current descriptions of Parkinsonian speech, where the reported monotonous voice is shown to be related not only to a narrowing of the pitch range but also to smaller variations in pitch range. Similarly, our results suggest that the reported monoloudness in PD speech could be due to smaller variations in voice intensity level. One may assume that because prosodic variations are smaller in PD speech, some communicative functions (e.g. topic changes, speakers' attitudes) may not be conveyed in read and in semispontaneous Parkinsonian speech.

Our study thus demonstrates that differences in prosodic features between consecutive interpausal runs better capture the prosodic characteristics of Parkinsonian speech than global measurements or interpausal runs' prosodic features. In the literature, most prosodic measurements have been made from entire recordings or from smaller temporal domains.

Investigating prosodic variations by calculating the differences between interpausal runs could therefore augment these established measurements, in particular when dealing with slightly dysarthric speech.

Table 1 .

 1 Mean of the differences (expressed in Hz) obtained for

	).
	f0-extrema	Mean of the diff	p-value
	minSR	9	3.389e-07
	minG	25	1.608e-12
	minD	33	5.939e-14
	maxSR	-44	6.882e-13
	maxG	-82	3.171e-11
	maxD	-243	< 2.2e-16
	the minima and maxima extracted with pitch ceiling and floors
	adjusted to default values (minD; maxD), adjusted according

to the speaker's gender (minG; maxG) and adjusted to the speaker's pitch range (minSR; maxSR).

Acknowledgements: This research was supported by ANR BLAN08-0125 of the French National Research Agency.