A Ammar

F Chinesta

A Falcó

On the Convergence of a Greedy Rank-One Update Algorithm for a Class of Linear Systems

In this paper we study the convergence of the well-known Greedy Rank-One Update Algorithm. It is used to construct the rank-one series solution for full-rank linear systems. The existence of the rank one approximations is also not new, but surprisingly the focus there has been more on the applications side more that in the convergence analysis. Our main contribution is to prove the convergence of the algorithm and also we study the required rank one approximation in each step. We also give some numerical examples and describe its relationship with the Finite Element Method for High-Dimensional Partial Differential Equations based on the tensorial product of one-dimensional bases. We illustrate this situation taking as a model problem the multidimensional Poisson equation with homogeneous Dirichlet boundary condition.

Introduction

In [START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in Kinetic Theory modelling Complex Fluids[END_REF][START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in Kinetic Theory modelling Complex Fluids. Part II: Transient simulations using space-time separated representations[END_REF], some of the authors of the present paper propose the use of a separated representation, which allows to define a tensor product approximation basis as well as to decouple the numerical integration of a high dimensional model in each dimension. The milestone of this methodology is the use of shape functions given by a tensorial based construction. This fact has advantages as the manipulation of only one dimensional polynomials and its derivatives, that provides a better computational performance and simplified implementation and use one-dimensional integration rules. Moreover, it makes possible the solution of models defined in spaces of more than hundred dimensions in some specific applications. This problem is closely related with the decomposition of a tensor as a sum of rank-one tensors, that it can be considered as a higher order extension of the matrix Singular Value Decomposition.

The purpose of this work is to formalize and analyze the above strategy in the framework of methods for solving linear systems by means tensor decompositions. As we will show the approximation given in [START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in Kinetic Theory modelling Complex Fluids[END_REF][START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in Kinetic Theory modelling Complex Fluids. Part II: Transient simulations using space-time separated representations[END_REF], is closely related with the best low-rank approximation problem for high order tensors (see [START_REF] De Silva | Tensor Rank and Ill-posedness of the best low-rank approximation problem[END_REF]). Unfortunately, in [START_REF] De Silva | Tensor Rank and Ill-posedness of the best low-rank approximation problem[END_REF] it has been proved that tensors of order 3 or higher can fail to have best rank-r approximation for r ≥ 2. Our strategy, with the perspective of [START_REF] De Silva | Tensor Rank and Ill-posedness of the best low-rank approximation problem[END_REF] in mind, is to use the fact that tensors of order 3 or higher have best rank-1 approximation.

In this context, we propose the use of a Greedy Rank-One Update Algorithm to construct, for a full rank linear system, a rank-r approximate solution. This approach is based in the so-called by the signal processing community as the Matching Pursuit Algorithm of Mallat and Zhang [START_REF] Mallat | Matching pursuit with time-frequency dictionaries[END_REF], also known as Projection Pursuit by the statistics community (see Friedman and Stuezle [START_REF] Friedman | Projection pursuit regression[END_REF] and Huber [START_REF] Huber | Projection pursuit[END_REF]) or as a Pure Greedy Algorithm (see the recent survey of Temlyakov [START_REF] Temlyakov | Greedy approximation[END_REF]) in the approximation theory community. Our main contribution, stated in Theorem 1, is to prove the convergence of this Greedy rank-one update algorithm. and characterize the speed of convergence in terms of a sequence of angles. This strategy depends strongly on the computation of the best rank-1 approximation of the residual obtained at each step of the proposed algorithm. To solve this we will propose the use of a Block Coordinate Descend Method, because it has global convergence. In particular, we will show that for the class of invertible matrices, this problem collapses, for each selected direction, to an ordinary least-squares problem. We remark that this strategy appears to be identical to so-called Alternating Least Squares (ALS) method proposed in [START_REF] Beylkin | Algorithms for numerical analysis in high dimensions[END_REF] for the class of separable matrices (see statement of Corollary 1 below).

In [START_REF] Kolda | Orthogonal tensor decompositions[END_REF] an orthogonal Greedy tensor decomposition has been used in order to compute a rank-r approximation. However, as the author points out, the computational difficulty of this approach arises in enforcing the constrains of the rank-one approximation needed at each step of the proposed algorithm. Zhang and Golub [START_REF] Zhang | Rank-one approximation to high order tensors[END_REF], also explores various computational techniques when the tensor has a completely orthogonal decomposition, in which case the problem is more simpler. All these methods use an ALS Approach for computing the rank-one approach. On the other hand, the approach following in this paper differs from the [START_REF] De Lathauwer | On the best rank-1 and rank-(R 1 , . . . , R N ) approximations of high order tensors[END_REF] based in the associated Lagrange Equations to the minimum least-squares cost function.

This paper is organized as follows. In the next section we introduce the notation used in this paper and give our main result, the convergence of the Greedy Rank-One Update Algorithm for solving full rank linear systems. In Sect. 3 we study the rank-one approach, in particular we prove that a Block Cyclic Coordinate Descend strategy implies an ALS Algorithm. Section 4.1 is dedicated to give some numerical examples of the above algorithms and describe its relationship with the Finite Element Method for High-Dimensional Partial Differential Equations based on the tensorial product of one-dimensional bases. We illustrate this situation taking as a model problem the multidimensional Poisson equation with homogeneous Dirichlet boundary condition. We conclude with some comments and remarks.

Definitions and Statement of Main Result

First at all we introduce some notation. We denote by R N ×M , the set of N × M-matrices and by A T the transpose of a given matrix A. As usual we use x, y = x T y = y T x to denote the Euclidean inner product in R N , and its corresponding 2-norm, by x 2 = x, x 1/2 . Let I N be the N × N -identity matrix and when the dimension is clear from the context, we simply denote it by I. Given a sequence {u j } ∞ j =0 ⊂ R N , we say that a vector u ∈ R N can be written as

u = ∞ j =0 u j if and only if lim n→∞ n j =0 u j = u
in the • 2 -topology. Now, we recall the definition and some properties of the Kronecker product. The Kronecker product of A ∈ R N 1 ×N 1 and B ∈ R N 2 ×N 2 , written A ⊗ B, is the tensor algebraic operation defined as

A ⊗ B = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ A 1,1 B A 1,2 B • • • A 1,N 1 B A 2,1 B A 2,2 B • • • A 2,N 1 B . . . . . . . . . . . . A N 1 ,1 B A N 1 ,2 B • • • A N 1 ,N 1 B ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ∈ R N 1 N 2 ×N 1 N 2 .
Also, it can be defined by

(A ⊗ B) (j 1 -1)N 2 +j 2 ;(i 1 -1)N 2 +i 2 = A j 1 ;i 1 B j 2 ;i 2 . ( 1
)
Assume that

A i ∈ R N i ×N i for 1 ≤ i ≤ d.
Proceeding inductively we show that

(A 1 ⊗ • • • ⊗ A d ) s,t = (A 1 ) j 1 ,i 1 • • • (A d ) j d ,i d
if and only if s and t satisfy

s = j d + d-1 l=1 (j l -1) d l =l+1
N l and

t = i d + d-1 l=1 (i l -1) d l =l+1
N l .

(

Finally, we list some of the well-know properties of the Kronecker product (see for example [START_REF] Graham | Kronecker products and matrix calculus with applications[END_REF] or [START_REF] Van Loan | The ubiquitous Kronecker product[END_REF]). The concept of separated representation was introduced by Beylkin and Mohlenkamp in [START_REF] Beylkin | Algorithms for numerical analysis in high dimensions[END_REF] and it is related with the problem of constructing the approximate solutions of some classes of problems in high-dimensional spaces by means a separable function. In particular, for a given map

1. A ⊗ (B ⊗ C) = (A ⊗ B) ⊗ C. 2. (A + B) ⊗ (C + D) = (A ⊗ C) + (B ⊗ C) + (A ⊗ D) + (B ⊗ D). 3. AB ⊗ CD = (A ⊗ C)(B ⊗ D). 4. (A ⊗ B) -1 = A -1 ⊗ B -1 . 5. (A ⊗ B) T = A T ⊗ B T . 6. If
u : [0, 1] d ⊂ R d -→ R,
we say that it has a separable representation if

u(x 1 , . . . , x d ) = ∞ j =1 u (j ) 1 (x 1 ) • • • u (j ) d (x d ).
(

Now, consider a mesh of [0, 1] in the x k -variable given by N k -mesh points, 1 ≤ k ≤ d, then we can write a discrete version of (3) by

u(x i 1 , . . . , x i d ) = ∞ j =1 u (j ) 1 (x i 1 ) • • • u (j ) d (x i d ), (4)
where

1 ≤ i k ≤ N k for 1 ≤ k ≤ d. Observe that for each 1 ≤ k ≤ d, if x j k ∈ R N k denotes the vector with components u (j ) k (x i k ) for 1 ≤ i k ≤ N k , then (4) it is equivalent to u = ∞ j =1 x j 1 ⊗ • • • ⊗ x j d .
(

We point out that ( 5) is an useful expression to implemented numerical algorithms using the MATLAB and OCTAVE function kron.

Suppose that for given a linear Partial Differential Equation, and after a discretization by means Finite Elements, we need to solve the linear system

Au = f, (6)
where

A is a (N 1 • • • N d ) × (N 1 • • • N d )-dimensional invert- ible matrix, for some N 1 , . . . , N d ∈ N.
Then from all said above a low rank approximation

A -1 f ≈ u n = n j =1 x j 1 ⊗ • • • ⊗ x j d
with sufficient approximation exists, for some n ≥ 1 and where x j i ∈ R N i for i = 1, 2, . . . , d and j = 1, 2, . . . , n. Moreover, we would to show that

lim n→∞ A -1 f -u n 2 = 0, that is, A -1 f = ∞ j =1 x j 1 ⊗ • • • ⊗ x j d .
Thus, in a first approach to solve it, we would to determine vectors x j 1 , . . . , x j d for j = 1, 2, . . . , n that minimizes

f -A n j =1 x j 1 ⊗ • • • ⊗ x j d 2 ,
or, in short

argmin rank ⊗ u≤n f -Au 2 , ( 7
)
by using the notation introduced in [START_REF] De Silva | Tensor Rank and Ill-posedness of the best low-rank approximation problem[END_REF]. Note that this problem is closely related with the multi-linear generalization of the best low rank approximation of high-order tensors (see also [START_REF] De Lathauwer | A multilinear singular value decomposition[END_REF][START_REF] De Lathauwer | On the best rank-1 and rank-(R 1 , . . . , R N ) approximations of high order tensors[END_REF] and references therein).

For each n ∈ N, we define the set

S n = {x ∈ R N 1 •••N d : rank ⊗ x ≤ n}, introduced in [4], in the following way. Given x ∈ R N 1 •••N d we say that x ∈ S 1 = S 1 (N 1 , N 2 , . . . , N d ) if x = x 1 ⊗ x 2 ⊗ • • • ⊗ x d , where x i ∈ R N i , for i = 1, . . . , d. For n ≥ 2 we define inductively S n = S n (N 1 , N 2 , . . . , N d ) = S n-1 + S 1 , that is, S n = x : x = k i=1 x (i) , x (i) ∈ S 1 for 1 ≤ i ≤ k ≤ n .
Note that S n ⊂ S n+1 for all n ≥ 1. Unfortunately, from Proposition 4.1(a) of [START_REF] De Silva | Tensor Rank and Ill-posedness of the best low-rank approximation problem[END_REF], we have that the set S n is not necessarily (or even usually) closed for each n ≥ 2. However, from Proposition 4.2 of [START_REF] De Silva | Tensor Rank and Ill-posedness of the best low-rank approximation problem[END_REF] it follows that S 1 is a closed set in any norm-topology. This fact implies, as we will see below in Lemma 2, that given an invertible matrix

A ∈ R N 1 N 2 •••N d ×N 1 N 2 •••N d , then for every b ∈ R N 1 •••N d we have argmin x∈S 1 b -Ax 2 = ∅. (8)
This allow to consider the following iterative scheme. Let u 0 = y 0 = 0, and for each n ≥ 1 take

r n-1 = f -Au n-1 , (9)
u n = u n-1 + y n where y n ∈ argmin y∈S 1 r n-1 -Ay 2 . ( 10
)
Note that for each vector

f ∈ R N 1 •••N d and each invert- ible matrix A ∈ R N 1 N 2 •••N d ×N 1 N 2 •••N d ,
we can construct for each n, by using ( 9)-( 10), a vector

u n = n j =1 y n ∈ S n \ S n-1 ,
here we assume that y j = 0 for 1

≤ j ≤ n, that is, rank ⊗ u n = n. Since u n ≈ A -1 f, we define the rank ⊗ for
A -1 f obtained by the Greedy Rank-One Update Algorithm ( 9)-( 10) as

rank G ⊗ (A -1 f) = ∞ if {j ≥ 1 : y j = 0} = ∅, min{j ≥ 1 : y j = 0} -1 otherwise.
The following theorem, which is the main result of this paper, gives the convergence of the Greedy Rank-One Update Approximation for solving linear systems with full rank matrix.

Theorem 1 Let f ∈ R N 1 N 2 •••N d and A ∈ R N 1 N 2 •••N d ×N 1 N 2 •••N d ,
be an invertible matrix. Then, by using the iterative scheme

(9)-(10), we obtain that the sequence { r n 2 } rank G ⊗ (A -1 f) n=0
, is strictly decreasing and

A -1 f = lim n→∞ u n = rank G ⊗ (A -1 f) j =0 y j . ( 11
)
Moreover, the rate of convergence is given by

r n 2 r 0 2 = n j =1 sin θ j (12) for 1 ≤ n ≤ rank G ⊗ (A -1 f) where θ j = arccos r j -1 , Ay j r j -1 2 Ay j 2 ∈ (0, π/2) ( 13
)
for 1 ≤ j ≤ n.
We remark that Grasedyck [START_REF] Grasdyck | Existence and computation of low Kroneckerrank approximations for large linear systems of tensor product structure[END_REF] provides the existence and computation of a low Kronecker rank approximate for a linear system where the matrix A possesses a particular tensor structure and f ∈ S 1 .

From [START_REF] De Lathauwer | On the best rank-1 and rank-(R 1 , . . . , R N ) approximations of high order tensors[END_REF] we obtain that if rank G

⊗ (A -1 f) < ∞, then r n 2 = 0 for all n > rank G ⊗ (A -1 f).
Thus, the above theorem allow to us to construct a procedure, that we give in the pseudo-code form in Algorithm 1, under the assumption that we have a numerical method in order to find a y solving (8) (see the step 5 in Algorithm 1) and that we introduce below.

In order to prove Theorem 1 we need the following two lemmas. From Proposition 4.2 of [START_REF] De Silva | Tensor Rank and Ill-posedness of the best low-rank approximation problem[END_REF] it follows the next result.

Lemma 1 S 1 is a closed set in the • 2 -topology of R N 1 •••N d .
The next lemma give us the existence (but not unicity) of a minimizer for the map

(x) = b -Ax 2 defined from S 1 to [0, ∞). Algorithm 1 Greedy Rank-One Update 1: procedure GROU(f, A, ε, tol, rank_max) 2: r 0 = f 3: u = 0 4:
for i = 0, 1, 2, . . . , rank_max do 5:

y = procedure (min rank ⊗ x≤1 r i -Ay 2 2 )
6: return u and r rank_ max 2 .

r i+1 = r i -Ay 7: u ← u + y 8: if r i+1 2 < ε or | r i+1 2 -r i 2 | <
12: break 13:

return u and r i+1 2 14: end procedure Lemma 2 Assume that A invertible. Then there exists

x * ∈ S 1 such that b -Ax * 2 = min x∈S 1 b -Ax 2 . ( 14
) Proof Let S = { b -Ax 2 : x ∈ S 1 } ⊂ R. Since it is bounded below, there exists γ = inf S = inf x∈S 1 b - Ax 2 ≥ 0. Moreover, for each y 0 ∈ R N 1 •••N d such that y 0 2 = γ, then x 0 = A -1 (b -y 0 ) satisfy that b -Ax 0 2 = γ = inf x∈S 1 b -Ax 2 .
Since γ is the infimum of the set S, for each n > 0 there exists x n ∈ S 1 such that

γ + 1 n ≥ b -Ax n 2 ≥ γ. ( 15
)
Now, let y n = b -Ax n , then for all n > 0,

y n ∈ {y : (γ + 1) ≥ y 2 ≥ γ } ,
which is a compact set in the • 2 -topology. Thus, there exists a subsequence {y n k } such that y n k → y * as n k → ∞ in the • 2 -topology. In consequence,

x n k = A -1 (b -y n k ) ∈ S 1 -→ x * = A -1 (b -y * ) ∈ S 1 , as n k → ∞, because S 1 is a closed set the • 2 -topology.
Finally, substituting in (15) n by n k and taking limits as n k → ∞, we obtain that b -Ax * 2 = γ and this ends the proof of lemma.

Proof of Theorem 1 It is clear that if f = 0 the theorem follows by consider u 0 = y 0 = 0. Thus, from now one we assume that f = 0.

Recall that u 0 = y 0 = 0, and

r n = f -Au n = r n-1 -Ay n ,
for all n ≥ 0. By using Lemma 2, the residual at the n-th step r n satisfy

r n 2 = r n-1 -Ay n 2 = min y∈S 1 r n-1 -Ay 2 . ( 16
)
Since, 0 ∈ S 1 it follows that

r n 2 ≤ r n-1 2 (17)
for all n ≥ 1. From (17) we have that either there exists k ∈ Z + such that

r k+1 2 = r k 2 < r j 2 for j < k, ( 18
) that is, k = rank G ⊗ (A -1 f) < ∞, or r n+1 2 < r n 2 , ( 19
)
for all n ≥ 0. The later condition implies that rank G ⊗ (A -1 f) = ∞. In particular we obtain that the sequence

{r n } rank G ⊗ (A -1 f) n=0 is strictly decreasing. Assume that rank G ⊗ (A -1 f) = k < ∞. Since if x i ∈ R N i for i = 1, 2, . . . , d, then λx 1 ⊗ x 2 ⊗ k • • • ⊗ x d ∈ S 1 for all λ ∈ R. Consequently, we have r k -A(λx 1 ⊗ x 2 ⊗ k • • • ⊗ x d ) 2 ≥ r k 2 (20)
for every λ ∈ R. This implies that r k and A( (24)

x 1 ⊗ x 2 ⊗ k • • • ⊗ x d ) are orthogonal. Since the vectors x i ∈ R N i for i = 1,
By using ( 16) and ( 22)-( 24) we have that for each 1 ≤ n ≤ rank G ⊗ (A -1 f) there exists w n ∈ D, such that

Ay n = r n-1 , w n w n . ( 25
)
Thus,

Ay n 2 = r n-1 , w n > 0,
and

r n-1 , w n ≥ | r n-1 , w |, (26)
for all w ∈ D.

From ( 22)-( 24) and ( 25) we obtain

r n 2 2 = r n-1 -Ay n 2 2 (27) = r n-1 2 2 -| r n-1 , w n | 2 (28) = r n-1 2 2 (1 -ρ 2 n ), ( 29
)
for all 1 ≤ n ≤ rank G ⊗ (A -1 f) and where

ρ n = r n-1 , w n r n-1 2 = r n-1 , Ay n r n-1 2 Ay n 2 = cos θ n ∈ (0, 1).
Proceeding inductively from ( 27)-( 29) follows [START_REF] Mallat | Matching pursuit with time-frequency dictionaries[END_REF]. From (28) we obtain that

r n 2 2 = f 2 2 - n j =1 | r j -1 , w j | 2 (30) for all 1 ≤ n ≤ rank G ⊗ (A -1 f). Since the sequence { r n 2 } ∞ n=0
is decreasing and bounded below, there exists

lim n→∞ r n 2 2 = R 2 = f 2 2 -lim n→∞ n j =1 | r j -1 , w j | 2 .
In consequence,

∞ j =1 | r j -1 , w j | 2 = rank G ⊗ (A -1 f) j =1 | r j -1 , w j | 2
is a convergent series. Thus lim n→∞ | r n-1 , w n | 2 = 0, and from (26) we have that

lim n→∞ | r n-1 , w | 2 = 0 (31) for all w ∈ D. Recall that if rank G ⊗ (A -1 f) < ∞ then r j = 0 for all j ≥ rank G ⊗ (A -1 f). Thus, to end the proof of theo- rem we need to show that if rank G ⊗ (A -1 f) = ∞, then the sequence {r n } ∞ n=0 is convergent in R N 1 N 2 •••N d with the • 2 - topology. Because if it is true, let g = lim n→∞ r n-1 .
From (31), we have g, w = 0 for all w ∈ D. Since S 1 contains a basis of R N 1 N 2 •••N d and A has full rank, then we obtain g = 0. Finally, lim n→∞ r n-1 = 0 and the theorem follows.

Thus, to conclude the proof we only need to show that the residuals sequence {r n } ∞ n=0 is a Cauchy sequence. To prove it the following three lemmas will be useful.

Lemma 3 For each m, n ≥ 1, | r n-1 , Ay m | ≤ r n-1 , w n r m-1 , w m . Proof Since | r n-1 , Ay m | = | r n-1 , r m-1 w m w m | = | r n-1 , w m | r m-1 w m ≤ r n-1 , w n r m-1 w m ,
the lemma follows.

Lemma 4 For every ε > 0 and N ∈ N, there exists τ ≥ N such that

r τ -1 , w τ τ k=1 r k-1 , w k ≤ ε. (32) Proof Since ∞ j =1 r k-1 , w k 2 < ∞, for a given ε > 0 and N ∈ N we choose n ≥ N such that ∞ j =n+1 r k-1 , w k 2 ≤ ε/2.
Since lim k→∞ r k-1 , w k = 0, we can construct a map τ : N -→ N, defined inductively by τ (1) = 1 and

τ (k + 1) = min m≥τ (k) r m-1 , w m ≤ r τ (k)-1 , w τ (k)
for all k ≥ 1, such that τ is strictly increasing and lim k→∞ τ (k) = ∞. Observe that for all k ≥ 1 and j satisfying that τ(k) < j < τ(k + 1) it follows

r τ (k+1)-1 , w τ (k+1) ≤ r τ (k)-1 , w τ (k) ≤ r j -1 , w j . Thus r τ (k+1)-1 , w τ (k+1) ≤ r j -1 , w j ,
for all 1 ≤ j < τ(k + 1). Now, we can choose τ = τ (k + 1) > n large enough, satisfying that

r τ -1 , w τ n k=1 r k-1 , w k ≤ ε/2, because lim k→∞ r τ (k)-1 , w τ (k) = 0. Then r τ -1 , w τ τ k=1 r k-1 , w k = r τ -1 , w τ n k=1 r k-1 , w k + r τ -1 , w τ τ k=n+1 r k-1 , w k ≤ ε/2 + τ k=n+1 r τ -1 , w τ r k-1 , w k ≤ ε/2 + τ k=n+1 r k-1 , w k 2 ≤ ε/2 + ∞ k=n+1 r k-1 , w k 2 ≤ ε.
This proves the lemma.

Lemma 5 For each M > N > 0 it follows that r N-1 -r M-1 2 2 ≤ r N-1 2 2 -r M-1 2 2 + r M-1 , w M M k=1 r k-1 , w k . Proof Since r N-1 = r M-1 + M-1 k=N Ay k we have r N-1 -r M-1 2 2 = r N-1 2 2 + r M-1 2 2 -2 r M-1 + M-1 k=N Ay k , r M-1 = r N-1 2 2 -r M-1 2 2 -2 M-1 k=N Ay k , r M-1 ≤ r N-1 2 2 -r M-1 2 2 + 2 M-1 k=N r M-1 , w M r k-1 , w k by Lemma 3, ≤ r N-1 2 2 -r M-1 2 2 + 2 M k=1 r M-1 , w M r k-1 , w k
by adding positive terms.

This ends the proof of the lemma.

Since lim n→∞ r n 2 2 = R 2 , and it is a decreasing sequence, for a given ε > 0 there exists k ε > 0 such that

R 2 ≤ r m-1 2 2 ≤ R 2 + ε 2 /8
for all m ≥ k ε . Assume that m > k ε . From Lemma 4, for each m + p there exists τ > m + p such that

r τ -1 , w τ τ k=1 r k-1 , w k ≤ ε 2 /8.

Now, we would to estimate

r m-1 -r m+p-1 2 ≤ r m-1 -r τ -1 2 + r τ -1 -r m+p-1 2 .
By using Lemma 5 with M = τ and N = m and m + p, respectively, we obtain that

r m-1 -r τ -1 2 2 ≤ R 2 + ε 2 /8 -R 2 + ε 2 /8 = ε 2 /4,
and

r m+p-1 -r τ -1 2 2 ≤ R 2 + ε 2 /8 -R 2 + ε 2 /8 = ε 2 /4,
respectively. In consequence {r n } ∞ n=0 is a Cauchy sequence and it converges to 0.

A Block Coordinated Descent Approach for the Rank-One Minimization Problem

In this section we study the Rank-One minimization problem min

x∈S 1 b -Ax 2 , (33)
that we can write as the following unconstrained optimization problem:

min (x 1 ,...,x d )∈R N 1 +•••+N d b -A (x 1 ⊗ • • • ⊗ x d ) 2 . ( 34
)
A popular method for minimizing a real-valued continuously differentiable function of N 1 + • • • + N d real variables, subject to bound constrains, is the (block) coordinate descend method. In this method, the coordinates are partitioned into N k blocks and, at each iteration, b is minimized with respect to one of the coordinate blocks while the others are held fixed (see Algorithm 2). These cyclic methods have the advantage of not requiring any information about the gradient to determine the descent directions. However, their convergence properties are poorer than steepest descend methods. Moreover, its are attractive because of their

Algorithm 2 A Block Coordinated Descent Algorithm 1: Given : R N 1 +•••+N d → R 2: Initialize x 0 i ∈ R N i for i = 1, 2 . . . , d. 3: for n = 1, 2, . . . do 4:
for k = 1, 2, . . . , d do 5:

x n k ∈ arg min x k ∈R N k (x n 1 , . . . , x n k-1 , x k , x n-1 k+1 , . . . , x n-1 d ) 6:
end for 7: end for easy implementation in some particular cases as we will see below.

The main result of this section is the following.

Theorem 2 Let b ∈ R N 1 N 2 •••N d and A ∈ R N 1 N 2 •••N d ×N 1 N 2 •••N d , be an invertible matrix. Assume that for each k ∈ {1, 2, . . . , d} the (N 1 • • • N d ) × N k -matrix Z k = A x 1 ⊗ • • • ⊗ x k-1 ⊗ I N k ⊗ x k+1 • • • ⊗ x d has linearly independent columns for every (x 1 , . . . , x d ) ∈ R N 1 +•••+N d satisfying b -A (x 1 ⊗ • • • ⊗ x d ) 2 ≤ b -A(x 0 1 ⊗ • • • ⊗ x 0 d ) 2 . ( 35
)
Then every accumulation point (x * 1 , . . . , x * d ) of the sequence

x n 1 , . . . , x n d ∞ n=0 ,
generated by Algorithm 2 using the map

(x 1 , . . . , x d ) = b -A(x 1 ⊗ • • • ⊗ x d ) 2 ,
satisfies the equation ∇ (x * 1 , . . . , x * d ) = 0. Moreover, assume that x 1 , . . . , x k-1 , x k+1 , . . . , x d , are fixed for some k ∈ {1, 2, . . . , d}, then

x k = (Z T k Z k ) -1 Z T k b,
is the global minimum of the directional minimization problem

min x∈R N k b -A (x 1 ⊗ • • • ⊗ x k-1 ⊗ x ⊗ x k+1 ⊗ • • • ⊗ x d ) 2 . ( 36
)
Proof To prove the theorem we will use the following lemma, that it can be proved by using the same argument as the proof of Theorem 5.32 in [START_REF] Ruszczyński | Nonlinear optimization[END_REF].

Lemma 6 Assume that the function : R N 1 +•••+N d → R is continuously differentiable and that the set Since

X 1 = (x 1 , . . . , x d ) ∈ R N 1 +•••+N d : (x 1 , . . . , x d ) ≤ (x 0 1 , . . . , x 0
x 1 ⊗ • • • ⊗ x k-1 ⊗ x ⊗ x k+1 • • • ⊗ x d = (x 1 ⊗ • • • ⊗ x k-1 ⊗ I k ⊗ x k+1 • • • ⊗ x d ) × (1 ⊗ • • • ⊗ 1 ⊗ x ⊗ 1 • • • ⊗ 1) = (x 1 ⊗ • • • ⊗ x k-1 ⊗ I k ⊗ x k+1 • • • ⊗ x d ) x, we can write b -A (x 1 ⊗ • • • ⊗ x d ) 2 = b -Z k x k 2 . ( 37
)
Then, for a fixed x 1 , . . . , x k-1 , x k+1 , . . . , x d the directional minimization problem (36) is equivalent to the standard Least Squares problem min

x∈R N k b -Z k x 2 . ( 38
)
Since if Z k has linearly independent columns then (38) has a unique solution

x k = (Z T k Z k ) -1 Z T k b,
from Lemma 6, the theorem follows.

Note that given a point (x 1 , . . . , x d ), descend with respect to the coordinate x k means, from Theorem 2, that we need to solve the standard least squares problem (38). In particular we minimize cyclically with respect to the coordinate variables Thus, Theorem 2 allow to us to solve the rank-one minimization problem (33) by means the Alternate Least Squares (ALS) Algorithm 3. We point out that Algorithm 3 ALS Algorithm 1: procedure ALS(b, A, iter_max, tol)

2:
Initialize x 0 i for i = 1, 2, . . . , d.

3:

iter = 1 4:

while iter < iter_max do 5: xk ← x 0 k , k = 1, . . . , d 6:
for k = 1, 2, . . . , d do 7:

Z = A(x 0 1 ⊗ • • • ⊗ x 0 k-1 ⊗ I N k ⊗ xk+1 ⊗ • • • ⊗ xd ) 8: x 0 k = (Z T Z) -1 Z T b 9:
end for 10: return x 0 = (x 0 1 , . . . , x 0 d ) 15: end procedure for high-dimensional problems the numerical implementation of solving the equation

if d k=1 x 0 k -xk 2 <
Z T k Z k x = Z T k x, (39)
can be a hardly task. However, if the matrix A can be represented also in separated representation form, then as the following corollary shows (39) can be implemented in a more easy way by using the properties of the Kronecker product.

Corollary 1 Assume that

A = r A j =1 A j 1 ⊗ • • • ⊗ A j d ,
where

A j i ∈ R N i ×N i for i = 1, 2, . . . , d and j = 1, 2, . . . , r A . Let k ∈ {1, 2, . . . , d} and assume that for x 1 , . . . , x k-1 , x k+1 , . . . , x d fixed, the (N 1 • • • N d ) × N k -matrix Z k = r A j =1 A (j ) 1 x 1 ⊗ • • • ⊗ A (j ) k-1 x k-1 ⊗ A (j ) k ⊗ A (j ) k+1 x k+1 • • • ⊗ A (j ) d x d ,
has linearly independent columns. Then

x * k = (Z T k Z k ) -1 Z T k b (40)
is the global minimum of the directional minimization problem (36).

This corollary implies that we can solve the minimization problem min

(x 1 ,...,x d ) b - r A i=1 A i 1 x 1 ⊗ • • • ⊗ A i d x d 2 , (41)
by means the ALS Algorithm 4.

Numerical Examples

In order to illustrate the Greedy Rank-One Update Algorithm 1 by using the Alternate Least Squares Algorithms 3 and 4 we give the following examples. First, we compute for a given vector f ∈ R 14 3 its approximation by a vector

n j =1 f 1 j ⊗ f 2 j ⊗ f 3 j .
To this end we consider

A = I ⊗ I ⊗ I
where I is the identity matrix of size 14 × 14. We take as a choice of f a random vector. Then for the parameter values iter_max = 10, rank_max = 1000 tol = 2.22e -16, and ε = 1.0e -08, we compute the sequences {θ j } n j =1 , Fig. 1 The sequences {θ j } 1000 j =1 , and {log 10 r j 2 } 1000 j =0 when A is the identity matrix Algorithm 4 ALS with A is separable form

1: procedure ALS(b, r A i=1 A i 1 ⊗ • • • ⊗ A i d , iter_max , tol) 2:
Initialize x 0 i for i = 1, 2 . . . , d.

3:

iter = 1 4:

while iter < iter_max do 5: xk ← x 0 k , k = 1, . . . , d 6:
for k = 1, 2, . . . , d do 7:

Z = r A j =1 A (j ) 1 x 0 1 ⊗ • • • ⊗ A (j ) k-1 x 0 k-1 ⊗ A (j ) k ⊗ A (j ) k+1 xk+1 • • • ⊗ A (j ) d xd 8: x 0 k = (Z T Z) -1 Z T b 9:
end for 10: return x 0 = (x 0 1 , . . . , x 0 d ) 15: end procedure given by [START_REF] Ruszczyński | Nonlinear optimization[END_REF], and { r j 2 } n j =0 . As can be seen in Fig. 1, the first one is an oscillating sequence around approximately 1.39, and, as Theorem 1 shows, the norm residuals is strictly decreasing. The algorithm stopped at n = 1000 after 8.9 seconds with a relative error equal to 2.86238020895006e -08. Now, in our second example we would to study the relationship between the sequence of angles {θ j } n j =1 ⊂ (0, π/2), given in the statement of Theorem1, and the sequence of the norm of the residuals { r j 2 } n j =1 . Here n is the rank G ⊗ u where u is a numerical approximation of A -1 f. Motivated by the above example we consider the mean angle, denoted by

if d k=1 x 0 k -xk 2 <
θ n = 1 n n j =1 θ j ,
of a finite sequence {θ j } n j =1 . Now, we perform the following numerical experiment. We fix the values N 1 = 3, N 2 = 4, tol = 2.22e -16, ε = 1.0e -08, iter_max = 10 and rank_max = 1000. We repeat the following procedure 500-times. First, we generated randomly a matrix A ∈ R N 1 N 2 ×N 1 N 2 , and a solution A -1 f, in order to construct the right-side hand f of the linear system. Finally, we run Algorithm 1 to obtain an approximate solution u with rank G ⊗ u and a mean angle θ rank G ⊗ u . In Fig. 2 we can see a sample obtained in this experiment, for this particular case the approximate solu-Fig. 2 The sequences {θ j } 62 j =1 , and {log 10 r j 2 } 62 j =0 when we take as a choice of A a 20 × 20 random matrix tion u had rank G ⊗ u = 62, with a relative error equal to 9.47228661927500e -08. Moreover, the matrix A was nearly singular with a condition number equal to 488. Note, that the angle sequence, as we can see in Fig. 2, oscillates around θ rank G ⊗ u = 0.949051119522545. In Fig. 3 we plot the mean angle as a function of the rank G ⊗ of the approximate solution. From this picture we can deduce that some weakly linear dependence exists. In this section we present the following example. Recall that the Poisson problem reads

-u = f in ⊂ R d u| ∂ = 0. ( 42
)
where f = f (x 1 , x 2 , . . . , x d ) is a given function and

= d i=1 ∂ 2 ∂x 2 i
is the Laplace operator. In order to find its variational formulation, we recall the following Green formula for the Laplacian:

- uv dx = ∇u • ∇v dx - ∂ ∂u ∂n v dγ. ( 43
)
Assume that = (0, 1) d and for d = 1, 2, . . . let be the bilinear form

a d (u; v) = ∇u • ∇v dx.
Take d = 3 in (42) and then we can easily deduce that u satisfies the following problem: Find u ∈ H 1 0 ( ) such that

a 3 (u; v) = f v dx for all v ∈ H 1 0 ( ). ( 44
)
The Galerkin approximation to (44) reads:

find u h ∈ V h : a 3 (u h ; v h ) = f v h dx ∀v h ∈ V h . ( 45
)
Assume V h = P ⊗ P ⊗ P where P = span{ω 1 , . . . , ω N } and ω 1 , . . . , ω N in H 1 0 (0, 1) are the following N -linearly independent maps. First, we partitioned the interval [0, 1] into N -parts Denote h i = x i+1x i and h = max 1≤i≤N h i . For i = 1, 2, . . . , N -1, let

0 = x 1 < x 2 < • • • < x N+1 = 1.
ω (i) (x) = ⎧ ⎪ ⎨ ⎪ ⎩ (x -x i )/ h i x i ≤ x ≤ x i+1 , ( x i+2 -x)/ h i+1 x i+1 ≤ x ≤ x i+2 , 0 otherwise.
These functions are continuous and piecewise linear. It is easy to see that they are linearly independent. The first order weak derivatives of the basis functions are piecewise constant. Indeed for i = 1, 2, . . . , N -1

d dx ω (i) (x) = ⎧ ⎪ ⎨ ⎪ ⎩ 1/h i x i ≤ x ≤ x i+1 , -1/h i+1 x i+1 ≤ x ≤ x i+2 , 0 otherwise.
We assume that we have an uniform partition, that is, h i = h for i = 1, 2, . . . , N. Then the following formulas are useful

1 0 d dx ω (i) (x) d dx ω (i-1) (x) dx = - 1 h for 2 ≤ i ≤ N -1, ( 46
) 1 0 d dx ω (i) (x) 2 dx = 2 h for 1 ≤ i ≤ N -1, ( 47
) 1 0 ω (i) (x)ω (i-1) (x) dx = h 6 for 2 ≤ i ≤ N -1, ( 48
) 1 0 ω (i) (x) 2 dx = 2h 3 for 1 ≤ i ≤ N -1. (49)
Now, let the stiffness multilinear matrix A 3 defined by

(A 3 ) j 1 ,j 2 ,j 3 ;i 1 ,i 2 ,i 3 = a 3 (ω (i 1 ) ⊗ ω (i 2 ) ⊗ ω (i 3 ) ; ω (j 1 ) ⊗ ω (j 2 ) ⊗ ω (j 3 ) ) = ∂ω (i 1 ) ∂x 1 ∂ω (j 1 ) ∂x 1 ω (i 2 ) ω (j 2 ) ω (i 3 ) ω (j 3 ) + ∂ω (i 2 ) ∂x 2 ∂ω (j 2 ) ∂x 2 ω (i 1 ) ω (j 1 ) ω (i 3 ) ω (j 3 ) × ∂ω (i 3 ) ∂x 2 ∂ω (j 3 ) ∂x 2 ω (i 1 ) ω (j 1 ) ω (i 2 ) ω (j 2 ) dx = A j 1 ;i 1 B j 2 ;i 2 B j 3 ;i 3 + B j 1 ;i 1 A j 2 ;i 2 B j 1 ;i 1 + B j 1 ;i 1 B j 2 ;i 2 A j 3 ;i 3 ,
where A j,i = a 1 (ω (i) ; ω (j ) ) and B j ;i =

1 0 ω (i) (x)ω (j ) (x)dx.
Then, by using ( 46)-(49), we have that

A = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 2 h -1 h -1 h 2 h -1 h . . . . . . . . . -1 h 2 h -1 h -1 h 2 h ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ and B = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 2h 3 h 6 h 6 2h 3 h 6 . . . . . . . . . h 6 2h 3 h 6 h 6 2h 3 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . Fig. 4 The relative error u 1 -A -1 f 2 / A -1 f 2 in logarithmic scale
Finally, by using (2), we obtain that

(A 3 ) j 1 ,j 2 ,j 3 ;i 1 ,i 2 ,i 3 = (A ⊗ B ⊗ B + B ⊗ A ⊗ B + B ⊗ B ⊗ A) s,t .
On the other hand we have that f (ω (j 1 ) ⊗ ω (j 2 ) ⊗ ω (j 3 ) ) dx = F j 1 ,j 2 ,j 3 .

In order to solve (45) we consider a map in V h given by

u h = i 1 ,i 2 ,i 3 U i 1 ,i 2 ,i 3 ω (i 1 ) ⊗ ω (i 2 ) ⊗ ω (j 3 ) , where U i 1 ,i 2 ,i 3 ≈ u( x i 1 +1 , x i 2 +1 , x i 3 +1 ). Thus, U k,i 2 ,i 3 = U i 1 ,k,i 3 = U i 1 ,i 2 ,k = 0 for k ∈ {0, N}. Substitute u h in (45) with v h = ω (j 1 ) ⊗ ω (j 2 ) ⊗ ω (j 3 ) . Then we obtain that N-1 i 1 ,i 2 ,i 3 =1 (A 3 ) j 1 ,j 2 ,j 3 ;i 1 ,i 2 .i 3 U i 1 ,i 2 ,i 3 = F j 1 ,j 2 ,j 3 (50)
must be hold for all j 1 , j 2 , j 3 ∈ {1, . . . , N -1}. This implies that (45) is equivalent to solve the following linear system

(A ⊗ B ⊗ B + B ⊗ A ⊗ B + B ⊗ B ⊗ A)u = f, (51)
here we consider that a general multi-index tensor can be represented by a standard vector by using that

V j 1 ,j 2 ,...,j d = v s if and only if s = j d + d-1 l=1 (j l -1) d l =l+1 N l .
For a general d ≥ 2 it can be shown that in order to solve numerically (42) we need to solve the following linear system:

d j =1 A (j ) 1 ⊗ • • • ⊗ A (j ) d u = f, (52)
where

A (j ) k = A if k = j, B if k = j.
Finally, it is not difficult to show that A and B are symmetric and definite positive matrices. Then, from the properties 7-8 listed at the end of Sect. 1, d j =1 A (j )

1 ⊗ • • • ⊗ A (j )
d is a symmetric and definite positive matrix. In consequence, it is invertible.

Some Numerical Examples

Example 1 Firstly, we consider the following problem in 3D: Solve for (x 1 , x 2 , x 3 ) ∈ = (0, 1) 3 : We used the separable representation algorithm given in Sect. 2 with parameter values iter_max = 5, rank_max = 1000 and ε = 0.001. The algorithm give us an approximated solution u 1 ∈ S 1 . In Fig. 4 we represent the relative error of the solution computed using the separable representation algorithm, using logarithmic scale, as a function of the number of nodes used in the discretization of the Poisson equation. All the computations were performed using the GNU software OCTAVE in a AMD 64 Athlon K8 with 2Gib of RAM. In Fig. 5 we represent the CPU time, in logarithmic scale, used in solving the linear system (52) against the separable representation algorithm. In both cases all the linear systems involved were solved using the standard linear system solver (A\b) of OCTAVE.

-u = (2π) 2 • 3 • sin(2πx 1 -π) × sin(2πx 2 -π) sin(2πx 3 -π), (53)
Example 2 Finally we are addressing some highly multidimensional models. To this end we solve numerically (42) for Here we consider the true solution u given by U i 1 ,...,i d = u( x i 1 +1 , . . . , x i d +1 ). For d = 10 we use the parameter values iter_max = 2, rank_max = 10 and ε = 0.001. In a similar way as above the algorithm give us an approximated solution u ∈ S 1 . In Fig. 6 we represent the absolute error uu 2 as a function of h = π/N for N = 5, 10, 20, . . . , 160 in log 10 -scale. By using similar parame-ters values the problem has been solved for d = 100 in about 20 minutes.

Concluding Remarks

In this paper we prove the convergence of the Greedy Rank-One Update Algorithm for solving linear systems with a full rank matrix. Moreover, we study the Rank-One Minimization Problem and show that a Block Cyclic Coordinate Descend strategy implies an Alternating Least Squares Algorithm. As we can show the method runs under very weak conditions, recall that we only use the assumption that the linear system has a an invertible matrix. However, its efficiency depends strongly on the matrix form (symmetric, tridiagonal, full, sparse, . . .).

  A and B are banded, then A ⊗ B is banded. 7. If A and B are symmetric, then A ⊗ B is symmetric. 8. If A and B are definite positive, then A ⊗ B is definite positive.

1 b 1 b

 11 2, . . . , d, are arbitrary and A has full rank, this means that r k = 0 and in consequence u k = A -1 f and r j = 0 for all j ≥ k. Let be D = {w : w = Ay, w 2 = 1 and y ∈ S 1 } , then min y∈S -Ax 2 = min w∈D,λ∈R bλw 2 . Note that bλw 2 2 = b 2 2 -2λ b, w + λ 2 , (21) then the minimum for bλw 2 2 when w ∈ D is obtained for λ = b, w > 0 and bλw 2 2 = b 2 2 -| b, w | 2 . Thus, min y∈S -Ay 2 2 = min λ∈R,w∈D b -

4. 1 A

 1 Model Problem: The Poisson Equation in (0, 1) d

Fig. 3

 3 Fig. 3 The mean angle θ rank G ⊗ u as a function of rank G ⊗ u. The discontinuous line is the least squares straight line y = (7.69351915752358e -05)x + 1.36946939179450

Fig. 5

 5 Fig.[START_REF] Graham | Kronecker products and matrix calculus with applications[END_REF] The CPU time, in seconds, used in solving the linear system as a function of the number of nodes employed in the discretization of the Poisson Equation

(x 1

 1 , . . . , x d ) ∈ = (0, π) d where f = d k=1 -(1 + k) sin (-1+k) (x k )(-k cos 2 (x k ) + sin 2 (x k )) d k =1,k =k sin (1+k ) (x k ),

Fig. 6

 6 Fig.[START_REF] Grasdyck | Existence and computation of low Kroneckerrank approximations for large linear systems of tensor product structure[END_REF] The absolute error uu 2 as a function of h = π/N for N = 5, 10, 20, . . . , 160 in log 10 -scale

tol then goto 14

	11:	end if
	12:	iter = iter + 1
	13:	end while
	14:


Audrey Bürki¹²
email: audrey.buerki@pse.unige.ch

Cécile Fougeron¹

Cédric Gendrot¹
email: cedric.gendrot@univ-paris3.fr

On the Categorical Nature of the Process Involved in Schwa Elision in French

Keywords: schwa elision, french, categorical process

This paper examines the nature of the process involved in optional schwa elision in French. More specifically, it aims at testing whether this process is gradual or categorical, on the basis of an analysis of the distribution of the duration of over 4000 schwas extracted from a large corpus of continuous speech. The distribution observed is bimodal, with absent schwas (at 0 ms duration) on one side, and realized schwas on the other side, the two groups being separated by a small gap in the distribution. Realized schwas present cases of strong temporal reduction, but this reduction does not show a continuous pattern toward zero duration, as would be predicted if schwa elision was the end-point of a gradual reduction process. Different interpretations of this distribution are discussed.

Introduction

Schwa elision in French is defined as the alternation between two variants of a same lexical entry containing or not a schwa. For example, the word 'fenêtre' (window) can be equally pronounced as a schwa-full form [fђnϯtల] or as a schwa-less form [fnϯtల]. Various factors are known to affect the realization of one or the other variant (speech style, speaker age and dialect, position in word, etc…), although the choice for one or the other variant is not fully predictable.

The alternation between schwa-less and schwa-full forms has been interpreted as either a vowel deletion process, or a vowel epenthesis process (see [START_REF] Côté | The nature of the schwa/zero alternation in French clitics: experimental and non-experimental evidence[END_REF] for other accounts). Be it a deletion or an epenthesis, traditional phonological accounts assume that the process is complete, and results from a rulebased segmental process.

As part of a larger project on the nature of schwa and schwa elision in French, the present study aims at determining whether experimental evidence can be found to assess the categorical nature of this process. Only few acoustic or articulatory studies have addressed this question. They have usually focused on the phonetic properties of the consonant sequences resulting from the elision of schwa, asking whether some residues attributable to schwa remain in cases where the vowel is apparently absent. Consonant sequences resulting from schwa elision are thus compared to similar consonant sequences not resulting from elision. Results usually agree on the fact that the former differ from the latter. However, while these acoustic or articulatory differences are always taken to reflect remnants of the vowel, their existence is not always interpreted in the same way.

In [START_REF] Rialland | Schwa et Syllabes en Français[END_REF] for example, the observed acoustic differences between sequences of the type 'le bar trouvé' vs 'le bas r'trouvé' (the apostrophe indicating the elision, the bar found vs the sock found again) are attributed to the fact that in the second case the sequence retains its phonological nucleus despite the loss of its segmental content. A palatographic study, [START_REF] Fougeron | Au Delà de la Syllabe: le Rôle des Informations Articulatoires Stockées dans le Lexique pour l'Analyse de la Chute de Schwa[END_REF] found articulatory differences when comparing consonant sequences resulting from schwa elision (e.g. 'pas d'rôle', no role) to underlying onset clusters (e.g. 'pas drôle', not funny) and consonant sequences spanning a word boundary (e.g. 'jade rose', rose jade). They do not take these results as evidence for the maintenance of the syllable, but rather interpret them as a preservation of the pre-vocalic underlying articulatory specification of the consonants, despite the deletion of schwa. In her formalization of these results, [START_REF] Steriade | Paradigm Uniformity and the Phonetics-Phonology Boundary[END_REF] considers that these characteristics of the surrounding consonants, inherited from the form with schwa, are extended to the form where the schwa is categorically elided following Paradigm Uniformity principles.

While these interpretations ( [START_REF] Rialland | Schwa et Syllabes en Français[END_REF][START_REF] Fougeron | Au Delà de la Syllabe: le Rôle des Informations Articulatoires Stockées dans le Lexique pour l'Analyse de la Chute de Schwa[END_REF][START_REF] Steriade | Paradigm Uniformity and the Phonetics-Phonology Boundary[END_REF]) have different implications, none of them raises doubts about the fact that the absence of schwa results from a categorical segmental process. This is not the case of other studies, which consider the phonetic differences observed in [START_REF] Rialland | Schwa et Syllabes en Français[END_REF] and [START_REF] Fougeron | Au Delà de la Syllabe: le Rôle des Informations Articulatoires Stockées dans le Lexique pour l'Analyse de la Chute de Schwa[END_REF] as evidence for a gradual reduction process. According to the analysis proposed for English schwa elision within the Articulatory Phonology framework [START_REF] Browman | Targetless" Schwa: an Articulatory Analysis[END_REF], schwa elision in French would be the consequence of a gradual process of overlap between gestures. According to this view, in apparent schwa-less French forms, the vowel is not acoustically or perceptualy visible, not because of segmental deletion, but because it is hidden by the overlap of the surrounding consonantal gestures. So-called schwa elision would thus be the end point of this gradual reduction mechanism. In that sense, it is similar to the other coarticulatory or allophonic variations, explained by the overlap of invariant underlying gestures. Support for this overlap account in French is found in the results of [START_REF] Smorodinsky | The phonology and phonetics of schwas in parisian French: an articulatory analysis[END_REF] who observed some lingual gestures attributable to schwa in apparent schwa-less forms, and in the results of [START_REF] Barnes | Phonetic analogy and schwa deletion in French[END_REF] who observed some remnants of the rounded articulation of the elided schwa on the surrounding consonants.

Although these results have received contradictory explanations or have failed to be replicated (see [START_REF] Côté | The nature of the schwa/zero alternation in French clitics: experimental and non-experimental evidence[END_REF]), their major shortcoming resides in their assumption that schwa elision in French is a 'casual speech' or 'fast speech' process. Their argument is that in casual or fast speech, overlap is known to increase, and thus schwa has more chance to be fully overlapped (thus hidden). However, even though it is well agreed upon that there are more schwa elisions in casual and fast speech style in French, schwa-less variants do occur frequently at slow and formal speech styles. It cannot thus be reduced to a continuous speech process.

In a recent paper, [START_REF] Davidson | Schwa elision in fast speech: segmental deletion or gestural overlap?[END_REF] examined the acoustic properties of consonant sequences resulting from the elision of schwa in English. Her results confirm the existence of acoustic residues attributable to schwa in apparent schwa-less forms. Following [START_REF] Browman | Targetless" Schwa: an Articulatory Analysis[END_REF], she also argued that these residues are better accounted for by a change in gestural overlap, rather than segmental deletion. More interestingly, although she manipulated speech rate to trigger schwa elision, she observed that some of her speakers also produced elision at normal rates. This leads her to interpret elision at fast and at normal rate by the same process of overlap, a categorical shift toward a more overlapping coordination pattern between gestures. Hence, vowel elision, although explained by overlap, would not only result from a gradual reduction process, but could be the consequence of a categorical change in the specification of the coordination of the gestures. According to Articulatory Phonology principles, this specification is assumed to be assigned by the phonological system. While this hypothesis has to be further elaborated to be applied to schwa alternation in French, it reconciles the views that elision may be due to overlap and that this process would be categorical.

This review of the literature shows that determining whether schwa elision is a gradual or categorical process in French is not straightforward. In this study, our goal is to question the nature of the schwa-zero alternating process, but our focus is not on the potential acoustic or articulatory residues of schwa in apparent schwa-less forms. Our approach is different. Rather than focussing on cases where schwa is acoustically and perceptively absent, as done in the studies presented above, we concentrate our attention on the different realizations of words containing a schwa (in their schwa-less and schwa-full forms). This is done by observing the production in continuous speech of a large variety of words uttered by a large number of different speakers. Our hypothesis is that variations in the realization of these words can inform us about the nature of the process. More specifically, we test whether the realizations of this vowel follow a continuous pattern of reduction ranging from a full vowel realization to a fully absent vowel (i.e. the end point of a reduction continuum), or whether two categories emerge: cases with full vowels on one side, and with completely deleted vowels on the other side. Support for a categorical process would be found in a bimodal distribution with a gap in the distribution between the 2 groups. Support for a gradual process would appear in a continuous distribution (without a gap).

Method

A repertory of 18553 French words containing a schwa (obligatory or optional, not occurring at word or clitic boundaries or in compounds) was created through the compilation of several lexical databases (Lexique [START_REF] New | Une base de données lexicales du français contemporain sur internet : LEXIQUE[END_REF], Brulex [START_REF] Content | Brulex, une base de données lexicales informatisée pour le français écrit et parlé[END_REF], IlPho [START_REF] Boula De Mareüil | A French Phonetic Lexicon with variants for Speech and Language Processing[END_REF], Grand Robert dictionary). In order to observe the realization of a large number of lexical forms produced by many speakers, we selected the ESTER corpus containing 24h of radio broadcasted news including productions of a total of 574 speakers [START_REF] Galliano | ESTER Phase II Evaluation Campaign for the Rich Transcription of French Broadcast News[END_REF]. We decided to discard schwas appearing before derivational suffixes or inflectional endings.

Our focus was on schwas that can alternate with ∅. Two types of schwa are described in French: optional (alternating) schwas and obligatory schwas (always produced). While some words contain unambiguous obligatory schwas (e.g. 'belette', weasel), for other it is in practice very difficult to assert if the schwa is optional or not. A first intuitive categorization of a subset of the corpus according to whether the word could be 'heard' without its schwa, revealed a great disparity in the judgements of the 3 authors. Thus, in order to select our test items, we decided to choose a more objective method based on the analysis of the actual realization of these words in our speech database. On the basis of the transcription given by the automatic recognition system developed at IRISA, we selected the words that were recognized by the system as being produced in the corpus sometimes in their schwa-less form, sometimes in their schwa-full form. After elimination of the production of 22 speakers having a strong foreign accent, and of 95 occurrences showing a clear recognition error or having a bad acoustic quality, 4320 occurrences (corresponding to 203 different words) were selected for analysis.

A manual examination of the output of automatic recognition system was undertaken for all the sound files in order to correct schwa detection errors (wrong omission / insertion). Since the shortest segment duration assigned by the system is 30 ms, and since the automatic alignments are far from being reliable with respect to the precision required in this study, duration measurements could not depend on this automatic segmentation. Automatic segmentations / alignments were thus corrected manually, according to the following criteria. We considered a schwa as being realized if it presented a voiced portion and a formant structure. Fully devoiced vowels were thus not considered as schwas. The onset and offset of the vowels were determined based on the signal periodicity and on apparition / disappearance of the second formant, as illustrated in the segmentation shown in Figure 3. These criteria for segmentation of the vowel intervals were the best compromise with regard to the wide variety of consonantal contexts. As can be seen in the example of 'semaine' in Figure 3, the duration of the segment corresponding to our definition is much shorter than if we had included voiceless portions of the vowel. The consequences of this type of segmentation are discussed in section 4. The duration of the schwa intervals was then automatically computed with Praat [START_REF] Boersma | Praat: doing phonetics by computer[END_REF].

In order to avoid potential confounds due to factors known to affect the frequency of schwa elision, the words present in the corpus were coded according to both their syllabic position in the word, and their location relative to a morpheme boundary. Schwas in word-initial syllables are known to be less frequently elided than non-initial schwa, while morpheme-internal schwas are more consensually considered as underlying lexical schwas compared to those located at morpheme boundaries [see[ [START_REF] Côté | Consonant Cluster Phonotactics: a Perceptual Approach[END_REF]].

Results

General distribution of duration

The distribution of the duration of the 4185 optional schwas of our corpus is illustrated in Figure 1. This distribution is divided into two sub-groups. It excludes 4% of the occurrences, whose duration was over 2 standard-deviations away from the mean. On one side, at 0 ms, 1215 occurrences (29%) of the words containing an optional schwa were produced with a complete absence of voicing and formant structure. According to our criteria, these forms are considered as schwa-less. As an illustration of these cases, a schwa-less occurrence of the word 'sera' (will be) is presented in figure 2. On the other side, between 16 and 86 ms, 71% of the occurrences present a vocalic portion corresponding to a schwa. Confirming the Schwa duration Nb of occurrences literature, these 2970 realized schwas in our corpus are relatively short vowels, with a mean duration of 49 ms (σ = 15). The distribution of the duration of these schwa-full forms is unimodal with close to zero skewness (0.13) and kurtosis (-0.52) coefficients. Figure 3 shows two examples of schwa-full occurrences of the words 'sera' (will be) and 'semaine' (week).

More interestingly, with regard to our acoustic criteria, no schwas were observed with a duration of less than 16 ms, thus a gap is observed in the distribution beween 0 and 16 ms. This gap in the duration distribution between the 2 subgroups is in favour of a categorical process. It has to be noted that the outliers we removed on the left-hand side of the distribution represent only 0.7 % of the schwa-full forms distribution. This removal cannot therefore be responsible for the gap we observe.

It also has to be mentionned that the distribution is similar when the words having a schwa in an /R/ context (words difficult to segment) are excluded previously from the analysis.

Figure 2: Schwa-less occurrences of the word 'sera', will be (note that the /R/ is devoiced) and 'appelé', called.

Distribution of the schwa duration according to its position in word and in morpheme

Out of our 4185 occurrences, 3445 words have an optional schwa in the word-initial syllable (e.g. 'regarde' look, 'devient' become). For these words, the distribution is similar to the one presented above. Two groups emerge with a gap in the duration distribution between 0 and 16 ms. The mean schwa duration of the 2815 schwa-full forms is 50 ms (σ = 15), with close to zero skewness (0.14) and kurtosis (-0.51) coefficients. This distribution is therefore close to normality while right-skewed.

The remaining 740 words have a non-initial schwa, placed either in the second or third syllable of the word. Their distribution is similar to that of initial schwas with 2 separate sub-groups and a gap in the distribution between 0 and 17 ms.

According to what is known in the literature [e.g. [START_REF] Côté | Consonant Cluster Phonotactics: a Perceptual Approach[END_REF]], these non-initial schwas have a higher frequency of elision than initial ones. Indeed, only 155 occurrences (21%) show acoustic vocalic information. The schwa duration of these schwa-full forms (43 ms, σ = 17) is shorter than that of initial schwas (t(2968) = 5.626, p<0.0001), the coefficients of skewness and kurtosis are respectively of 0.32 and -0.84.

A separate analysis was also done according to the position of the schwa relative to morphological boundaries. Schwas were coded as being either morpheme-internal (3218 occurrences, e.g. 'cheval', horse) or at a morpheme boundary (938 occurrences, e.g. 'retrouver', found again), according to a rather conservative morphological analysis. Three words ('devenir', 'remarque' and 'remarquer', corresponding to 29 occurrences) for which categorisation was uncertain were excluded. The distributions of the duration of both morpheme-internal and non-internal schwas are similar to the ones presented above: schwa-less forms on one side (31% of morpheme internal schwas; 24% of non-internal schwas) and schwa-full forms on the other side, the two subgroups being separated by a gap in the distribution between 0 and 16 ms.

Taken together, these results show that the typical bimodal distribution of the duration of the schwas in our corpus is not dependent on the position of the schwa in the word or of the morpheme boundary. All the types of schwas analyzed are distributed over two separate subgroups: schwa-full forms with visible vocalic information and schwa-less forms defined by a lack of voicing and formant structure.

Figure 3: Schwa-full occurrences of the words 'sera' (will be) and 'semaine' (week). Note that the devoiced portion of the vowel in 'semaine' is not included in our segmentation, according to our acoustic criteria.

Discussion and conclusion

Optional schwa elision in French is defined as an alternation between the schwa-less and the schwa-full form of a word. In our corpus of radio broadcasted speech, we indeed observed that 29% of the occurrences of the words are realized with no voicing and formant structure corresponding to a schwa vowel. This process of schwa elision is considered in traditional phonological accounts as a categorical process of segmental deletion. In this study of the production of 4185 optional schwas, we have shown that there is a gap in the duration distribution of these schwas. Thus, two groups emerge with forms realized without a vocalic portion on one side, and forms realized with a vocalic portion of some duration on the other side.

Recall that according to our hypotheses, a continuous distribution, without a gap between schwa-full and schwa-less forms would have supported the view of a gradual reduction process, in which the schwa-less forms is the end-point. Our results, showing a bimodal distribution separated by a gap are on the contrary in support of the categorical nature of the process.

However, our results do not allow us to be assertive on this interpretation. Indeed, a categorical account is challenged by the fact that the gap we observe between the two groups is rather small. No schwas occur with a duration lesser than 16 ms. It corresponds roughly to 2 periods for a fundamental frequency of 125Hz. Even though we have been very strict on our segmentation criteria, we cannot exclude a measurement error rate of 1 or 2 periods.

Along the same line, it could be argued that the gap we observe in the distribution is a by-product of the criteria we chose to define the acoustical interval corresponding to the vowel. While voicing and formant structure are the most salient cues defining a vowel (and the clearest attributes to consider when dealing with such a large corpus), we cannot ignore the fact that vowels can reduce greatly in running speech and can surface without voicing or a clear formant structure. As stated earlier devoiced vowels have not been considered here as schwa-full realizations. Nonetheless, we are reluctant to consider that modifying our segmentation criteria to include such realizations would have change the distribution to support a continuous gradual process. The spectrogram presented in figure 4 illustrates such a case. The change in intensity between the end of the [v] and the [n] could be attributed to a devoiced vowel. The portion of noise with increased energy is around 25ms. Hence, if this realization was categorized as schwa-full, it would not have inflated our distribution in the gap region. Similarly, if our shortest vowels, such as the one of 'semaine' in figure 3, had been segmented to include some potential devoiced or formant less portions, their duration would have been longer. This would have changed the shape of our distribution by increasing the gap between the 2 groups, but not by decreasing it. An alternative explanation to consider that gives a full account of our data, is that the distribution we observed results from both a categorical and a gradual process. On the one hand, a phonological categorical process would account for the leftmost mode of the distribution. This process corresponds to the traditional optional schwa elision phenomena typical of French, in which the lexical variant produced is schwa-less. Whether this form results from a segmental deletion at one representation level or whether it corresponds to a stored lexical variant is beyond the scope of this paper. On the other hand, a gradual phonetic process of reduction would account for the gradual temporal reduction we observe within the group of the schwa-full forms. As do other segments, schwa is subject to reduction in continuous speech. We are not addressing here the question of whether this reduction is due to a mechanism of overlap, gestural reduction or target undershoot.

This dual interpretation has to be confirmed by the thorough analysis of the 1215 occurrences categorized as being schwa-less in our study. Of great interest will be the distribution of schwa-less forms and the potentially reduced schwa-full forms of specific lexical items and particular speaker productions. Further analysis is also needed to verify whether the dispersion in the realization of the schwa vowels mirrors that of other French vowels in the corpus.

Figure 1 :

 1 Figure 1 : Distribution of the schwa durations

Figure 4 :

 4 Figure 4 :Realization of the word 'avenir' (future)

The authors wish to thank G. Gravier, R. Ridouane and U. Frauenfelder for their help and suggestions.


Amine Ammar
email: amine.ammar@ujf-grenoble.fr

Francisco Chinesta
email: francisco.chinesta@paris.ensam.fr

David Ryckelynck
email: david.ryckelynck@paris.ensam.fr

Deterministic particle approach of Multi Bead-Spring polymer models

Keywords: Fokker-Planck equation, Multi Bead Spring models, method of particles. MOTS-CLÉS : équation de Fokker-Planck, modèles Multi-Haltères, méthodes particulaires

. Dans un travail antérieur (Ammar and Chinesta, 2005), une approche déterministe a été proposée pour traiter le modèle cinétique des suspensions de fibres courtes où, après introduction du terme de diffusion dans celui d'advection, l'équation de convection pure résultante est discrétisée à l'aide d'une technique de particules. Dans ce travail, cette approche déterministe est appliquée à la simulation du modèle Multi-Haltère pour des écoulements homogènes afin d'évaluer sa pertinence ainsi que son impact sur le nombre de particules pour de grandes dimensions de l'espace des conformations.

Introduction

As a consequence of the increasing use of composite materials, micro-macro simulation comes back into fashion nowadays with the development of computer performance. However, it remains crucial that simplicity and performance of numerical techniques affect the usability of microscopic models.

In the case of dilute polymer, simple dumbbell models, such as the FENE model, are not representative of the behaviour of a polymeric chain. In fact, both, the flexibility of the chain and the description of its complex configurations state are better accounted using a Multi Bead-Spring (MBS) representation.

Simulation of MBS has been widely performed by using stochastic techniques, whose accuracy is related to the number of realizations considered. In the context of deterministic techniques, the use of standard discretization techniques such as finite differences or finite elements are forbidden, because of the extremely large number of degrees of freedom involved due to the multidimensional character of the associated Fokker-Planck equation. Thus, the direct resolution of the Fokker-Planck equation in the finite element framework only concerned simple dumbbell models [START_REF] Lozinski | A fast solver for Fokker-Planck equation applied to viscoelastic flows calculations: 2D FENE model[END_REF][START_REF] Chauviere | Simulation of dilute polymer solutions using a Fokker-Planck equation[END_REF]). An alternative to the standard discretization techniques of the Fokker-Planck equation consists of the use of the method of particles. However, this technique has never been used until now for solving models defined in a multidimensional space, as in the case of the MBS polymer models. The main aim of this work is the analysis of the application of this technique to those models.

Results of such simulation highlight the flexibility effects of the chain (especially in shearing flows), which cannot be observed in the case of simple dumbbell models. However, a number of difficulties persist, some of them restricting seriously its application in some simulation conditions.

Mechanical model

As shown in Figure 1, the MBS chain consists of N+1 beads connected by N springs. The bead serves as an interaction point with the solvent and the spring contains the local stiffness information depending on local stretching (see Bird et al., 1987, for more details).

Figure 1. MBS model for polymer chain

The dynamics of the chain is governed by hydrostatic, Brownian and connector forces. If we denote by j r the velocity of the bead 'j' and by j q the velocity of the spring connector, then we obtain:

N j r r q j j j ,.., 1 , 1 [1]
The dynamics of each bead can be written as:

0 1 Hydrostatic effects Interactions Forces Brownian effects ( grad ) ln( ) 0 c c j j b j j j r v v r k T F F r [2]
where is the drag coefficient, v is the velocity field, v 0 is an average velocity, k b is the Boltzman constant, T is the absolute temperature and is the probability distribution function. From Equations [1] and [2] it results:

1 1 grad ln( ) N c j j j k b k k k q v q A k T F q [3]
where A jk is the Rouse Matrix defined by:

1 if 1 if 2 k j k j A jk
In the Rouse model connector force F c is a linear function of the connector orientation vector. But in our case a Finitely Extensible Non Linear spring will be used so that each spring branch behaves as in FENE model whose connector force will be detailed later.

Using now the evolution of the distribution function given by N j N j j N t q q q q t t q q 1 1 1 ) , ,.., ( ) , ,.., ( [4] and introducing Equation [3] it results:

1 1 2 1 1 1 grad + N N c j j k k j k j N N b jk j k j k v q A F t q k T A q q [5]
Moreover, Equation [3] can be rewritten excluding Brownian effects which are retained in Equation [5]. Then, Equations [3] and [5] can be expressed, in a dimensionless form (see Bird et al., 1987, for more details), as:

1 grad ( ) 4 c q vq AF q [6] 1 ( ) 4 q A t q q q =7]
Where q denotes the vector containing the N spring connector vectors j q , with 1 j N , , .

The problem [6]- [7] is defined in a MxN-dimensional space, being M=1,2 or 3 (the dimension of the physical space) and N the number of springs that compose the polymeric chain. In the present work the connector force is given by:

j j j c q b q q F / 1 1 ) ( 2 [8]
where b is the maximum stretching of each spring connector of the chain.

Deterministic discretization approach using the method of particles

Advection form of the Fokker-Planck equation

For applying the method of particles to discretize the Fokker-Planck equation, we firstly introduce its diffusion term into the advection one (Ammar and Chinesta 2005;[START_REF] Chinesta | On the solution of the Fokker-Planck equations in steady recirculating flows involving short fiber suspensions[END_REF][START_REF] Chaubal | Smoothed particle hydrodynamics techniques for the solution of kinetic theory problems[END_REF]. Thus, it results the pseudoadvection field:

~1 4 q q q A [9]
being the associated Fokker-Planck equation:

( ) q t q [10] Problem [9]-[10] is equivalent to the initial problem [6]-[7].

Smoothed particle approximation

In order to discretize Equations [9]-[10] one must define the approximation of the unknown field which is, in our case, the distribution function. The method of particles uses a Dirac representation, consisting of ( )

Cn i i i
q q , where the superscript refers to each chain of the population involved in the simulation that consists of Cn chains (particles). However, from a computational point of view the derivative of Dirac masses cannot be defined in a proper way. So, usually, a smoothed approximation of the Dirac mass is introduced according to: where denotes a smoothing parameter and 'd' denotes a distance in the MxNdimensional space of molecule configurations.

Thus we can write: ( )

i d i i q q d [12]
where i i d q q is the Euclidean norm in the MxN-dimensional space.

The smoothed distribution can then be written as: 1 ( )

Cn i i i q d [13]
Where i is the weight related to each chain, which, in this work, is considered the same overall the population ( = 1/Cn for satisfying the normality condition of the distribution function).

Equilibrium conditions

Before starting simulations, an equilibrium state eq can be defined, and then used for defining the initial condition for start-up problems, as well as for adjusting the value of the smoothing parameter by comparing eq with the solution computed from the steady state related to the transient solution of the Fokker-Planck equation without the advection term.

According to Bird et al., (1987) and after the application of the normality condition in the 2xN dimensional case considered in this work, the expression of eq (for the connector force defined in Equation [8] is given by:

N N j b j eq b b b q q 2 2 1 ) ( 1 2 / 2
[14]

Simulation algorithm

The simulation algorithm proceeds as follows: -Consider Cn particles (molecules) consisting of N connectors whose orientation must describe the initial known distribution: (0) , 1, , ; 1, , i j q i Cn j N that is, (0) i q is assumed known.

Index i refers to the chain (particle) involved in the simulation, j to the connector in the chain, and (n) indicates the time step.

-For each time step: -Compute the advection field according to Equation [9]:

( ) ( ) ( ) ~( ) ( ) ( ) 1 4 ( ) i n n i n q i n i n n q q q A q where ( ) ( ) 1
( )

Cn n i n i i q d , being ( ) ( ) i n i n d q q
. Update the molecular conformation:

( ) ( 1) ( ) i n i n i n q q q t
Update the distribution function:

( 1) ( 1) 1 ( )

Cn n i n i i q d

Numerical test

In order to validate the approach just proposed we consider firstly a problem whose solution can be obtained analytically. This is for example the case of a single FENE dumbbell in a one-dimensional physical space. This problem was deeply analysed in [START_REF] Keunings | On the peterlin approximation for finitely extensible dumbbells[END_REF]. The number of particles involved in the simulation is in the order of 2400, being the other simulation parameters the following: 0.001 t [15]

that is equivalent to Equation [14] after normalization. We prefer using in this numerical example Equation [15] as initial condition instead of Equation [14] because in this manner the computed solution can be directly compared with the results reported in [START_REF] Keunings | On the peterlin approximation for finitely extensible dumbbells[END_REF].

The evolution of the distribution function is depicted in Figure 2 which represents the solutions at times: t=0,0.05,0.1,0.15,0.2,0.3,0.4,0.5,1.

Figure 2. Evolution of the distribution function related to a single 1D FENE model

Figure 3 zooms the region located around the maximum molecule stretching. These distributions are in perfect agreement with the exact solution as well as with the numerical solutions computed using the stochastic technique that can be found in [START_REF] Keunings | On the peterlin approximation for finitely extensible dumbbells[END_REF]. Figure 4 depicts the histogram related to the particle position at t=1 (we recall that the portion of the particle indicates the position of one of its ends, the other one remains at the origin x=0). Simular results were obtained by using only around 600 particles. Now, after proving the accuracy of the proposed strategy, we consider higher dimensional problems. For this purpose we consider a MBS -FENE model defined in a one-dimensional physical space composed of two dumbbells. In this case a reference solution can be easily obtained by using a more experienced numerical technique as is the finite element method. We consider the problem defined by the following parameters: Cn=500 (number of particles), b=10 ( b being the dimensionless parameter giving the maximum stretching), N = 2 (number of branches), M=1 (dimension of the physical space) and 0.35 du We dx

.

Figure 5 depicts the reference steady solution as well as the solution computed by using the method of particles with a smoothing parameter 0.4 . The steady distribution function depends on the stretching of both branches 1 q and 2 q .

Figure 5. Distribution function related to a 2 branches FENE MBS molecule model conformation: (left) Reference solution; (right) Solution computed using the method of particles

The distribution depicted in Figure 5(right) is related to the particle distribution depicted in Figure 6. In order to conclude about the incidence of the smoothing parameter on the computed solution we show in Figure 7 the distributions computed by using different values of that parameter, in particular 0.1,0.2,0.3, 0.5 ( 0.4 being the optimum value whose associated distribution is depicted in Figure 5(right)). We can notice the significant impact of the smoothing parameter on the particles distribution and then on the conformation distribution. Moreover this parameter must be identified in each particular application because it depends on the number of particles considered in the simulation; on the solution itself (more localized is the solution smaller must be the smoothing parameter, but at the same time it must be large enough for coupling particles in the regions with low particles density); and on the dimension of the physical and conformation spaces. The definition of an appropriate adaptive procedure for evaluating the local value of the smoothing parameter is not an easy task. 0.1 (a), 0.2 (b), 0.3 (c) and 0.5 (d) Finally, to analyse the capabilities of the approach just described, for treating more realistic models, a MBS -FENE model composed of many dumbbells in a 2D physical space, is considered. This test is carried out in a simple shear flow with the following simulation conditions: b = 9, 1 du dy , 0.01 t , N = 5, M=2, Cn = 52, = 0.8 and being the initial condition given by Equation [14]. In all the simulations the flow kinematics is assumed known and uncoupled with the polymer presence.

The solutions at times t1 = 3 and t2 = 10 (at t2 the steady state is reached) are depicted in Figure 8. To simplify the solution post-processing, a smoothed distribution involving the connectors belonging to the same level (having the same index in the chain) is defined according to:

1 1 1 , 1, , Cn Cn i i j i j i j i i
q q q q j N Cn [16] where i j q refers the connector j of chain i. From the computed results, we notice that symmetric levels (in our case levels 1 and 5, or 2 and 4) exhibit the same behaviour. Connectors are more stretched in the centre of the chain than on its ends. Figures 8(a) and 8(b) depicts 1 and 3 at time t1 and Figures 8(c) and 8(d) the same distributions at time t2. We can also notice the perfect symmetry of the computed distributions since the initial distributions respect symmetry. Figure 9 depicts the chain configurations at time t2 of some of the chains involved in the simulation. In this figure the origin of the first spring connector q 1 of all chains is located at the centre of the figure.

Conclusions

Despite of its simplicity, the deterministic particle technique just described seems to be robust and accurate when the chain model consists of a reduced number of connectors. However a number of difficulties persist, being the most important:

-the number of particles increases significantly with the dimension of the conformation space as was also found in the context stochastic simulations;

-the smoothing parameter has a strong impact on the computed results, even if it can be identified from the equilibrium conditions, when the distribution becomes highly localized an adaptation of this value is required to guaranty the solution accuracy. From some numerical experiments we have conclude that each particle must have a minimum number of neighbours that seems to be around 3;

-it is well known that in highly multidimensional spaces the use of a Euclidean distance results delicate, because in a high dimensional space the measure of the surface of a hyper-sphere becomes higher than the measure of its volume.

.

  At the initial time the particles are distributed for verifying the equilibrium initial condition given by:

Figure 3 .

 3 Figure 3. Zoom of the distribution functions in the region close to the maximum moldeule stretching

Figure 4 .

 4 Figure 4. Histogram related to the particle position at t=1

Figure 6 .Figure 7 .

 67 Figure 6. Particles distribution related to the solution depicted in figure 5(right)

Figure 8 .

 8 Figure 8. MBS conformation in a simple shear flow

Figure 9 .

 9 Figure 9. Example of chain representation at time t2

  Ammar A., Chinesta F., A Particle Strategy for Solving the Fokker-Planck Equation Governing the Fiber Orientation Distribution in Steady Recirculating Flows Involving Short Fiber Suspensions, Lectures Notes on Computational Science and Engineering, Springer, 43, 2005, p. 1-16. Bird R.B., Crutiss C.F., Armstrong R.C., Hassager O., Dynamic of polymeric liquid, vol. 2, Kinetic Theory, John Wiley & Sons, 1987.