Guillaume Chèze
email: guillaume.cheze@math.univ-toulouse.fr

André Galligo
email: galligo@math.unice.fr

Bernard Mourrain
email: bernard.mourrain@inria.fr

Jean-Claude Yakoubsohn
email: jean-claude.yakoubsohn@math.univ-toulouse.fr

A Subdivision Method for Computing Nearest Gcd with Certification

A new subdivision method for computing the nearest univariate gcd is described and analyzed. It is based on an exclusion test and an inclusion test. The exclusion test in a cell exploits Taylor expansion of the polynomial at the center of the cell. The inclusion test uses Smale's α-theorems to certify the existence and unicity of a solution in a cell.

Under the condition of simple roots for the distance minimization problem, we analyze the complexity of the algorithm in terms of a condition number, which is the inverse of the distance to the set of degenerate systems.

We report on some experimentation on representative examples to illustrate the behavior of the algorithm.

Introduction

Computing an approximate gcd of two polynomials is a fundamental and difficult problem in symbolic-numeric computation. It has been extensively studied in the Computer Algebra community, following different strategies: [START_REF] Brown | On Euclid's algorithm and the computation of polynomial greatest common divisors[END_REF][START_REF] Schönhage | Quasi-gcd computations[END_REF][START_REF] Noda | Approximate GCD and its application to illconditioned algebraic equations[END_REF][START_REF] Corless | The singular value decomposition for polynomial systems[END_REF][START_REF] Hribernig | Detection and validation of clusters of polynomial zeros[END_REF][START_REF] Emiris | Certified approximate univariate GCDs[END_REF][START_REF] Beckermann | When are two numerical polynomials relatively prime?[END_REF][START_REF] Emiris | Numerical Univariate Polynomial GCD[END_REF][START_REF] Bini | Polynomial and matrix computations[END_REF][START_REF] Bini | Structured matrix-based methods for polynomial -gcd: analysis and comparisons[END_REF][START_REF] Victor | Computation of approximate polynomial GCDs and an extension[END_REF][START_REF] Zeng | The approximate GCD of inexact polynomials. Part I: a univariate algorithm[END_REF][START_REF] Karmarkar | On approximate GCDs of univariate polynomials[END_REF][START_REF] Kaltofen | Structured low rank approximation of a Sylvester matrix[END_REF][START_REF] Nie | Global minimization of rational functions and the nearest GCDs[END_REF][START_REF] Terui | An Iterative Method for Calculating Approximate GCD of Univariate Polynomials[END_REF]. . . . The problem is usually formulated as finding a perturbation of a pair of polynomials (f, g) within a given ball of radius ε, such that the perturbed pair (p, q) has a non-trivial exact gcd (of highest possible degree) in this ball.

In our approach, we reformulate the problem as an optimization problem: Given the polynomial pair (f, g) of a given degree, we look at the nearest polynomial pair (p, q) of the same degree which has a non-trivial exact gcd. The ε-approximate gcd problem has a solution iff this nearest pair with a non-trivial gcd is within a distance ε.

This paper is an extended version of the conference paper [START_REF] Chèze | Computing Nearest Gcd with Certification[END_REF]. In the first section, we present the optimization problem that we are going to solve, we state our complexity main result and mention the related approaches in the literature for approximate gcd computation. In section 2, we describe the main ingredients of our subdivision method. In section 3, we analyze its complexity using α-theory. Finally, in section 4, we report of some experimentation to illustrate the behavior of the algorithm.

A minimization problem

As usual, we denote by R d [z] and C d [z] the vector spaces of univariate polynomials of degree less or equal to d, with coefficients in R and C. We denote by f = (f 0 , . . . , f d ) T the coefficients of the polynomial f (z) =

d k=0 f k z k ∈ R d [z]. The 2-norm of f is ||f || 2 = d k=0 f k fk .
We address the nearest gcd problem, i.e. the following minimization problem: Problem 1.1 Given two degree d polynomials f (z) and g(z), find two degree d polynomials p(z) and q(z) with a non trivial gcd, which are solutions of the minimization problem

min p,q∈C d [z] Resultant d (p,q)=0 ||f -p|| 2 + ||g -q|| 2 , ( 1
)
Here Resultant d (p, q) is the resultant of p and q when they are considered as polynomials of degree d. We recall that Resultant d (p, q) = 0 if and only if gcd(p, q) is non-trivial or deg (p) < d and deg (q) < d. In the last case p and q have a common root at infinity, see Example 2 in Section 4.

In our formulation of the nearest gcd problem, we look for the nearest pair (p, q) of degree at most d with a non-trivial gcd. In Example 1 in Section 4 we will see that if deg f = 7 and deg g = 8 then we can find p and q with degree 8. Thus the degree of p can be bigger than the degree of f . N.K. Karmarkar and Y.N. Lakshman in [START_REF] Karmarkar | On approximate GCDs of univariate polynomials[END_REF] reduced that problem to another minimization problem :

Problem 1.2 min z∈P(C) f (z)f (z) + g(z)g(z) d k=0 z k z k . ( 2
)
Moreover, a global minimum z 0 of (1.2) is a root of a nearest gcd and determines the polynomials p and q of (1.1) by the formulas:

p(z) = f (z) - f (z 0 ) d k=0 z 0 k z k 0 d k=0 z 0 k z k q(z) = g(z) - g(z 0 ) d k=0 z 0 k z k 0 d k=0 z 0 k z k .
We observe that with the change of variables z = X + iY , the problem amounts to minimize a homogeneous rational bivariate function

F(z) = F (X, Y ) = N (X, Y ) D(X, Y ) .
So we focus on the resolution of this task. 2

Notations

We denote by G(X, Y ) = (G 1 (X, Y ), G 2 (X, Y )) the couple of numerators of the gradient ∇F (X, Y ). We denote by P µ = N -µD where µ ∈ R and F = N D .

For

G(X, Y ) = (G 1 (X, Y ), G 2 (X, Y )) with G k (X, Y ) = deg G k l=0 i+j=l G k;i,j X i Y j ∈ C[X, Y ], k = 1, 2,
we define the (Bombieri) norm G B by:

G 2 B = k=1,2 deg G k l=0 i+j=l G k;i,j G k;i,j i!j! deg G k ! . A F := {(x i , y i ) | i = 1, .., n } is the set of global minima of F (X, Y ). We assume that this set is finite. Z H is the set of zeroes of a polynomial function H(X, Y ) and d(x, y, Z H ) is the Euclidean distance from (x, y) ∈ R 2 to Z H . The function H(X, Y ) will be specialized to P µ (X, Y ) or to G i (X, Y ), i = 1, 2.
Let S or S(x, y, r) be a square centered at (x, y) and of radius r.

Our approach

We use a bisection algorithm, based on an exclusion test and an inclusion test (see below), applied simultaneously to:

1-the polynomial P µ (X, Y ) = N (X, Y ) -µD(X, Y ) 2-the system G(X, Y ) = (G 1 (X, Y ), G 2 (X, Y )).
Our bisection algorithm will iteratively update a list of retained squares and a list of approximate global minima of F . We prove a complexity result in Section 3. We use G 1 and G 2 to compute local extrema. The polynomial P µ is used to exclude "bad" squares. Indeed, during our algorithm µ is a candidate for the minimum of F and we set µ = N (x 0 , y 0 )/D(x 0 , y 0 ) where (x 0 , y 0 ) is a center of a square. Thus, if P µ (x, y) > 0 in a square S then F (x, y) > µ, and we cannot find a global minimum in S.

We follow the approach initiated by S. Smale and his co-worker in a series of papers (see e.g. [5] and the references therein) relying on their celebrated α and γ-theorems.

In section 2 we will briefly review some of these notions and recall the definitions of α, β, γ, we also set γ(G, A F ) := max (x,y)∈A F γ(G; x, y). After that, we provide a precise quantitative definition for a point (x, y) to be, in our setting, an approximate global minimum of F (x, y). Then thanks to the γ-theorem, the approximation is sufficiently good to imply that a Newton iteration converges quadratically towards a global minimum.

To state our result on the complexity analysis of the subdivision method, we need the following notations: ) where d F is the Frobenius distance and Σ is the set of singular matrices.

• σ F := 1 min (x,y)∈A F d F (DG(x,y),Σ
• V (G) is the tubular neighborhood of the zero set G(x, y) = 0;

• N is the maximal number of connected components of U := V (G 1 ) ∩ V (G 2 ), for all > 0;

• for all > 0, K is bounding the radius of a connected component of U .

• J (G, F, r 0 ) := log 2 ( K r0γ(G,A F ) δ 0

) where r 0 > 0 and δ 0 is the smallest positive root of

(13 -3 √ 17) (2u 2 -4u + 1) 2 -4u = 0.
Theorem 1.3 (Global minimum) Assume that all the elements of A F are regular points of the system G(x, y) = 0 and are contained in some square S 0 := S(x 0 , y 0 , r 0 ). Denote by R j the set of retained squares at step j. Then, for j ≥ J = J (G, F, r 0 ), all the points of R j are approximate global minima of F (x, y). That is to say the Newton iteration applied to G from a point in R j converges quadratically to a global minimum of F .

Theorem 1.4 (Complexity analysis) For the same notations and assumptions as in above theorem, we have the following properties:

1-The number of exclusion tests is bounded by 1 + J N K 2 . Remark:

2-When

1. σ F = 1 min (x,y)∈A F d F (DG(x,y),Σ
) is related to the condition number of the system. That is why we keep this constant in the big O notation.

Unfortunately, we do not know how to bound this constant in terms of the degree of f and g. The computation of this kind of bound is a difficult problem.

1. A F ≤ (p,q) solution of Problem 1.1 deg gcd(p, q).

Thus if Problem 1.1 has a unique solution then

A F ≤ deg gcd(p, q) ≤ d.
2. In this paper, in order to apply α and γ-theorems, we suppose that all the elements of A F are regular points of the system G(x, y) = 0. We could avoid this hypothesis, instead of using α and γ-theorems, by relying on the tools developed in [START_REF] Giusti | On location and approximation of clusters of zeros of analytic functions[END_REF].

Other approaches in the literature

The nearest GCD problem has been studied with different approaches and with other formulations by many authors.

Algebraic approach via Euclid's algorithm, resultant and subresultant

The first papers, see for example Brown [START_REF] Brown | On Euclid's algorithm and the computation of polynomial greatest common divisors[END_REF], on the complexity of Euclid's algorithm only works with exact coefficients. Then the numerical case has been successively considered by Schönhage [START_REF] Schönhage | Quasi-gcd computations[END_REF], Noda and Sasaki [START_REF] Noda | Approximate GCD and its application to illconditioned algebraic equations[END_REF], Corless, Gianni, Trager and Watt [START_REF] Corless | The singular value decomposition for polynomial systems[END_REF], Hribernig and Stetter [START_REF] Hribernig | Detection and validation of clusters of polynomial zeros[END_REF], Emiris, Galligo and Lombardi [START_REF] Emiris | Certified approximate univariate GCDs[END_REF], Beckermann, Labahn [START_REF] Beckermann | When are two numerical polynomials relatively prime?[END_REF]. These authors consider the so-called near GCD or the -GCD problem. The singular value decomposition was first applied to the Sylvester matrix, in [START_REF] Corless | The singular value decomposition for polynomial systems[END_REF], and later applied to the subresultant matrices, in [START_REF] Emiris | Numerical Univariate Polynomial GCD[END_REF][START_REF] Emiris | Certified approximate univariate GCDs[END_REF] in order to get a certification when a condition (depending on the level of accuracy) is satisfied. An efficient implementation is done in [START_REF] Rupprecht | Eléments de Géométrie Algébrique Approchée: étude du PGCD et de la factorisation[END_REF].

Padé approximation and structured matrices approach

See the book of Bini and Pan [START_REF] Bini | Polynomial and matrix computations[END_REF] and the bibliography therein, and more recently Boito, Bini [START_REF] Bini | Structured matrix-based methods for polynomial -gcd: analysis and comparisons[END_REF].

Rootfinding and cluster root approach

Pan [START_REF] Victor | Computation of approximate polynomial GCDs and an extension[END_REF] and, more recently, Zeng [START_REF] Zeng | The approximate GCD of inexact polynomials. Part I: a univariate algorithm[END_REF] use root finding and least squares methods.

Optimization approach

The resolution of problem 1.2 is the main and the most "time-consuming" step of the method propose in [START_REF] Karmarkar | On approximate GCDs of univariate polynomials[END_REF] that we aim to improve. The authors rely on techniques from Arnon-McCallum [START_REF] Arnon | A polynomial-time algorithm for the topological type of a real algebraic curve[END_REF] and Manocha-Demmel [START_REF] Manocha | Algorithms for intersecting parametric and algebraic curves: I simple intersections[END_REF]. The second paper is based on resultant for expressing the intersection of two curves and on numerical computation of eigenvalues and eigenvectors by QR iterations. Therefore the expected running time is in O(p 3 ) where p is the product of degrees of two curves, hence the complexity of the algorithm is at least in O(d 6 ).

Kaltofen and his co-workers [START_REF] Kaltofen | Structured low rank approximation of a Sylvester matrix[END_REF] determine approximate GCDs from methods based on structured total least square (STLN). The STLN is an iterative method of the family of Gauss-Newton methods. The authors describe its application to the case of the Sylvester matrix associated to the input polynomials f and g, then show its interest and efficiency by producing the results of experiments. However, since the starting point of their Gauss-Newton like method is not precised, the method can diverge. This is an important drawback.

Nie-Demmel-Gu [START_REF] Nie | Global minimization of rational functions and the nearest GCDs[END_REF] uses a sum of squares (SOS) technique, the resolution relies on semi definite programming (SDP). For the nearest GCD problem, this yields linear matrix inequalities (LMI) whose size s is O((d + 2)(d + 1)) whose complexity of resolution by projective algorithm is in O(s 3 ) [START_REF] Nesterov | Interior-point polynomial algorithms in convex programming[END_REF]. So, one ends up with a complexity in O(d 6 ).

More recently, Terui [START_REF] Terui | An Iterative Method for Calculating Approximate GCD of Univariate Polynomials[END_REF] uses a generalization of the gradient-projection method proposed by Tanabe [START_REF] Tanabe | A geometric method in nonlinear programming[END_REF]. Then he obtains a fast algorithm for solving a constrained minimization problem related to the approximate gcd problem. Unfortunately, there exists no complexity results for this approach. Furthermore, in this method the degree of the approximate gcd is given in input.

The proposed bisection method

Principles

We propose a bisection method, described below, to approximate the global minima of the function F (X, Y ) defined in the introduction. In the next section, a procedure for computing an initial square is given ; so we suppose here that a square S(x 0 , y 0 , r 0 ) containing all the global minima is known.

The bisection method is based on an exclusion test and an inclusion test. An exclusion test, denoted hereafter by E(F, S) or E + (F, S), is defined on the set of squares and returns true if the function F has no zero in the square S or false if it might have a zero. Examples of exclusion tests are provided in section 3.

An inclusion test I(G, S) is needed to numerically prove the existence of a local minimum: it takes G = (G 1 , G 2 ) and a square S then returns true if there exists a ball B(x * , y * , r * ), containing S, which contains one zero of G, in that case it also returns (x * , y * , r * ), otherwise it returns false.

Our definition of approximate minima is based on γ-theorem and α-theorem of Smale [START_REF] Shub | Complexity of Bézout's theorem. I. Geometric aspects[END_REF], [5], [START_REF] Giusti | On location and approximation of clusters of zeros of analytic functions[END_REF], applied to the system G(x, y). We first introduce some quantities and the corresponding classical notations:

• β := β(G; x, y) = ||DG(x, y) -1 G(x, y)|| • γ := γ(G; x, y) = sup k≥2 `1 k! ||DG(x, y) -1 D k G(x, y)|| ´1/(k-1) • α := α(G; x, y) = βγ. • γ(G; A F ) = max (x,y)∈A F γ(G; x, y).
Smale's γ-theorem [5] states:

Theorem 2.1 Let (x * , y * ) be a zero of G and suppose that DG(x * , y * ) is invertible. If (x, y) -(x * , y * ) γ(G, x * , y * ) ≤ 3 - √ 7 2
then the Newton iteration from (x, y) converges quadratically to (x * , y * ).

This leads to the following definition and quantitative results on the convergence of Newton scheme.

Definition 2.2 The point (x, y) is an approximate global minimum of F (x, y) if d((x, y), A F ) ≤ 3 - √ 7 2γ(G, A F ) .
Under this condition, the γ-theorem asserts that the Newton iteration from any point in the ball B(x * , y * , 3-

√ 7 2γ(G,A F )
) where (x * , y * ) ∈ A F , converges quadratically towards (x * , y * ). The following result gives a sufficient condition for a point to be a good starting point for the Newton iteration. This result is called α-theorem.

σ(x, y) = 1 + α - √ 1 -6α + α 2 4γ ≤ 2 - √ 2 2γ .
and the Newton iteration from (x, y) converges to (x * , y * ). Furthermore, we have

(x * , y * ) -N k G (x, y) ≤ 5 - √ 17 4γ 1 2 2 k -1
, where N k G (x, y) is the k-th iterate of the Newton iteration applied to G with starting point (x, y). The following result gives a radius r * of a ball centered in (x * , y * ) such that every point in B(x * , y * , r * ) satisfies the hypothesis of the α-theorem. Remark: 0.07 ≤ δ 0 ≤ 0.08.

Proof. We use the following inequality, see [START_REF] Shub | Complexity of Bézout's theorem. I. Geometric aspects[END_REF] α(G, x, y) ≤ u Ψ(u) 2 , where u = (x, y) -(x * , y * ) γ(G, x * , y * ). Furthermore, we have u/Ψ(u) 2 ≤ (13 -3 √ 17)/4 for u ∈ [0; δ 0 ] and this gives the desired result.

In our complexity study we will need a bound on γ(G, A F ). The following proposition will be useful.

Proposition 2.5 Let deg G = max(deg G 1 , deg G 2 ), (x, y) ∈ S(0, 0, r 0 ), Then γ(G; x, y) ≤ G B DG(x, y) -1 (deg G) 2 (1 + 2r 2 0 ) (deg G-2)/2 .
Proof. We set some notations: (x, y) 1 = (1 + x 2 + y 2 ) 1/2 , ∆(a i ) is the diagonal matrix with coefficients a i . By [23, Proposition 3], we have the following bound:

γ(G; x, y) ≤ C(G; x, y)deg G 3/2 2 (x, y) 1 , where C(G; x, y) = G B DG(x, y) -1 ∆(deg (G i ) 1/2 (x, y) deg Gi-1 1
) . We use the following classical inequality AB ≤ A B to conclude.

Let a threshold > 0 be given, the output of this bisection method will be a set (eventually empty)

Z = {(x * i , y * i , r * i , µ * i )}, such that : 1-the (x * i , y * i )'s are approximate global minima of F (x, y) and µ * i = F (x * i , y * i ).

2-the ball B(x *

i , y * i , r * i ) contains one and only one zero of G(x, y) where

r * i = σ(x * i , y * i ). 3-If i = j then B(x * j , y * j , r * j ) ∩ B(x * i , y * i , r * i ) = ∅ and |µ * i -µ * j | < .

Computation of an initial square

Lemma 1 in [START_REF] Karmarkar | On approximate GCDs of univariate polynomials[END_REF] gives us a bound for the initial square:

Lemma 2.6 Let f = d 1 k=0 f k z k , g = d 2
k=0 g k z k , where f d 1 = 0 and g d 2 = 0. Let F be the rational bivariate function corresponding to the approximate gcd problem. Let (x, y) ∈ A F then

(x, y) ≤ 5 max f 2 f 2 d1 , g 2 g 2 d2 .

Sketch of algorithm

The algorithm consists of an initialization followed by a while loop with an internal for loop. We call step k, the k-th step of the while loop. Proposition 2.7 Assume there exists > 0 such that for all square S(x, y, ) ⊂ S 0 we have : 1-the inclusion test is true if S(x, y, ) contains a zero of G.

2-the exclusion test is true if S(x, y, ) does not contain a zero of G.

Then we have :

Algorithm 2.1: Approximate gcd Input: F = N/D, G = (G 1 , G 2 )
, an initial square S 0 and a threshold > 0, as described above. Create a set of squares L := {S 0 }, a set of solutions Z := ∞ and a value (the minimum to be updated) µ := F(∞).

While L is not empty do

• Compute δ = minF (x i , y i ) where the (x i , y i )'s are the centers of squares of L.

• If δ < µ then µ := δ.

• For each square S of L perform the exclusion tests E + (P µ , S), E(G 1 , S) and E(G 2 , S).

If at least one of these exclusion tests is true then remove S from L ; else, perform an inclusion test I(G, S).

• If it returns false then divide S in 4 equal squares ; else an approximate local minimum (x * , y * ) is provided, it is the unique zero of G(x, y) in the ball B(x * , y * , r * ).

-If µ * = F (x * , y * ) > µ + then remove S ; else update the set Z as follows:

* If µ * = F (x * , y * ) < µ then put µ := µ * . For each element (x * i , y * i , r * i , µ * i ) of Z do • If µ < µ * i and |µ -µ * i | > then remove the solution (x * i , y * i , r * i , µ * i ) from Z ; • If for any element (x * i , y * i , r * i , µ * i ) of Z, |µ -µ * i | < and (x * , y * ) is not in the ball B(x * i , y * i , r * i ) then add (x * , y * , r * , µ * ) to Z. Output:
The set Z of approximate global minima of F .

1-The algorithm stops.

2-Let µ( ) be the value of µ at the end of the algorithm with input . Then lim →0 µ( ) = min (x,y)∈S 0 F (x, y).

Proof. The point 1 holds by construction under these assumptions. For point 2, the radius r( ) of the retained squares in the last iteration decreases. Hence from the continuity of F , we obtain lim →0 µ( ) = min (x,y)∈S0 F (x, y).

Complexity analysis

Exclusion test

Let H be a polynomial in R[X, Y ] and denote by D k H(x, y) the homogeneous part of degree k of the Taylor expansion of H at the point (x, y).

Let S(x 0 , y 0 , r 0 ) be a square. To define an exclusion function E(H, S), we rely on the following expression and lemma:

M H (x 0 , y 0 , r 0 ) = |H(x 0 , y 0 )| - k≥1 ||D k H(x 0 , y 0 )|| k! r k 0 .
Lemma 3.1 If M H (x 0 , y 0 , r 0 ) > 0 then the closed square S(x 0 , y 0 , r 0 ) does not contain any zero of the polynomial H(x, y).

Proof. The proof follows from Taylor formula and a simple inequality.

A key to analyze the complexity of the algorithm of section 2 is the following lemma, see [START_REF] Dedieu | Computing the distance from a point to an algebraic hypersurface[END_REF].

Lemma 3.2 Let H ∈ R[X, Y
] be a polynomial of degree e, consider the associated algebraic variety

Z H = {(x, y) ∈ R 2 : H(x, y) = 0}.
Let L e = 2 1/e -1 √ 2 and m H (x, y) be the function implicitly defined by

M H (x, y, m H (x, y)) = 0.
Then m H (x, y) is related to the distance d(x, y, Z H ) by the following inequalities:

L e .d(x, y, Z H ) ≤ m H (x, y) ≤ d(x, y, Z H ).
The exclusion tests to be used for the algorithm of section 2 are defined for a polynomial P by:

1-E(P, S) is true if M P (x, y, r) > 0,

2-E + (P, S) is true if M P (x, y, r) > 0 and P (x, y) > 0.

Since the degree of P µ is 2d and the degree of the G i 's is 4d -2, we get: if

P µ (x, y) > 0 then m P µ (x, y) ≥ L 2d .d(x, y, Z P µ ) and m G i (x, y) ≥ L 4d-2 .d(x, y, Z G i ). Remark: deg G i ≤ 4d -2 because the coefficient of the term of degree 4d -1 is 2df 2 d + 2dg 2 d - 2df 2 d -2dg 2 d = 0.
Putting these facts together, we obtain the following proposition:

Proposition 3.3 Let H(X, Y ) = P µ (X, Y ) or H(X, Y ) = G i (X, Y ), i = 1, 2. Let S := S(x, y,

r). The logical relation E(H, S

) is true means m H (x, y) > r, so S does not contain any zero and is excluded in the bisection algorithm. Otherwise

E(H, S) is false means m H (x, y) ≤ r,
and S may contain zeros. In this case S will be divided into four squares each of them with a radius r/2.

Inclusion test

Our inclusion test is based on Smale's α-theory. The test is true if :

1-α(G; x, y) < 13 -3 √ 17 4
, and 2-S(x, y, r) ⊂ B x, y, σ(x, y) , in other words r √ 2 < σ(x, y), see Theorem 2.3.

Proof of theorem 1.3

Consider a retained square S := S(x l , y l , r 0 2 k ) at step k. We have µ := µ k and

P µ (x l , y l ) ≥ 0. If E(P µ , S) is false and E(G i , S) is false, it means, by Proposition 3.3, that P µ (x l , y l ) = 0 or m Pµ (x l , y l ) ≤ r 0 2 k and m Gi (x l , y l ) ≤ r 0 2 k . From Lemma 3.2, we deduce that if P µ (x l , y l ) > 0 then L 2d .d(x l , y l , Z P µ ) ≤ m(x l , y l ) ≤ r 0 2 k .
Hence if P µ (x l , y l ) > 0, then for each (x, y) ∈ S we have

d(x, y, Z Pµ ) ≤ ||(x, y) -(x l , y l )|| + d(x l , y l , Z P µ ) ≤ σ k r k
where σ k = (1 + 1/L 2d ) and r k := r 0 2 k . In the same way,

d(x, y, Z G i ) ≤ ||(x, y) -(x l , y l )|| + d(x l , y l , Z G i ) ≤ ω k r k , i = 1, 2.
where

ω k = (1 + 1/L 4d-2 ). For H(x, y) ∈ R[x, y] and > 0, let V (H) = {(x, y) ∈ R 2 , st. ∃(u, v) ∈ C 2 , H(u, v) = 0, d((x, y), (u, v)) < }.
V (H) is called the tubular neighborhood of the solution set of H(x, y) = 0 at distance ε. By the previous inequalities, if P µ (x l , y l ) > 0 then we have

S ⊂ U k := V σ k (P µ ) ∩ V ω k (G 1 ) ∩ V ω k (G 2 ). Notice that A F ⊂ ∩ k≥0 U k . Let κ k r k be half the maximal diameter of a connected component U k . The area of U k is bounded by π ν k κ 2 k r 2
k where ν k is the number of connected components of U k . For k big enough, this number of connected components is the number of real roots of G 1 (x, y) = G 2 (x, y) = 0 with P µ (x, y) > 0. We denote by N < ∞ the maximum of all the ν k .

When k tends to infinity, the connected components of U k tend to the real roots of G 1 (x, y) = G 2 (x, y) = 0 with P µ (x, y) > 0 and κ k tends to 2 (1 + cos(α)) where α is the angle between the tangents of the curves G 1 (x, y) = 0 and G 2 (x, y) = 0, at the root. Let K bound all the possible κ k .

Let J be the first J such that K r J ≤ δ 0 γ(G, A F ) , then for k ≥ J, all the points of the retained squares are approximated zeros of the set A F , by Proposition 2.4 and Theorem 2.3. An upper bound for J is given in the theorem.

Proof of theorem 1.4

To compute an upper bound for the number of exclusion tests, we notice that at step k, U k contains the union of the kept squares (whose exclusion test is false), their number is denoted by q k . Since the area of this union is q k r 2 k and must be less or equal to the area of U k , we get

q k ≤ π ν k δ 2 k ≤ π N K 2 .
Now, let p k be the number of excluded squares at step k. As we know the relation 4q k-1 = p k +q k with p 0 = 0, q 0 = 1 holds, the number of exclusion tests until step j is bounded by

j k=0 p k + q k = 1 + 4 j-1 k=0 q k ≤ 1 + j π N K 2 .
The last part of Theorem 1.4 comes from Proposition 2.5 which gives: for (x, y)

∈ A F , log γ(G; x, y) belongs to O log d 2 r 4(d-1) 0 ||G|| B ||DG -1 (x, y)|| or to O log(d) + d log(r 0 ) + log(||G|| B ) + log(σ F ) .

Examples

Our algorithm have been implemented in Matlab 7 (we use an Intel Xeon 3 Ghz processor) and we have perform some experiments to show its behavior. Currently, our algorithm is slower than other methods (e.g. [START_REF] Bini | Structured matrix-based methods for polynomial -gcd: analysis and comparisons[END_REF][START_REF] Kaltofen | Structured low rank approximation of a Sylvester matrix[END_REF][START_REF] Terui | An Iterative Method for Calculating Approximate GCD of Univariate Polynomials[END_REF]) due to the bound on the initial square given in Lemma 2.6. Nevertheless, our algorithm is certified. In the following we give examples and timings.

Increasing the degree

We consider the following example taken from the work of D. Rupprecht [START_REF] Rupprecht | Eléments de Géométrie Algébrique Approchée: étude du PGCD et de la factorisation[END_REF]. In his Ph.D. thesis, he developed a technique which allowed him to certify the degree of an approximate gcd but only if the required precision ε belongs to some intervals. There are small gaps between these intervals and the numerical computation of an approximate gcd is sensitive.

In our formulation of the nearest gcd problem, we look for the nearest pair (p, q) of degree d with a non-trivial gcd. The order of the perturbation is 10 -8 and the degree of the corresponding gcd is 1. Note that in this example, even if f is of degree 7, we allow a perturbation of degree 8 (which is the degree of g).

Solution at infinity

Here, we show that the global minimum can be reached at infinity. We consider:

f = z 2 -4z + 3, g = z 2 + 4z + 3.
In this situation, we have:

F (X, Y ) = N (X, Y ) D(X, Y ) = 2 9 + 22X 2 + 10Y 2 + (X 2 + Y 2 ) 2 1 + X 2 + Y 2 + (X 2 + Y 2 ) 2 .
Figure 1 shows the graph of F . Notice that there are no global minima of F in C. Indeed, we have lim

(x,y) →∞ F (x, y) = F(∞) = 2, N (X, Y ) -2D(X, Y ) = 16 + 42X 2 + 18Y 2 > 0.
Then the global minimum is reached at infinity. In this situation the initial square given by Lemma 2.6 is [-130, 130] 2 , and our algorithm excludes this square after 10 iterations. That is to say, the size of the smallest excluded square is 130/2 10 . The cpu-time in this example is 0.06 seconds.

In this situation, we have F(∞) = 2. Furthermore, with z 0 = ∞ and the formula recalled in Subsection 1.1, we obtain p = -4z+3 and q = 4z+3. Notice that deg (p) < 2 and deg (q) < 2.

Remark: We compute F(∞) in the following way: F(∞) = F r (0), where F r (z) = z deg (F ) F(1/z).

Example 3

Here, we consider:

f = 9z 3 -18z 2 + 3z + 1, g = -37z 3 + 63z 2 -30z + 3.
We illustrate the algorithm step by step. The smallest perturbation that we obtain is: p(z) = 9.059z 3 -17.875z 2 + 3.26z + 1.544, q(z) = -36.916z 3 + 63.174z 2 -29.636z + 3.758, and the common root is z = 0.479.

Step 1. We begin in Figure 2 (a) with a big square and we divide it in four squares.

( Step 2. As we cannot exclude any squares and all the inclusion tests are false we divide each square in four squares (Figure 2 (b)).

During the algorithm we compute F for each center of each square. The point in the figure corresponds to the center where F is minimum. Furthermore, the value of F at this point is equal to the real number µ defined in our algorithm.

Step 3. In Figure 2 (c), we cannot exclude any squares and the inclusion tests are false.

Step 4. At this step, see Figure 2 (d), we can exclude a lot of squares.

Step 5. In Figure 2 (e), we see that we can exclude a lot of squares and one inclusion test is true. So with the Newton's method applied to G we can compute the global minima. That is to say we apply Newton's method with the black point as a starting point and we get the red point (on the x-axis) as the global minimum.

Step 6. In this last step, can exclude last squares and we find the global minimum, see Figure 2 (f).

Random examples

We have carried out the following tests:

• f and g are polynomials of degree d with random coefficients. Our results are given in Figure 3, and cpu denotes the cpu-time needed in seconds to performs our algorithm.

• f (z) = (z -1)f 1 (z) + 10 -3 1 (z), g(z) = (z -1)g 1 (z) + 10 -3 2 (z), where f 1 and f 2 are random polynomials with degree d -1, 1 and 2 are random polynomials with degree d. In these examples our algorithm gives an approximate gcd near z -1. Our results are given in Figure 4, and cpu denotes the cpu-time needed in seconds to performs our algorithm.

Here random means that we use the Gaussian distribution with mean 0 and variance 1.

Theorem 2 . 3

 23 If α < (13 -3 √ 17)/4 then G has one and only one zero in the open ball B(x, y, σ(x, y)) with

Proposition 2 . 4 3 √ 17 4 Ψ

 2434 Let (x * , y * ) be a zero of G and suppose that DG(x * , y * ) is invertible. Let δ 0 be the smallest positive root of u -13 -(u) 2 , where Ψ(u) = 2u 2 -4u + 1. If (x, y) -(x * , y * ) γ(G, x * , y * ) ≤ δ 0

Figure 1 :

 1 Figure 1: F has no global minima in C.

Figure 2 :

 2 Figure 2: Subdivision steps

Figure 4 :

 4 Figure 3: Random examples

  The nearest perturbed pair (p, q) with a non-trivial gcd that we compute is:

			p	=	(0.000000001537219 + 0.000000006148875 i) z	8
						+(0.999999993851125 + 0.000000001537219 i) z	7
						+(1.999999998462781 -0.000000006148875 i) z	6
						+(-3.000998983851125 -0.000000001537219 i) z	5
						+(0.998000021537219 + 0.000000006148874 i) z	4
						+(-0.998000047158864 + 0.000000001537219 i) z	3
						+(-2.003000011557218 -0.000000006148874 i z	2
						+(3.002999035178864 -0.000000001537219 i) z
						+(-1.001000008472781 + 0.000000006148874 i)
			q	=	(1.000000000000000 + 0.000000008719619 i) z	8
						+(-0.000000008719619 -0.000000000000000 i) z	7
						+(0.001000030000000 -0.000000008719619 i) z	6
						+(3.000000008719619 + 0.000000000000000 i) z	5
						+(-1.999000229970001 + 0.000000008719619 i) z	4
						+(0.003000081280381 -0.000000000000000 i) z	3
						+(-0.001000030200006 -0.000000008719619 i) z	2
						+(-2.997000081190381 + 0.000000000000000 i) z
						+(0.999000229770005 + 0.000000008719619 i)
	f	=	(z	2 -1.001) (z	2 + 1.00000001) (z	3 + 2 z	2 -2.999999 z + 1)
		=	1.000000000000000 z	7 + 2.000000000000000 z	6 -3.000998990000000 z	5
			+0.998000020000000 z	4 -0.998000041009990 z	3 -2.003000010020000 z	2
			+3.002999029029990 z -1.001000010010000
	g	=	(z	2 -0.999)(z	2 + 1.00000003)(z	4 + 3z -1.0000002)
		=	1.000000000000000z	8 + 0.001000030000000 z	6 + 3.000000000000000 z	5
			-1.999000229970000 z	4 + 0.003000090000000 z	3 -0.001000030200006 z	2
			-2.997000089910000 z + 0.999000229770006


Nicolas Obin
email: nobin@ircam.fr

Volker Dellwo
email: v.dellwo@ucl.ac.uk

Anne Lacheret
email: anne.lacheret@u-paris10.fr

Xavier Rodet
email: rodet@ircam.fr

Expectations for Discourse Genre Identification: a Prosodic Study

Keywords: discourse genre, speaking style, prosody, perception, speech synthesis

Speech can be divided into discourse genres based on the contextual environment it occurs in (e.g. political speech, sport commentary speech, etc.). The present study investigated whether listeners can distinguish between speech from different discourse genres on the basis of acoustic prosodic cues only 1 . In a perception experiment with delexicalized speech 70 listeners with varying experience in French (native speakers, nonnative speakers, and non-speakers) were asked to identify four different types of discourse genres (church service, political, journal, and sport commentary). Results revealed a fair identification ability with a significant increase in performance with increasing experience in French. Identification confusion was used to cluster discourse genres according to their perceptual similarity.

Introduction

The concept of discourse genre (DG) has been studied widely in rhetoric and literature and more recently extended to the oral domain ( [START_REF] Halliday | Spoken and written Language[END_REF][START_REF] Biber | Variation Across Speech and Writing[END_REF]). It is generally hypothesized ( [START_REF] Benvéniste | Problème de linguistique générale[END_REF][START_REF] Bakhtine | Esthétique de la création verbale[END_REF])that each situation and each given social context correspond with a specific mode of production which is associated to specific formal markers that bear the traces at all levels (semantic, syntactic, phonological) of the DG. Following these studies, research in textual typology aims to : 1) describe the diversity of discourses (literary , legal, political, religious, etc.) ; 2) understand their classification into genres ( [START_REF] Rastier | Sens et textualité[END_REF]) ; 3) estimate their formal markers, in particular the co-occurrence of specific cues that can be considered as being typical of a genre. The challenge is to provide a robust and shared DG's typology. However, it remains very difficult to go further than conventional generic types (private, professionnal, public speech) and subdivisions (face to face conversation, phone conversation, public debates, radio and TV broadcasts, unprepared vs. planned speech, etc. [START_REF] Biber | Variation Across Speech and Writing[END_REF]). In the absence of a comprehensive representation of discourse genres and classes, each domain defines specific classification criteria which best account for its purpose (social, language activity, and formal spaces [START_REF] Branca-Rosoff | Types, modes et genres : entre langue et discours[END_REF][START_REF] Koch | Langage parlé et langage écrit[END_REF]). In strictly phonetic terms, studies focus on the description of phonostyles ( [START_REF] Obin | A method for automatic and dynamic estimation of discourse genre typology with prosodic features[END_REF][START_REF] Simon | Les phonostyles : une description prosodique des styles de parole en franc ¸ais[END_REF]). In particular, public discourse (such as po-litical, religious, journalistic and sport), considered as cultural stereotypes, are related to expressive strategies that act as markers of a phonostyle ( [START_REF] Lacheret-Dujour | La prosodie du franc ¸ais[END_REF]). Methods have been proposed recently to relate achievements in phonostylistic ( [START_REF] Fonagy | La vive voix[END_REF]) to text segmentation and classification in semantic in order to provide a unified interpretation framework : thematic content and informational sequence, semanticsyntactic structures and prosodic patterns (for recent work on French, see [START_REF] Lacheret | Schématisation discursive et schématisation intonative : question de genre ?[END_REF][START_REF] Degand | Mapping prosody and syntax as discourse strategies : how basic discourse units vary across genres[END_REF]).

This study aims at : 1) assessing whether the acoustic prosodic patterns of DGs are perceptually salient ; 2) estimating a classification of DG's on the basis of perceptual proximity ; 3) providing a reference for the evaluation of a speaking style speech synthesis system. Methodologically, a perceptual experiment was carried out in which listeners of different native language backgrounds had to identify different French DGs based on prosodic patterns. The paper is organized as follows : section 2 presents the design of speaking style corpus from which the stimuli for the present experiment were derived ; section 3 presents the experimental setup ; finally, results are presented and discussed in sections 4 and 5 respectively.

Speaking Style Corpus Design

Corpus Design

For the purpose of speaking style speech synthesis, a 4 hours French multi-media corpus was designed from which the stimuli for the present perceptual experiment were selected. The corpus consists of four different DG's : catholic mass ceremony, political, journalistic, and sport commentary. In order to reduce the DG intra-variability, the different DGs were restricted to specific discourse contexts (see list below) and to male speakers only. The following is a description of the four selected DG's :

• mass : christian church sermon (pilgrimage and sunday highmass sermons) ; single speaker monologue, no interaction. TAB. 1 -Description of the speaking style speech corpus.

The corpus was enriched in particular with the following linguisitic annotations : orthographical transcription ; phonemic alignement and breath detection using ircamAlign [START_REF] Lanchantin | Automatic phoneme segmentation with relaxed textual constraints[END_REF], then manually corrected ; syllabification on inter-pausal groups.

Experimental Design

The experiment consists in a multiple choice DG's identification task based on prosodic cues. It was not possible to control the linguistic content of the speech utterances which is an evident cue for DG's identification (a single keyword may sometimes be sufficient to identify a DG). For this reason the stimuli were delexicalized using low pass-filtering. Stimuli were utterances that were extracted from the speaking style speech corpus (see section 2), filtered to remove semantic access, and then presented as a multiple choice identification experiment to multi-language speakers in a crowd-sourcing framework.

Subjects

70 subjects participated in this experiment : 37 native French speakers, 20 French speaker, 15 non-French speakers ; 46 expert subjects, 24 naïve subjects. Expert subjects had a variety of backgrounds : speech and audio technologies, linguistic, musicians. 7 subjects were removed because they produced mistakes in using the experimental interface. In the case of multiple participation of a subject, his first participation was used for analysis only. Subjects were aged from 20 to 65 years, with a strong proportion (65%) within the 20-35 year range.

Stimuli

40 speech utterances (10 per DG) were selected from the speaking style corpus. Such selection was conducted to provide various and representative prosodic patterns for each DG. Firstly, segmentation into speech utterances was accomplished according to the prosodic period. Prosodic period is defined as being a sequence of inter-pausal groups which ends with a conclusive frontier (combination of three acoustic criteria : low pitch accent, long syllable duration, followed by a long pause). It is in particular supposed to be related to discursive speech units (a prosodic object as a marker of a communicative and discursive object). However analysis of speech samples reveals that the definition of the prosodic period fails to account for some specific speech sequences or more generally to some DG. In particular, it was regularly observed that journalistic utterances do not end with a major frontier (low pitch accent which is not followed by a pause, or even high-pitch accent which is followed or not by a pause, or inconsiderable long prosodic periods). This was especially observed for spontaneous speech such as the sport commentary where prosodic periods could not be directly associated to discursive units. Consequently, segmentation into speech utterances was finally decided by an expert linguist according to discursive-prosodic cues in such conflictual cases. In particular, speech segmentation was chosen as being prosodic objects which were not necessarily formally defined by prosodic constraints only but combining prosodic cues to weak discursive, syntactic and semantic dependencies to the discursive context. Secondly, the selection criteria were derived from a classification of speech utterances into discursive sequences as well as prosodic structure and complexity. In particular, archetypal speech utterances were selected depending on the DG (mass : "au nom du père et du fils, et du Saint-Esprit, ainsi soit-il", in the name of the Father, the Son, and the Holy Spirit, Amen. ; political : "mes chers compatriotes, vive la République et vive la France", my fellow countrymen, long live the republic ! Long live France ! ; sport commentary : "oh le but de Babel ! le but de Babel ! le but de Babel !", What a goal by Babel ! Goal by Babel ! Goal by Babel ! ; journal : no specific speech utterance was observed). Then, speech utterances were classified into discursive sequences depending on the DG. For instance, journalistic chronicles can be formally described as a sequence of topic sequences with punctual interaction with a lead speaker during topic changes. Speech utterances were thus classified into global introduction from a lead speaker ("l'Eco du jour" : c'est l'actualité économique de ce lundi 26 octobre 2009 vue par Philippe Lefébure, "Eco du Jour" is the economic news for today, Monday, October 26, 2009 hosted by Philippe Lefébure), and initial, medium, terminal and transitional sequences for each topic (initial : "MacDo va quitter l'Islande : conséquence directe de la crise, raconté dans L 'Humanité" , MacDonald's is quiting Iceland : a direct result of the financial crisis, story in "l'Humanité". ; conclusive : "les franc ¸ais disent qu'on va dans la mauvaise direction, une seule réponse du côté de l'UMP : il faut y aller plus vite", The French say we have taken the wrong path, the only reply from the UMP is : "We must go faster". ). Sport commentary sequences were classified according to actional, situational (current action, past action, off-line comments), and emotional (more or less intense) criteria. Other DG's speech utterances were classified in the same manner. Such classification interestingly relates to specific prosodic patterns. Then various prosodic sequences were chosen for each DG (in particular : low, medium, and high terminal pitch accent, intermediate lexical pitch accents as well as hesitations). Thirdly, as it was observed that speech utterance's duration strongly depend on the DG, speech utterances were classified into short (4 ± 0.5s.) and long (10 ± 1s.) utterances that were homogenously distributed for each DG. Finally, 2 speech utterances were selected for each speaker to remove any identification based on speaker recognition.

Then, speech samples were processed as follows : a) background noise and reverberation were minimized using a noise cancelation algorithm ( [START_REF] Bogaards | An interface for analysis-driven sound processing[END_REF]) ; b) semantic access was removed using a band-pass filter. A pass-band was chosen that insured that the lowest frequency of the fundamental frequency and the highest frequency of its first harmonic was included ( [START_REF] Bogaards | An interface for analysis-driven sound processing[END_REF]). This was done to extract speech prosodic characteristics only c) active speech mean level normalization at -20dBov [START_REF]Measuring speech activity[END_REF] ; d) speech samples were compressed in mp3 format at 192Kb/s.

Procedure

The experiment consists of a multiple choice identification task from speech prosody perception. It was conducted according to source-crowding technique using web social networks 3 . Subjects were given a brief description of the different speaking styles without being exposed to actual speech examples. This approach was adopted in order to focus subjects on their own mental representation of the different speaking styles and their expected prosodic cues. They were asked to associate a speaking style to each of the speech samples 4 . For this purpose, subjects were given three options :

• total confidence : select only one speaking style when certain of the choice ; • confusion : select two different speaking styles when two speaking styles are possible ; • total indecision : select "indecision" when completely unsure. Subjects were asked to use this possibility only as a very last resort. Additional informations was gleaned from the participants : speech expertise (expert, naïve), language (native French speaking, French speaking, non-French speaking), age, and listening condition (headphones or not). Subjects were encouraged to use headphones.

Results

Identification performance was estimated using a newly developed measure based on Cohen's Kappa statistic. Cohen's Kappa provides agreement between two raters in the case of caterogical rating [START_REF] Cohen | A coefficient of agreement for nominal scales[END_REF]. Our measure monitors the agreement between subjects' ratings and the correct answer. The resulting Kappa values are considered as a measure of identification performance. The measure varies from -1 to 1 : -1 is perfect disagreement ; 0 is chance ; 1 is perfect agreement.

Overall score reveals fair identification performance (K=0.45). Figure 1 shows that there are differences according to native language background of the listeners. French natives perform better than non-native French-speakers ; non-French speakers perform only slightly. ANOVA analysis (one-way analysis of variance [START_REF] Hair | Multivariate data analysis[END_REF]) was conducted to assess whether identification performance depends on the language of the subjects. Analysis reveals a significant effect of the language (F (2, 59) = 15; p < 0.001). Post-hoc analysis reveals significant difference between native French speakers and the others (F (1, 52) = 13; p < 0.001 , F (1, 43) = 24; p < 0.001) but no effect between non-native French speakers and non-French speakers (F (1, 23) = 3; p = 0.07) (fig. 1). Investigating the results in finer detail reveals that identification performance significantly depends on the DG. Substantial identification performance was observed for sport commentary (K=0.7) ; fair identification for the journalistic discourse (K=0.54) ; and only slight identification for mass discourse and political discourse (K = 0.38 and 0.34 respectively). 3 Interestingly, identification performance reveals different configurations of the language effect depending on the DG (fig. 2). For the mass, a significant difference between native French speakers and the others was observed (F (1, 52) = 6.9; p = 0.01, F (1, 43) = 5.8; p = 0.02), but none among these (F (1, 23) = 0.19; p = 0.7). For the political discourse, a significant difference between all language pairs was observed : native French vs. French speaking (F (1, 52) = 11; p = 0.001), native French vs. non-French speaking (F (1, 43) = 30.5; p < 0.001 ) as well as French speaking vs. non-French speaking (F (1, 23) = 7.8; p = 0.01). For the journalistic discourse, no significant difference between native French and French speaking (F (1, 52) = 1; p = 0.3) and a significant between these and non-French speaking (F (1, 43) = 13.5; p < 0.001,F (1, 23) = 5.1; p = 0.03) were observed. For the sport commentary, there was a significant difference between native French speakers and the other (F (1, 52) = 14; p < 0.001,F (1, 43) = 17; p < 0.001), but none between these (F (1, 23) = 0.8; p = 0.37). Furthermore, it was noticed that non-French speaking participants provides slight and random identification performance in the case of journalistic and political discourse respectively.

Ircam Analysis and Synthesis Perceptual

Finally, multi-dimensional scaling (MDS, [START_REF] Borg | Modern Multidimensional Scaling : theory and applications[END_REF]) and hierarchical clustering ( [START_REF] Hastie | The Elements of Statistical Learning[END_REF]) methods were used to represent and estimate DG's similarity according to the observed perceptual confusion. For each speech sample, the observed confusion matrix was used to define speech utterance coordinates. A similarity distance between speech utterances was then estimated according to the city-block metric. The set of pairwise speech utterances similarity was then used to represent speech utterances into a 2-dimensional space according to multi-dimensional scaling. Finally, DGs were clustered using the complete linkage method. In parallel, the 4 DG's have been discussed and rated by two expert linguists after the Koch's conceptual scale ( [START_REF] Koch | Langage parlé et langage écrit[END_REF]) with 3 degrees. Then DGs were clustered according to this conceptual description in the same manner. Figure 3 represents speech utterances into the resulting 2-dimensional space and the comparison of the resulting DG's typology.

Discussion

The experiment reveals fair identification performance from prosodic perception. This confirms evidence for the hypothesis that DGs relate to prosodic patterns that do not depend on a specific speaker and that are shared by a group of listeners. Not suprisingly, these prosodic objects dramatically depend on the language backgrounds. This shows that DG's abstract prosodic representation depends at least on the language. It could be hypothesized that such representation could more generally depend on language classes or even on culture 5 . More precisely, the experiment reveals that the language factor significantly depends on the DG. This suggests that DG relate to prosodic objects that could be more or less shared regardless of the language. In particular, prosodic cues related to sport commentary appeared clearly almost common across all languages while prosodic cues related to political discourse and mass discourse were language specific. Interesting intermediate identification performance of nonnative French speakers was observed : if identification performance is systematically situated between native French speakers and non-French speakers, different grouping were observed, depending on the DG (grouping with non-French speakers for mass discourse and sport commentary, grouping with native French speakers for the journalistic discourse, and no grouping for the political discourse). This suggests a native language effect as well as a cultural background dependency. A comparison of DG clusters as estimated from the conceptual classification and from prosodic perception revealed a similar cluster structure (fig. 3). This confirms that discourse context (situational, spatio-temporal, ... context) consistently relate to prosodic strategies. Moreover, prosodic clusters precise the perceptual distance that in particular clearly distinguishes journalistic discourse from political discourse and mass discourse on the prosodic dimension. This result supports the hypothesis that prosodic strategies act as markers of a specific speech act ( [START_REF] Searle | Speech Acts[END_REF]) (for instance : neutrally describing an event with distanciation 5 such hypothesis is supported by non-French speaking participants which comment that they could not represent themselves "how sounds" a Christian sermon (religious dependency) nor political new year's speech (cultural dependency) for the journalistic discourse vs. arguing and persuading for the political discourse and mass discourse). Sport commentary stands significantly apart from the other DGs. This confirms previous studies on the very specific nature of the sport commentary ( [START_REF] Deulofeu | Les commentaires sportifs constituent-t-ils un genre ? au sens linguistique du terme ?[END_REF]), in particular in its iconical dimension : sportscasters do not only describe but vocally mimic the action being observed. This is even more true in the case of radio sport commentary, where sportcaster must supply the absence of the image media.

Conclusion

A perceptual experiment on the identification of discourse genres on a speech prosodic perception basis was proposed. Identification performance confirms evidence for the hypothesis that DGs relate to prosodic patterns that do not depend on a specific speaker and that are shared among listeners. When overall factorial analysis reveals a significant language effect, it is clearly dependent on the DG. A comparison between DG clusters obtained from a conceptual description and perceptual confusion indicates that discourse context consistently relates to specific prosodic strategies. DG perceptual clusters even precise and suggest other discursive effects to explain observed differences of prosodic configuration. In a further study, these results will be used as a reference identification performance to evaluate a speaking style speech synthesis system.

FIG. 2 -

 2 FIG.2-Identification performance according to the subject's language background for each discourse genre. From top to down and left to right : mass, political, journalistic, sport.

FIG. 3 -

 3 FIG. 3 -(a) Representation of the speech utterances according to their similarity according to Multi Dimensional Scaling. (b) DG's clusters from : conceptual scale (top) and perceptual confusion (down)

This study was supported by ANR Rhapsodie 07 Corp-030-01 ; reference prosody corpus of spoken French ; French National Agency of research

; 2008-2012.

This was reached with the exception of the sport commentary which has half duration than the other DG's


Nicolas Obin
email: nobin@ircam.fr

Pierre Lanchantin
email: lanchantin@ircam.fr

Mathieu Avanzi
email: mathieu.avanzi@unine.ch

Anne Lacheret-Dujour

Xavier Rodet
email: rodet@ircam.fr

TOWARD IMPROVED HMM-BASED SPEECH SYNTHESIS USING HIGH-LEVEL SYNTACTICAL FEATURES

Keywords: HMM-based speech synthesis, Prosody, High-Level Syntactical Analysis

A major drawback of current Hidden Markov Model (HMM)-based speech synthesis is the monotony of the generated speech which is closely related to the monotony of the generated prosody. Complementary to model-oriented approaches that aim to increase the prosodic variability by reducing the "over-smoothing" effect, this paper presents a linguistic-oriented approach in which high-level linguistic features are extracted from text in order to improve prosody modeling. A linguistic processing chain based on linguistic preprocessing, morpho-syntactical labeling, and syntactical parsing is used to extract high-level syntactical features from an input text. Such linguistic features are then introduced into a HMM-based speech synthesis system to model prosodic variations (f0, duration, and spectral variations). Subjective evaluation reveals that the proposed approach significantly improve speech synthesis compared to a baseline model, event if such improvement depends on the observed linguistic phenomenon.

INTRODUCTION

Research on speech synthesis has brought significant improvements over the past decade that makes possible to generate natural speech from text [START_REF] Yoshimura | Simultaneous modeling of spectrum, pitch and duration in hmm-based speech synthesis[END_REF][START_REF] Black | Statistical parametric speech synthesis[END_REF]. However, if the synthesized speech sounds acoustically natural, it is often considered poor according to the way of speaking (prosodic artefacts and monotony). Now, modeling the variability in the way of speaking (variations of prosodic parameters) is required to provide natural expressive speech in many applications of high-quality speech synthesis such as multi-media (avatar, video game, story telling) and artistic (cinema, theater, music) applications. Despite growing attention to prosody modeling over the past few years, one of the major drawbacks of actual prosody models remains the monotony of the generated prosodic parameters. The prosody monotony is both related to poor dynamics ("averaging problem") and poor variability of the generated prosodic parameters. Poor variability basically results in the generation of stereotypical prosody and is mainly due to the lack of linguistic knowledge extracted from the text . In order to model the variability of prosodic parameters, it is necessary to determine an appropriate representation of prosodic parameters and then to extract and estimate the effects of high-level linguistic features (syntactic, semantic, discursive) on the observed prosodic parameters. At the signal level many approaches have been proposed to represent the variations of the prosodic parameters ( [START_REF] Yoshimura | Simultaneous modeling of spectrum, pitch and duration in hmm-based speech synthesis[END_REF][START_REF] Latorre | Multilevel parametric-base f0 model for speech synthesis[END_REF] for fundamental frequency variations and [START_REF] Zen | Hidden semi-markov model based speech synthesis[END_REF][START_REF] Gao | Duration refinement by jointly optimizing state and longer unit likelihood[END_REF][START_REF] Obin | A multi-level contextdependent prosodic model applied to durational modeling[END_REF] for durations) based on statistical parametric modeling. At the symbolic level, prosodic variations are affected in many ways by linguistic constraints (phonologic, syntactic, semantic, discursive,...) occurring on a set of more or less linguistically well-defined units. Strictly speaking of syntactic-prosodic relationships, syntactic structure does not affect prosodic structure [START_REF] Black | Assigning intonation elements and prosodic phrasing for english speech synthesis from high level linguistic input[END_REF][START_REF] Sridhar | Exploiting acoustic and syntactic features for automatic prosody labeling in a maximum entropy framework[END_REF] only, but acoustic variations as well, as shown in recent studies (for instance: oral parenthesis [START_REF] Dehé | Clausal parentheticals, intonational phrasing, and prosodic theory[END_REF][START_REF] Gachet | Les parenthèses en franc ¸ais : Etude prosodique[END_REF]). This paper presents an exploratory study that aims to incorporate high-level syntactical features extracted from text within a HMM-based speech synthesis framework in order to improve prosodic model variability. The paper is organized as follows: section 2 presents the linguistic processing chain and the syntactical features extracted from text; section 3 introduces the HMM-based synthesis framework; thus perceptual experiment and results are presented in sections 4 and 5.

HIGH-LEVEL SYNTACTICAL ANALYSIS

Linguistic Processing Chain

At the symbolic level, the input text (sentence, set of sentences or raw text) is processed by an automatic linguistic parser (Alpage Linguistic Processing Chain 1 ) in order to extract high-level linguistic features (morpho-syntactical and syntactical) within the sentence level. The Alpage Linguistic Processing Chain is a full linguistic processing chain for French based on a sequence of processing modules: a lexer module, a parse module, and a post-processing module. The input text is first preprocessed by the lexer module that output Direct Acyclic Graphs (DAGs) combined with lexical information retrieved from a morphological and syntactic lexicon. Then deep parsing is performed by the parser (FRMG), a symbolic parser based on a compact Tree Adjoining Grammar (TAG) for French that is automatically generated from a meta-grammar. The parsing result is then enriched by a series of post-processing modules whose role is to organize all the informations retrieved along those steps. The output of the parser is a shared derivation forest that represents all derivation structures that the grammar can build for the input sentence, and indicates which TAG operation (substitution, adjunction, anchoring) took place on a given node of a given tree for a given chunk. This forest is then transformed into a shared dependency forest: anchors of trees related to a given node label are put into a dependency relationship with this label. A dependency forest is represented into a DEP XML format that incorporates the following items:

• clusters that are associated to the forms of the sentence;

• nodes that point on a given cluster and associated to a lemma, a syntactical category and a set of derivations;

• edges that connect a source node with a target node with a given label. More precisely, a given edge is associated with a set of derivations related to this edge and the related source and target chunk operations.

At last the forest is disambiguated by a heuristic-based module that outputs a unique dependency tree. An example of a disambiguated dependency graph as provided by the parser is shown in figure 1.

Syntactical Features Extraction

From the output of the linguistic process described in 2.1, a set of more or less high-level linguistic features were extracted.

morpho-syntactical form features are extracted from the surface processing.

• form segment;

• form lexical category and class (function vs. content form);

form dependencies are extracted from the deep parse processing. This set basically encodes the relationship between forms.

• {governor, current, governee} form lexical category and class;

• edge type and label between current form and {governor, gov-ernee} form;

• signed dependency distance between current form and {governor, governee} forms (in forms and in chunks);

recursive chunk are retrieved in a top-down process according to the operations and associated derivations. As a chunk of a given level can have several derivations (i.e governee chunks), we chose to stack governee chunk of a given level from left to right in order to provide a binary tree chunk representation. For our example sentence (cf. fig. 1), the transformed recursive chunks are:

(S (AdvP Longtemps ) ( (VP je me suis couché ) (NP de bonne heure ) ) )

Recursive chunks were finally transformed into non-recursive chunks by extracting only leaves of the transformed chunk tree (and encoding only partially recursivity through chunk depth estimation).

• {governor, current, governee} chunk category;

• edge type and label between current chunk and {governor, gov-ernee} chunks;

• signed dependency distance between current chunk and {governor, governee} chunks (in forms and in chunks);

• chunk depth;

adjunction represents a specific type of syntactical phenomena. In particular, adjunctions could concern different text spans (from a single form to a full sentence). Adjunction covers a large amount of syntactical phenomena such as incises, parentheses, subordinate and coordinate clauses, enumerations, ...

In the FRMG parser formalism, adjunctions can be easily extracted according to some specific pattern matching (including modified, introducer, and modifier nodes). Full adjunction is thus extracted by retrieving the full dependency descendance from the introducer. Modifier, introducer and modified category as well as edges type and label are then used to identify adjunction's type (such as relative clauses, adjective clauses, adverb clauses, ...).

• {modified, introducer, modifier} form lexical category and class;

• adjunction type retrieved from an adjunction dictionary;

• edge type and label between modified and introducer nodes and between introducer and modifier nodes; • signed dependency distance between the adjunction's introducer and the modified node (in forms and in chunks);

As adjunctions have recursive property (a given adjunction can be embedded within another adjunction, or within an adjunction of the same type), features were extracted separately for each known adjunction type (with related pattern). Then, in case of recursivity of a given adjunction's type, adjunction with the larger span was extracted only.

INTEGRATION IN THE HMM-BASED SPEECH SYNTHESIS

Over the last decade, HMM-based speech synthesis system has grown in popularity [START_REF] Yoshimura | Simultaneous modeling of spectrum, pitch and duration in hmm-based speech synthesis[END_REF][START_REF] Black | Statistical parametric speech synthesis[END_REF][START_REF] Lanchantin | A HMM-based speech synthesis system using a new glottal source and vocal-tract separation method[END_REF]. This system models spectrum, excitation, and durations in a unified HMM framework. Compared to unit-selection systems it offers the ability to model different voices, speaking styles or emotional speech without requiring the recording of large databases. The implementation of our HMM-based speech synthesis system is based on the HTS Toolkit.

During the training, both spectrum and excitation parameters vectors (including f0) and their dynamic vectors are extracted from the speech corpus and used to train the context-dependent HMMs according to the alignment. In this way, a mapping is performed between speech acoustics and linguistic features. Due to the large amount of context-dependent models, the models are hierarchically clustered into acoustically similar models using decision trees estimated according to Maximum-Likelihood Minimum Description Length criterion (ML-MDL). Multi-space probability distributions (MSD) are used to model variable dimensional parameter sequence such as log f0 with unvoiced regions properly. Each contextdependent HMM has state duration probability density functions (PDFs) to capture the temporal structure of speech [START_REF] Zen | Hidden semi-markov model based speech synthesis[END_REF][START_REF] Gao | Duration refinement by jointly optimizing state and longer unit likelihood[END_REF].

During the synthesis, the text to be synthesized is first converted to a context-dependent label sequence. An utterance HMM is then built by concatenating the most appropriate context-dependent HMMs according to the label sequence and the trained decision trees. State durations of the utterance HMM are then determined based on the state duration PDFs. Finally, the speech parameter are generated using a maximum likelihood algorithm including dynamic features [START_REF] Tokuda | An algorithm for speech parameter generation from continuous mixture HMMs with dynamic features[END_REF] from which a speech waveform is synthesized using a speech synthesis filter.

It is important to note that by adopting this HMM framework, all prosodic dimensions could be estimated according to the linguistic contextual features (i.e. not only f0 variations and duration but voice quality as well). The following linguistic feature sets were introduced as contextual features:

• baseline linguistic units: phonem, syllable and prosodic group, with phonem and syllable phonological features (phonem phonological features, syllabic structure, and prosodic structure: prosodic frontiers and syllable prominence) [START_REF] Obin | A multi-level contextdependent prosodic model applied to durational modeling[END_REF]; Fig. 1. Disambiguated dependency graph for the sentence: "Longtemps, je me suis couché de bonne heure." ("For a long time I used to go to bed early"). Squares represent clusters (with associated form), circles represent nodes (with associated lemma and lexical category), arrows represents edges (with associated dependency label) going from source node (governor) to target node (governee).

• morpho-syntactical linguistic unit: form;

• dependency linguistic unit: form;

• chunk linguistic unit: chunk;

• adjunction linguistic unit: adjunction;

For each feature set, low-level linguistic features were extracted on each linguistic unit with left-to-right contexts (second order for the phoneme features, and first order for the other linguistic units with the exception of adjunction that does not necessarily have direct adjacent context): locational (position of a given unit within higher level units) and weight (number of observations of a given linguistic unit within higher level units) features.

EXPERIMENT

Experimental setup

For the purpose of that experiment, two models were compared: a baseline model including baseline and morpho-syntactical linguistic features; and a rich linguistic model including all linguistic features. Such linguistic sets were used to train all speech synthesis parameters (spectral and f0 variations and durations). Models were trained on a 1 hour (956 sentences) laboratory corpus spoken by a French non professional male speaker and recorded at 16kHz in an anechoic room.

The following processing chain was applied to the training corpus: phonemic segmentation using ircamAlign [START_REF] Lanchantin | Automatic phoneme segmentation with relaxed textual constraints[END_REF]; syllabation on interpausal groups; automatic syllable prominence detection using ircam-Prom [START_REF] Obin | A syllable-based prominence model based on discriminant analysis and context-dependency[END_REF]. Linguistic feature were extracted from text with the Alpage Linguistic Processing Chain. All analysis were conducted within the ircamCorpusTools [START_REF] Veaux | Ircamcorpustools: an extensible plateform for speech corpora exploitation[END_REF] framework. Due to the high complexity of the linguistic structure as well as their dependencies with speech prosody, estimating precisely the influence of each of the extracted linguistic features on the synthesized speech and prosody is unreachable. It is more likely to reduce the evaluation to a subset of well-defined specific syntactical phenomena. For that reason the evaluation corpus was designed according to syntactical adjunctions only. This choice is motivated by several reasons: adjunctions concerns long-term speech scope (oral parenthesis) which makes it easier to evaluate variations across different models; adjunctions could be described with limited vocabulary which makes easier to generate a limited sentence corpus which is needed in a subjective evaluation experiment; adjunctions has been shown in a previous study to be the most relevant syntactical feature for prosody modeling [START_REF] Obin | Using high-level syntactical features to improve prosodic modeling[END_REF].

The text corpus used for the evaluation has been designed in the following maneer: 10 baseline sentences were chosen without any specific syntactic property (ex: Le chat a mangé la souris, The cat ate the mouse). These sentences were then enriched according to a set of adjunction type (subordinate participial and relative clauses, coordinate clauses, incises and enumerations). For each adjunction type, original sentences were enriched according to two control parameters: position (initial, medial, final) and complexity (presence or not of adjunctions within the current adjunction) of the introduced adjunction. This finally lead to a 54 sentences evaluation corpus.

For the feasibility of the study, only 20 were selected for the subjective evaluation. Finally speech parameters were synthesized for each sentence according to the considered model and contextual features sequence extracted from text (prosodic structure being inferred from a common model then equal for both compared models).

Subjective evaluation

The evaluation consists of a subjective comparison between the 2 models. A comparison category rating (CCR) test was used to compare the quality of the synthetic speech generated by baseline and rich linguistic models. The evaluation were conducted according to source-crowding using social networks 2 . 50 French native speakers (including 17 experts and 33 naive listeners) compared a total of 20 speech sample pairs. They were asked to attribute a preference score according to the quality of each of the sample pairs on the comparison mean opinion score (CMOS) scale. Listening test was performed with headphones.

RESULTS & DISCUSSION

Analysis of variance (ANOVA) has been estimated on the resulting preference distributions in order to asses perceptive differences be- ). A closer look to the results (fig. 2) reveals noticeable differences among the different sentence sets. There is a slight but significant difference (F (1, 499) = 8.4, pvalue = 4.10 -3 ) for the baseline sentences. This is actually in agreement to what could be expected since such sentences were defined without any particular syntactical property. Scores obtained for the enriched sentences reveals the strongest preference to the rich linguistic model, this result is not observed for all types of adjunction. If there is a clear preference to the rich linguistic model for the last three adjunction types (coordinate: F (1, 299) = 50, pvalue < 10 -5 ; incise: F (1, 399) = 38, pvalue < 10 -5 ; enumeration: F (1, 99) = 78.4, pvalue < 10 -5 ), there is however a preference to the baseline model in case of subordinate sentences (clear preference in case of participial subordinate: F (1, 399) = 17.58, pvalue < 10 -5 and no preference in the case of relative subordinate (F (1, 299) = 1.5 with pvalue = 0.2). The prosodic improvements of the synthesized speech can be listed as follows:

• accentuation: more prosodic variations observed on the accented syllables; accent form modification (conclusive, continuative accent); • prosodic phrasing: local improvement of the prosodic phrasing (pitch variations and local speech rate).

These are encouraging preliminary results even if such improvements are not systematic: this is probably due to sparse observations (number of occurrences for each extracted syntactical phenomenon) for which parameters of the model could not be robustly estimated. Further experiments with more training observations will be carry out.

CONCLUSION & FURTHER WORK

We have presented a speech synthesis system using high-level linguistic features in order to improve prosody modeling. Perceptual evaluation shows that the proposed approach improve speech synthesis compared to the state-of-the-art approach. However this improvement still depends on the observed linguistic phenomenon which is probably due to sparse observations. In further work such approach will be evaluated with more linguistic observations and according to the linguistic complexity of an observed corpus. A joint model approach will be used in order to improve the estimation of the dependencies between high-level linguistic features and prosodic variations by jointly modeling prosodic variations over different time scales. Furthermore, one needs to define evaluation measurement that could more precisely account for prosodic variability.

Fig. 2 .

 2 Fig. 2. CMOS preference score for the different adjunction types

  2 this perceptual evaluation group is available on http://www. facebook.com/group.php?gid=150354679034&ref=ts and the test on http://recherche.ircam.fr/equipes/ analyse-synthese/lanchant/index.php/Main/TestSP.

http://alpage.inria.fr/alpc.en.html

ACKNOWLEDGMENTS

This study was partially supported by: • ANR Rhapsodie 07 Corp-030-01; reference prosody corpus of spoken French; French National Agency of research (ANR); 2008-2012.; • Programmes Exploratoires Pluridisciplinaires (PEPS), CNRS/ST2I, 2008-2010; • Swiss National Science Foundation, grants n • PBNEP1-127788 and n • 100012-113726/1, Neuchâtel University.