Alessandra Bernardi

Jérôme Brachat

Pierre Comon
email: pcomon@i3s.unice.fr

B Mourrain

Bernard

GENERAL TENSOR DECOMPOSITION, MOMENT MATRICES AND APPLICATIONS

Keywords: tensor, decomposition, multihomogeneous polynomial, rank, Hankel operator, moment matrix, flat extension

The tensor decomposition addressed in this paper may be seen as a generalisation of Singular Value Decomposition of matrices. We consider general multilinear and multihomogeneous tensors. We show how to reduce the problem to a truncated moment matrix problem and give a new criterion for flat extension of Quasi-Hankel matrices. We connect this criterion to the commutation characterisation of border bases. A new algorithm is described. It applies for general multihomogeneous tensors, extending the approach of J.J. Sylvester to binary forms. An example illustrates the algebraic operations involved in this approach and how the decomposition can be recovered from eigenvector computation.

Introduction

Tensors are objects that appear in various contexts and applications. Matrices are tensors of order two, and are better known than tensors. But in many problems, higher order tensors are naturally used to collect information which depend on more than two variables. Typically, these data could be observations of some experimentation or of a physical phenomenon that depends on several parameters. These observations are stored in a structure called tensor, whose dimensional parameters (or modes) depend on the problem.

The tensor decomposition problem consists of decomposing a tensor (e.g. the set of observations) into a minimal sum of so-called decomposable tensors (i.e. tensors of rank 1). Such a decomposition which is independent of the coordinate system allows to extract geometric or invariant properties associated with the observations. For this reason, the tensor decomposition problem has a large impact in many applications. The first well known case is encountered for matrices (i.e. tensors of order 2), and is related to Singular Value Decomposition with applications e.g. to Principal Component Analysis. Its extension to higher order tensors appears in Electrical Engineering [START_REF] Swami | Multichannel ARMA processes[END_REF], in Signal processing [START_REF] Comon | Blind identification of under-determined mixtures based on the characteristic function[END_REF], [START_REF] Cichocki | Adaptive Blind Signal and Image Processing[END_REF], in Antenna Array Processing [START_REF] Dogan | Applications of cumulants to array processing .I. aperture extension and array calibration[END_REF] [START_REF] Chevalier | On the virtual array concept for higher order array processing[END_REF] or Telecommunications [START_REF] Van Der Veen | An analytical constant modulus algorithm[END_REF], [START_REF] Chevalier | Optimal separation of independent narrow-band sources -concept and performance[END_REF], [START_REF] Sidiropoulos | Blind PARAFAC receivers for DS-CDMA systems[END_REF], [START_REF] Ferreol | On the behavior of current second and higher order blind source separation methods for cyclostationary sources[END_REF], [START_REF] De Lathauwer | Tensor-based techniques for the blind separation of ds-cdma signals[END_REF], in Chemometrics [START_REF] Bro | Parafac, tutorial and applications[END_REF] or Psychometrics [START_REF] Kiers | An efficient algorithm for Parafac of three-way data with large numbers of observation units[END_REF], in Data Analysis [START_REF] Comon | Independent Component Analysis[END_REF], [START_REF] Cardoso | Blind signal separation: statistical principles[END_REF], [START_REF] Donoho | Uncertainty principles and ideal atomic decompositions[END_REF], [START_REF] Jiang | Kruskal's permutation lemma and the identification of CANDE-COMP/PARAFAC and bilinear models[END_REF], [START_REF] Smilde | Multi-Way Analysis[END_REF], but also in more theoretical domains such as Arithmetic complexity [START_REF]Three-way arrays: Rank and uniqueness of trilinear decompositions[END_REF] [8] [START_REF] Strassen | Rank and optimal computation of generic tensors[END_REF] [START_REF] Landsberg | Geometry and the complexity of matrix multiplication[END_REF]. Further numerous applications of tensor decompositions may be found in [START_REF] Cichocki | Adaptive Blind Signal and Image Processing[END_REF] [START_REF] Comon | Handbook of Blind Source Separation, Independent Component Analysis and Applications[END_REF] [START_REF] Smilde | Multi-Way Analysis[END_REF].

From a mathematical point of view, the tensors that we will consider are elements of T := S δ1 (E 1 ) ⊗ • • • ⊗ S δ k (E k ) where δ i ∈ N, E i are vector spaces of dimension n i + 1 over a field K (which is of characteristic 0 and algebraically closed), and S δi (E i ) is the δ th i symmetric power of E i . The set of tensors of rank 1 form a projective variety which is called the Veronese variety when k = 1 or the Segre variety when δ i = 1, i = 1, . . . , k. We will call it hereafter the Segre-Veronese variety of P(T ) and denote it Ξ(T ). The set of tensors which are the linear combinations of r elements of the Segre-Veronese variety are those which admit a decomposition with at most r terms of rank 1 (i.e. in Ξ(T )). The closure of this set is called the r-secant variety and denoted Ξ r (T ). More precise definitions of these varieties will be given in Sec. 2.3.

Decomposing a tensor T consists of finding the minimal r such that this tensor is a sum of r tensors of rank 1. This minimal r is called the rank of T . By definition, a tensor of rank r is in the secant variety Ξ r (T ). Thus analysing the properties of these secant varieties and their characterisation helps determining tensor ranks and decompositions.

The case where k = 2 and δ 1 = δ 2 = 1 corresponds to the matrix case, which is well known. The rank of a matrix seen as a tensor of order k = 2 is its usual rank. The case where k = 1 and δ 1 = 2 corresponds to the case of quadratic forms and is also well understood. The rank of the symmetric tensor is the usual rank of the associated symmetric matrix. The case where k = 1, δ 1 ∈ N and n 1 = 1 corresponds to binary forms, which has been analyzed by J.J. Sylvester in [START_REF] Sylvester | Sur une extension d'un théorème de Clebsch relatif aux courbes du quatrième degré[END_REF]. A more complete description in terms of secant varieties is given in [START_REF] Landsberg | Equations for secant varieties via vector bundles[END_REF].

On our knowledge if k > 1 and if at least one of the δ i 's is larger than 1, then there is no specific result in the literature on the defining ideal of secant varieties of Segre-Veronese varieties Ξ(S δ1 (E 1 ) ⊗ • • • ⊗ S δ k (E k )) except for [START_REF] Catalisano | On the ideals of secant varieties to certain rational varieties[END_REF] where the authors conjecture that when Ξ r (S δ1 (E 1 ) ⊗ S δ2 (E 2 ))) is a defective hypersurface, then its defining equation is a determinantal equation.

In the case of the secant varieties of Veronese varieties (i.e. if k = 1 and δ 1 > 1), the knowledge of their ideal is sparse. Beside the classical results (see one for all [START_REF] Iarrobino | Power sums, Gorenstein algebras, and determinantal loci[END_REF]) we quote [START_REF] Landsberg | Equations for secant varieties to veronese varieties[END_REF] as the most up-to-date paper on that subject. We also quote [START_REF] Landsberg | Equations for secant varieties via vector bundles[END_REF] for a modern approach to equations of secant varieties in general using representation theory.

About the case of secant varieties of Segre varieties (i.e. δ i = 1 for i = 1, . . . , k) the only obvious case is the 2 factors Segre. For some of the non trivial cases in which equations of secant varieties of Segre varieties are known we refer to [START_REF] Landsberg | On the ideals of secant varieties of Segre varieties[END_REF], [START_REF] Landsberg | On the ideals and singularities of secant varieties of Segre varieties[END_REF], [START_REF] Allman | Phylogenetic ideals and varieties for the general Markov model[END_REF], [START_REF] Catalisano | On the ideals of secant varieties to certain rational varieties[END_REF].

The first method to compute such a decomposition, besides the case of matrices or quadratic forms which may go back to the Babylonians, is due to Sylvester for binary forms [START_REF] Sylvester | Sur une extension d'un théorème de Clebsch relatif aux courbes du quatrième degré[END_REF]. Using apolarity, kernels of catalecticant matrices are computed degree by degree until a polynomial with simple roots is found. See also [START_REF] Comas | On the rank of a binary form[END_REF], [START_REF] Iarrobino | Power sums, Gorenstein algebras, and determinantal loci[END_REF]. An extension of this approach for symmetric tensors has been analyzed in [START_REF] Iarrobino | Power sums, Gorenstein algebras, and determinantal loci[END_REF], and yields a decomposition method in some cases (see [START_REF] Iarrobino | Power sums, Gorenstein algebras, and determinantal loci[END_REF][p. 187]). Some decomposition methods are also available for specific degrees and dimensions, e.g. using invariant theory [START_REF] Comon | Decomposition of quantics in sums of powers of linear forms[END_REF]. In [START_REF] Bernardi | Computing symmetric rank for symmetric tensors[END_REF], there is a simplified version of Sylvester's algorithm, which uses the mathematical interpretation of the problem in terms of secant varieties of rational normal curves. The same approach is used in [START_REF] Bernardi | Computing symmetric rank for symmetric tensors[END_REF] to give algorithms for the decompositions of symmetric tensors belonging to Ξ 2 (S d (E)) and to Ξ 3 (S d (E)). In [START_REF] Ballico | Stratification of the fourth secant variety of veronese variety via symmetric rank[END_REF] a complete rank stratification of Ξ 4 (S d (E)) is given.

In [START_REF] Brachat | Symmetric tensor decomposition[END_REF], Sylvester's approach is revisited from an affine point of view and a general decomposition method based on a flat extension criterion is described. The main contribution of the current paper is to extend this method to more general tensor spaces including classical multilinear tensors and multihomogeneous tensors. In particular we give a new and more flexible criterion for the existence of a decomposition of a given rank, which extends non trivially the result in [START_REF] Laurent | A Sparse Flat Extension Theorem for Moment Matrices[END_REF] and the characterization used in [START_REF] Brachat | Symmetric tensor decomposition[END_REF]. This criterion is a rank condition of an associated Hankel operator. Moreover we use that criterion to write a new algorithm which checks, degree by degree, if the roots deduced from the kernel of the Hankel operator are simple. This allows to compute the rank of any given partially symmetric tensor.

This paper is an extended version of [START_REF] Bernardi | Multihomogeneous polynomial decomposition using moment matrices[END_REF], with the complete proofs and with detailed examples.

In Sec. 2, we recall the notations, the geometric point related to secants of Segre and Veronese varieties, and the algebraic point of view based on moment matrices. In Sec. 3, we describe the algorithm and the criterion used to solve the truncated moment problem. In Sec. 4, an example of tensor decompositions from Antenna Array Processing illustrates the approach.

2. Duality, moment matrices and tensor decomposition 2.1. Notation and preliminaries. Let K be an algebraically closed field (e.g. K = C the field of complex numbers). We assume that K is of characteristic 0. For a vector space E, its associated projective space is denoted P(E). For v ∈ E -{0} its class in P(E) is denoted v. Let P n be the projective space of E := K n+1 . For a subset F = {f 1 , . . . , f m } of a vector-space (resp. ring) R, we denote by F (resp. (F )) the vector space (resp. ideal) generated by F in R.

We consider hereafter the symmetric δ-th power S δ (E) where E is a vector space of basis x 0 , . . . , x n . An element of S δ (E) is a homogeneous polynomial of degree δ ∈ N in the variables x = (x 0 , . . . , x n ). For x 1 = (x 0,1 , . . . , x n1,1 ) , . . . ,

x k = (x 0,k , . . . , x n k ,k ), S δ1 (E 1 ) ⊗ • • • ⊗ S δ k (E k ) (with E i = x 0,i , . . . , x ni,i
) is the vector space of polynomials multihomogeneous of degree δ i in the variables x i .

Hereafter, we will consider the deshomogeneisation of elements in

S δ1 (E 1 )⊗• • •⊗S δ k (E k
), obtained by setting x 0,i = 1 for i = 1, . . . , k. We denote by R δ1,...,δ k this space, where R = K[x 1 , . . . , x k ] is the space of polynomials in the variables x 1 = (x 1,1 , . . . , x n1,1 ), . . . ,

x k = (x 1,k , . . . , x n k ,k ).
For

α i = (α 1,i , . . . , α ni,i ) ∈ N ni (i = 1, . . . , k), let x αi i = ni j=1 x αj,i j,i , |α i | = ni j=1 α j,i
, and

x α = ni j=1 x αi i . An element f of R δ = R δ1,...,δ k is represented as f = α=(α1,...,α k );|αi|≤δi f α x α1 1 • • • x α k k . The dimension of R δ := R δ1,...,δ k is n δ1,...,δ k ;n1,...,n k = k i=1 ni+δi δi . For δ ∈ N, α ∈ N n with |α| ≤ δ, let δ α = δ! α1!•••αn!(δ-|α|)! .
We define the apolar inner product on R δ1,...,δ k by

f |g = |αi|≤δi f α g α δ1 α1 • • • δ k α k .
The dual space of a K-vector space E is denoted E * = Hom K (E, K). It is the set of Klinear forms from E to K. A basis of the dual space R * δ , is given by the set of linear forms that compute the coefficients of a polynomial in the monomial basis (x α ) α∈N n 1 ו••×N n k ;|αi|≤δi . We denote it by (d α ) α∈N n 1 ו••×N n k ;|αi|≤δi . We identify R * with the (vector) space of formal power series K

[[d]] = K[[d 1 , . . . , d k ]] = K[[d 1,1 , . . . , d n1,1 , . . ., d 1,k , . . . , d n k ,k ]]. Any element Λ ∈ R * can be decomposed as Λ = α∈N n 1 ו••×N n k Λ(x α ) d α .
Typical elements of R * are the linear forms that correspond to the evaluation at a point

ζ = (ζ 1 , . . . , ζ i ) ∈ K n1 × • • • × K n k : 1 ζ : R → K p → p(ζ)
The decomposition of

1 ζ in the basis {d α } α∈N n 1 ו••×N n k is 1 ζ = α∈N n 1 ו••×N n k ζ α d α = α∈N n 1 ו••×N n k k i=1 ζ αi i d αi i .
We recall that the dual space R * has a natural structure of R-module [START_REF] Elkadi | Introduction à la résolution des systèmes polynomiaux[END_REF] which is defined as follows: for all p ∈ R, and for all Λ ∈ R * consider the linear operator p ⋆ Λ : R → K q → Λ(pq).

In particular, we have

x i,j ⋆ d α1 1 • • • d αj j • • • d α k k = d α 1 1 • • • d α j-1 j-1 d α 1,j 1,j • • • d α i-1,j i-1,j d α i,j -1 i,j d α i+1,j i+1,j • • • d α n j ,j n j ,j d α j+1 j+1 • • • d α k k if α i,j
> 0 and 0 otherwise.

Tensor decomposition.

In this section, we present different formulations of the tensor decomposition problem, that we consider in this paper. We will consider hereafter a partially symmetric tensor

T of S δ1 (E 1 )⊗• • •⊗S δ k (E k
) where E i = x 0,i , . . . , x ni,i . It can be represented by a partially symmetric array of coefficients [START_REF] Abo | Induction for secant varieties of Segre varieties[END_REF] [T ] = (T α1,...,α k ) αi∈N n i +1 ;|αi|=δi .

For α i ∈ N ni with |α i | ≤ δ i , we denote α i = (δ i -|α i |, α 1,i , . . ., α ni,i
) and, with an abuse of notation, we identify T α1,...,α k := T α1,...,α k . Such a tensor is naturally associated with a (multihomogeneous) polynomial in the variables x 1 = (x 0,1 , . . . , x n1,1 ), . . ., x k = (x 0,k , . . . , x n k ,k )

T (x) = α=(α 1 ,...,α k )∈N n 1 ו••×N n k ; |α i |≤δ i T α x α1 1 • • • x α k k .
or to an element T (x) ∈ R δ1,...,δ k obtained by substituting x 0,i by 1 in T (x) (for i = 1, . . . , k):

T (x) = α∈N n 1 ו••×N n k ; |α i |≤δ i T α x α1 1 • • • x α k k . An element of R * = K[[d]]
can also be associated naturally with T :

T * (d) = α∈N n 1 ו••×N n k ; |α i |≤δ i δ 1 α 1 -1 • • • δ k α k -1 T α d α1 1 • • • d α k k .
so that for all T ′ ∈ R δ1,...,δ k ,

T (x)|T ′ (x) = T * (d)(T ′ (x)).
The decomposition of tensor T can be stated as follows:

Tensor decomposition problem. Given T (x) ∈ S δ1 (E 1 ) ⊗ • • • ⊗ S δ k (E k )
, find a decomposition of T (x) as a sum of products of powers of linear forms in x j :

(2)

T (x) = r i=1 γ i l 1,i (x 1 ) δ1 • • • l k,i (x k ) δ k
where γ i = 0, l j,i (x j ) = l 0,j,i x 0,j + l 1,j,i x 1,j + • • • + l nj ,j,i x j,nj and r is the smallest possible integer for such a decomposition.

Definition 2.1. The minimal number of terms r in a decomposition of the form (2) is called the rank of T .

We say that T (x) has an affine decomposition if there exists a minimal decomposition of T (x) of the form (2) where r is the rank of T and such that l 0,j,i = 0 for i = 1, . . . , r. Notice that by a generic change of coordinates in E i , we may assume that all l 0,j,i = 0 and thus that T has an affine decomposition. Suppose that T (x) has an affine decomposition. Then by scaling l j,i (x j ) and multiplying γ i by the inverse of the δ th j power of this scaling factor, we may assume that l 0,j,i = 1. Thus, the polynomial

T (x) = r i=1 γ i |αi|≤δi δ 1 α 1 • • • δ k α k ζ α1 1,i • • • ζ α k k,i x α1 1 • • • x α k k with T α1,...,α k = r i=1 γ i |αi|≤δi δ1 α1 • • • δ k α k ζ α1 1,i • • • ζ α k k,i
. Equivalently, we have

T * (d) = r i=1 γ i |αi|≤δi ζ α1 1,i • • • ζ α k k,i d α1 1 • • • d α k k so that T * (d) coincides on R δ1,...,δ k with the linear form r i=1 γ i 1 ζ1,i,...,ζ k,i = r i=1 γ i 1 ζi with ζ i := (ζ 1,i , . . . , ζ k,i ) ∈ K n1 × • • • K n k .
The decomposition problem can then be restated as follows: Interpolation problem. Given T * ∈ R * δ1,...,δ k which admits an affine decomposition, find the minimal number of non-zero vectors

ζ 1 , . . . , ζ r ∈ K n1 × • • • × K n k and non-zero scalars γ 1 , . . . , γ r ∈ K -{0} such that (3) T * = r i=1 γ i 1 ζi on R δ1,...,δ k .
If such a decomposition exists, we say that Λ =

r i=1 γ i 1 ζi ∈ R * extends T * ∈ R * δ1,...,δ k .
2.3. Decomposable tensors. In this section, we analyze the set of tensors of rank 1, also called decomposable tensors [START_REF] Abo | Induction for secant varieties of Segre varieties[END_REF]. They naturally form projective varieties, which we are going to describe using the language of projective geometry. We begin by defining two auxiliary but very classical varieties, namely Segre and Veronese varieties.

Definition 2.2. The image of the following map

s k : P(E 1 ) × • • • × P(E k ) → P(E 1 ⊗ • • • ⊗ E k ) (v 1 , . . . , v k ) → v 1 ⊗ • • • ⊗ v k
is the so called Segre variety of k factors. We denote it by

Ξ(E 1 ⊗ • • • ⊗ E k ).
From Definition 2.1 of the rank of a tensor and from the Interpolation Problem point of view [START_REF] Atkinson | On the maximal multiplicative complexity of a family of bilinear forms[END_REF] we see that a Segre variety parametrizes projective classes of rank 1 tensors

T = v 1 ⊗ • • • ⊗ v k ∈ E 1 ⊗ • • • ⊗ E k for certain v i ∈ E i , i = 1, . . . , k.
Definition 2.3. Let (J 1 , J 2 ) be a partition of the set {1, . . . , k}. If J 1 = {h 1 , . . . , h s } and

J 2 = {1, . . . , k}\J 1 = {h ′ 1 , . . . , h ′ k-s }, the (J 1 , J 2 )-Flattening of E 1 ⊗• • •⊗E k is the following: E J1 ⊗ E J2 = (E h1 ⊗ • • • ⊗ E hs ) ⊗ (E h ′ 1 ⊗ • • • ⊗ E h ′ k-s ). Let E J1 ⊗ E J2 be any flattening of E 1 ⊗ • • • ⊗ E k as in Definition 2.3 and let f J1,J2 : P(E 1 ⊗• • •⊗E k ) → P(E J1 ⊗E J2 ) be the obvious isomorphism. Let [T ] be an array associated with a tensor T ∈ E 1 ⊗ • • • ⊗ E k ; let T ′ = f J1,J2 (T) ∈ P(E J1 ⊗ E J2 ) and let [A J1,J2 ] be the matrix associated with T ′ . Then the d-minors of the matrix [A J1,J2 ] are said to be d-minors of [T ]. An array [A] = (x i1,...,i k ) 0≤ij ≤nj , j=1,...,k is said to be a generic array of indeterminates of R = K[x 1 , . . ., x k ] if the entries of [A] are the independent variables of R.
It is a classical result due to R. Grone (see [START_REF] Grone | Decomposable tensors as a quadratic variety[END_REF]) that a set of equations for a Segre variety is given by all the 2-minors of a generic array. In [START_REF] Hà | Box-shaped matrices and the defining ideal of certain blowup surfaces[END_REF]

it is proved that, if [A] is a generic array in R of size (n 1 + 1) × • • • × (n k + 1) and I d ([A]) is the ideal generated by the d-minors of [A] , then I 2 ([A]
) is a prime ideal, therefore:

I(Ξ(E 1 ⊗ • • • ⊗ E k )) = I 2 ([A]).
We introduce now the Veronese variety. Classically it is defined to be the d-tuple embedding of P n into P ( n+d d )-1 via the linear system associated with the sheaf O(d) with d > 0. We give here an equivalent definition.

Let E be an n+ 1 dimensional vector space. With the notation S d (E) we mean the vector subspace of E ⊗d of symmetric tensors.

Definition 2.4. The image of the following map

ν d : P(E) → P(S d (E)) v → v ⊗d
is the so called Veronese variety. We indicate it with Ξ(S d (E)).

With this definition it is easy to see that the Veronese variety parametrizes symmetric rank 1 tensors.

Observe that if we take the vector space E to be a vector space of linear forms x 0 , . . . , x n then the image of the map ν d above parametrizes homogeneous polynomials that can be written as d-th powers of linear forms.

The Veronese variety Ξ(S d (E)) ⊂ P(S d (E)) can be also viewed as

Ξ(S d (E)) = Ξ(E ⊗d ) ∩ P(S d (E)). Let [A] = (x i1,...,i d ) 0≤ij ≤n, j=1,.
..,d be a generic symmetric array. It is a known result that:

(4)

I(Ξ(S d (E))) = I 2 ([A]).
See [START_REF] Wakeford | On canonical forms[END_REF] for the set theoretical point of view. In [START_REF] Pucci | The Veronese variety and catalecticant matrices[END_REF] the author proved that I(Ξ(S d (E))) is generated by the 2-minors of a particular catalecticant matrix (for a definition of "Catalecticant matrices" see e.g. either [START_REF] Pucci | The Veronese variety and catalecticant matrices[END_REF] or [START_REF] Geramita | Catalecticant varieties[END_REF]). A. Parolin, in his PhD thesis ( [START_REF] Parolin | Varietà secanti alle varietà di segre e di veronese e loro applicazioni, tesi di dottorato[END_REF]), proved that the ideal generated by the 2-minors of that catalecticant matrix is actually

I 2 ([A]).
We are now ready to describe the geometric object that parametrizes partially symmetric tensors

T ∈ S δ1 (E 1 ) ⊗ • • • ⊗ S δ k (E k ).
Let us start with the rank 1 partially symmetric tensors.

Definition 2.5. Let E 1 , . . . , E k be vector spaces of dimensions n 1 +1, . . . , n k +1 respectively. The Segre-Veronese variety Ξ(S δ1 (E 1 )⊗• • •⊗S δ k (E k )) is the embedding of P(E 1 )⊗• • •⊗P(E k ) into P N -1 ≃ P(S δ1 (E 1 ) ⊗ • • • ⊗ S δ k (E k )), where N = Π k i=1 ni+δi di
, given by sections of the sheaf O(δ 1 , . . . , δ k ).

I.e. Ξ(S δ1 (E 1 ) ⊗ • • • ⊗ S δ k (E k ))
is the image of the composition of the following two maps:

P(E 1 ) × • • • × P(E k ) ν δ 1 ו••×ν δ k -→ P ( n 1 +δ 1 δ 1 )-1 × • • • × P ( n k +δ k δ k
)-1

and

P ( n 1 +δ 1 δ 1 )-1 × • • • × P ( n k +δ k δ t )-1 s -→ P N -1
, where each ν δi is a Veronese embedding of

P(E i ) as in Definition 2.4, then Im(ν δ1 × • • • × ν δ k ) = Ξ(S δ1 (E 1 )) × • • • × Ξ(S δ k (E k )) and
Im(s) is the Segre variety of k factors. Therefore the Segre-Veronese variety is the Segre re-embedding of the product of k Veronese varieties.

If (δ 1 , . . . , δ k ) = (1, . . . , 1) then the corresponding Segre-Veronese variety is nothing else than the classical Segre variety of

P(E 1 ⊗ • • • ⊗ E k ).
If k = 1 then the corresponding Segre-Veronese variety is nothing else than the classical Veronese variety of P(S δ1 (E 1 )).

Observe that Ξ(S δ1 (E

1 ) ⊗ • • • ⊗ S δ k (E k )) can be viewed as the intersection with the Segre variety Ξ(E ⊗δ1 1 ⊗ • • •⊗ E ⊗δ k k
) that parametrizes rank one tensors and the projective subspace

P(S δ1 (E 1 ) ⊗ • • • ⊗ S δ k (E k )) ⊂ P(E ⊗δ1 1 ⊗ • • • ⊗ E ⊗δ k k ) that parametrizes partially symmetric tensors: Ξ(S δ1 (E 1 ) ⊗ • • • ⊗ S δ k (E k )) = Ξ(E ⊗δ1 1 ⊗ • • • ⊗ E ⊗δ k k ) ∩ P(S δ1 (E 1 ) ⊗ • • • ⊗ S δ k (E k )).
In [START_REF] Bernardi | Ideals of varieties parameterized by certain symmetric tensors[END_REF] it is proved that if [A] is a generic array of indeterminates associated with the multihomogeneous polynomial ring

S δ1 (E 1 ) ⊗ • • • ⊗ S δ k (E k ) (i.e. it is a generic partially symmetric array), the ideal of the Segre-Veronese variety Ξ(S δ1 (E 1 ) ⊗ • • • ⊗ S δ k (E k )) is I(Ξ(S δ1 (E 1 ) ⊗ • • • ⊗ S δ k (E k ))) = I 2 ([A]) with δ i > 0 for i = 1, . . . , k.
Now if we consider the vector spaces of linear forms

E i ≃ S 1 (E i ) for i = 1, . . . , k, we get that the Segre-Veronese variety Ξ(S δ1 (E 1 ) ⊗ • • • ⊗ S δ k (E k )) parametrizes multihomogenoeus polynomials F ∈ S δ1 (E 1 ) ⊗ • • • ⊗ S δ k (E k ) of the type F = l δ1 1 • • • l δ k k where l i are linear forms in S 1 (E i ) for i = 1, . . . , k.
From this observation we understand that the tensor decomposition problem of finding a minimal decomposition of type (2) for an element T ∈ S δ1 (E 1 ) ⊗ • • • ⊗ S δ k (E k ) is equivalent to finding the minimum number of elements belonging to the Segre-Veronese variety

Ξ(S δ1 (E 1 ) ⊗ • • • ⊗ S δ k (E k )) whose span contains T ∈ P(S δ1 (E 1 ) ⊗ • • • ⊗ S δ k (E k )).
The natural geometric objects that are associated with this kind of problems are the higher secant varieties of the Segre-Veronese varieties that we are going to define. Definition 2.6. Let X ⊂ P N be any projective variety and define X 0 s := P1,...,Ps∈X P 1 , . . . , P s .

The s-th secant variety X s ⊂ P N of X is the Zariski closure of X 0 s . Observe that the generic element of X s is a point P ∈ P N that can be written as a linear combination of s points of X, in fact a generic element of X s is an element of X 0 s . Therefore if X is the Segre-Veronese variety, then the generic element of Ξ s (S δ1 (E

1 ) ⊗ • • • ⊗ S δ k (E k )) is the projective class of a partially symmetric tensor T ∈ S δ1 (E 1 ) ⊗ • • • ⊗ S δ k (E k )
that can be written as a linear combination of s linearly independent partially symmetric tensors of rank 1. Unfortunately not all the elements of Ξ s (S δ1 (E

1 ) ⊗ • • • ⊗ S δ k (E k )) are of this form. In fact if T ∈ Ξ s (S δ1 (E 1 ) ⊗ • • • ⊗ S δ k (E k )) \ Ξ 0 s (S δ1 (E 1 ) ⊗ • • • ⊗ S δ k (E k )) then the rank of T is strictly bigger than s. Definition 2.7. The minimum integer s such that T ∈ P(S δ1 (E 1 ) ⊗ • • • ⊗ S δ k (E k )) belongs to Ξ s (S δ1 (E 1 ) ⊗ • • • ⊗ S δ k (E k )) is called the border rank of T .
In order to find the border rank of a tensor

T ∈ S δ1 (E 1 ) ⊗ • • • ⊗ S δ k (E k ) we should need a set of equations for Ξ s (S δ1 (E 1 ) ⊗ • • • ⊗ S δ k (E k )) for s > 1.
The knowledge of the generators of the ideals of secant varieties of homogeneous varieties is a very deep problem that is solved only in very particular cases (see eg. [START_REF] Ottaviani | An invariant regarding Waring's problem for cubic polynomials[END_REF], [START_REF] Landsberg | On the ideals and singularities of secant varieties of Segre varieties[END_REF], [START_REF] Landsberg | Generalizations of Strassen's equations for secant varieties of Segre varieties[END_REF], [START_REF] Landsberg | Equations for secant varieties to veronese varieties[END_REF], [START_REF] Buczynski | Determinental equations for secant varieties and the Eisenbud-Koh-Stillman conjecture[END_REF], [START_REF] Landsberg | Equations for secant varieties via vector bundles[END_REF]).

From a computational point of view, there is a very direct and well known way of getting the equations for the secant variety, which consists of introducing parameters or unknowns for the coefficients of l i,j and γ i in (2), to expand the polynomial and identify its coefficients with the coefficients of T . Eliminating the coefficients of l i,j and γ i yields the equations of the secant variety.

Unfortunately this procedure is far from being computationally practical, because we have to deal with high degree polynomials in many variables, with a lot of symmetries. This is why we need to introduce moment matrices and to use a different kind of elimination.

2.4. Moment matrices. In this section, we recall the algebraic tools and the properties we need to describe and analyze our algorithm. Refer e.g. to [START_REF] Brachat | Symmetric tensor decomposition[END_REF], [START_REF] Elkadi | Introduction à la résolution des systèmes polynomiaux[END_REF], [START_REF] Mourrain | Multivariate Polynomials, Duality, and Structured Matrices[END_REF] for more details.

For any Λ ∈ R * , define the bilinear form

Q Λ , such that ∀a, b ∈ R, Q Λ (a, b) = Λ(ab). The matrix of Q Λ in the monomial basis, of R is Q Λ = (Λ(x α+β )) α,β , where α, β ∈ N n .
Similarly, for any Λ ∈ R * , we define the Hankel operator H Λ from R to R * as

H Λ : R → R * p → p ⋆ Λ.
The matrix of the linear operator H Λ in the monomial basis, and in the dual basis, {d α }, is

H Λ = (Λ(x α+β )) α,β
, where α, β ∈ N n . The following relates the Hankel operators with the bilinear forms. For all a, b ∈ R, thanks to the R-module structure, it holds

Q Λ (a, b) = Λ(ab) = a ⋆ Λ(b) = H Λ (a)(b).
In what follows, we will identify H Λ and Q Λ .

Definition 2.8.

Given B = {b 1 , . . . , b r }, B ′ = {b ′ 1 , . . ., b ′ r ′ } ⊂ R, we define H B,B ′ Λ : B → B ′ * ,
as the restriction of H Λ to the vector space B and inclusion of From the definition of the Hankel operators, we can deduce that a polynomial p ∈ R belongs to the kernel of H Λ if and only if p ⋆ Λ = 0, which in turn holds if and only if for all q ∈ R, Λ(pq) = 0. Proposition 2.9. Let I Λ be the kernel of H Λ . Then, I Λ is an ideal of R.

R * in B ′ * . Let H B,B ′ Λ = (Λ(b i b ′ j )) 1≤i≤r,1≤j≤r ′ . If B ′ = B,
Proof. Let p 1 , p 2 ∈ I Λ . Then for all q ∈ R, Λ((p 1 + p 2 )q) = Λ(p 1 q) + Λ(p 2 q) = 0. Thus, p 1 + p 2 ∈ I Λ .
If p ∈ I Λ and p ′ ∈ R, then for all q ∈ R, it holds Λ(pp ′ q) = 0. Thus pp ′ ∈ I Λ and I Λ is an ideal.

Let A Λ = R/I Λ be the quotient algebra of polynomials modulo the ideal I Λ , which, as Proposition 2.9 states is the kernel of H Λ . The rank of H Λ is the dimension of A Λ as a K-vector space.

Definition 2.10. For any B ⊂ R, let Λ is invertible, this implies that p = 0 and p = 0. As a consequence, we deduce that b 1 ⋆ Λ, . . . , b r ⋆ Λ are linearly independent elements of R * . This is so, because otherwise there exists m = [µ 1 , . . . , µ r ] ⊤ = 0, such that µ 1 (b

B + = B ∪ x 1 B ∪ • • • ∪ x n B and ∂B = B + \ B. Proposition 2.11. Assume that rank(H Λ ) = r < ∞ and let B = {b 1 , . . . , b r } ⊂ R such that H B Λ is invertible. Then b 1 , . . . , b r is a basis of A Λ . If 1 ∈ B the ideal I Λ is
1 ⋆ Λ) + • • • + µ r (b r ⋆ Λ) = (µ 1 b 1 + • • • + µ r b r ) ⋆ Λ = 0. As {b 1 , . . . , b r } ∩ Ker(H Λ ) = {0},
p i ∈ K[∂ 1 , . . . , ∂ n ], such that (5) Λ = d i=1 1 ζi • p i (∂)
Moreover the multiplicity of ζ i is the dimension of the vector space spanned the inverse system generated by 1 ζi • p i (∂).

Proof. Since rank(H Λ ) = r, the dimension of the vector space A Λ is also r. Thus the number of zeros of the ideal I Λ , say {ζ 1 , . . . , ζ d } is at most r, viz. d ≤ r. We can apply the structure Theorem [START_REF] Elkadi | Introduction à la résolution des systèmes polynomiaux[END_REF]Th. 7.34,p. 185] in order to get the decomposition.

In characteristic 0, the inverse system of 1 ζi • p i (∂) by p i is isomorphic to the vector space generated by p i and its derivatives of any order with respect to the variables ∂ i . In general characteristic, we replace the derivatives by the product by the "inverse" of the variables [START_REF] Mourrain | Multivariate Polynomials, Duality, and Structured Matrices[END_REF], [START_REF] Elkadi | Introduction à la résolution des systèmes polynomiaux[END_REF]. Conversely, if rank H Λ = r, by Proposition 2.12 Λ = r i=1 1 ζi • p i (∂) with a polynomial of degree 0, since the multiplicity of ζ i is 1. This concludes the proof of the equivalence.

In the binary case this also corresponds to the border rank of T * , therefore the r-th minors of the Hankel operator give equations for the r-th secant variety to the rational normal curves [START_REF] Iarrobino | Power sums, Gorenstein algebras, and determinantal loci[END_REF].

In order to compute the zeroes of an ideal I Λ when we know a basis of A Λ , we exploit the properties of the operators of multiplication in A Λ : M a : A Λ → A Λ , such that ∀b ∈ A Λ , M a (b) = a b and its transposed operator M t a : A * Λ → A * Λ , such that for ∀γ ∈ A * Λ , M ⊤ a (γ) = a ⋆ γ. The following proposition expresses a similar result, based on the properties of the duality. Proposition 2.15. For any linear form Λ ∈ R * such that rank H Λ < ∞ and any a ∈ A Λ , we have

(6) H a⋆Λ = M t a • H Λ Proof. By definition, ∀p ∈ R, H a⋆Λ (p) = a p ⋆ Λ = a ⋆ (p ⋆ Λ) = M ⊤ a • H Λ (p).
We have the following well-known theorem: Theorem 2.16 ([27, 26, 31]). Assume that A Λ is a finite dimensional vector space. Then

Λ = d i=1 1 ζi • p i (∂) for ζ i ∈ K n and p i (∂) ∈ K[∂ 1 , . . . , ∂ n ] and
• the eigenvalues of the operators M a and M t a , are given by {a(ζ 1 ), . . . , a(ζ r )}. • the common eigenvectors of the operators (M t xi ) 1≤i≤n are (up to scalar) 1 ζi . Using the previous proposition, one can recover the points ζ i ∈ K n by eigenvector computation as follows. Assume that B ⊂ R with |B| = rank(H Λ ), then equation ( 6) and its transposition yield

H B a⋆Λ = M t a H B Λ = H B Λ M a
, where M a is the matrix of multiplication by a in the basis B of A Λ . By Theorem 2.16, the common solutions of the generalized eigenvalue problem [START_REF] Bernardi | Computing symmetric rank for symmetric tensors[END_REF] (

H a⋆Λ -λ H Λ )v = O
for all a ∈ R, yield the common eigenvectors H B Λ v of M t a , that is the evaluation 1 ζi at the roots. Therefore, these common eigenvectors H B Λ v are up to a scalar, the vectors

[b 1 (ζ i ), . . . , b r (ζ i )] (i = 1, . . . ,

r).

Notice that it is sufficient to compute the common eigenvectors of (H xi⋆Λ , H Λ ) for i = 1, . . . , n If Λ = d i=1 γ i 1 ζi (γ i = 0), then the roots are simple, and one eigenvector computation is enough: for any a ∈ R, M a is diagonalizable and the generalized eigenvectors H B Λ v are, up to a scalar, the evaluation 1 ζi at the roots.

Coming back to our problem of partially symmetric tensor decomposition, T * ∈ R * δ1,...,δ k admits an affine decomposition of rank r iff T * coincide on R δ1,...,δ k with

Λ = r i=1 γ i 1 ζi , for some distinct ζ 1 , . . . , ζ r ∈ K n1 × • • • × K n k and some γ i ∈ K -{0}.
Then, by theorem 2.14, H Λ is of rank r and I Λ is radical.

Conversely, given H Λ of rank r with I Λ radical which coincides on R δ1,...,δ k with T * , by proposition 2.12, Λ = r i=1 γ i 1 ζi and extends T * , which thus admits an affine decomposition.

Therefore we can say that if the border rank of T is r then also rank(H Λ ) = r. Conversely if rank(H Λ ) = r, we can only claim that the border rank of T is at least r.

We say that Λ ∈ R * extends

T * ∈ R * δ1,...,δ k , if Λ * |R δ 1 ,...,δ k = T * .
The problem of decomposition of T * can then be reformulated as follows:

Truncated moment problem. Given T * ∈ R * δ1,...,δ k , find the smallest r such that there exists Λ ∈ R * which extends T * with H Λ of rank r and I Λ a radical ideal.

In the next section, we will describe an algorithm to solve the truncated moment problem.

Algorithm

In this section, we first describe the algorithm from a geometric point of view and the algebraic computation it induces). Then we characterize which conditions T * can be extended to Λ ∈ R * with H Λ is of rank r. The algorithm is described in 3.1. It extends the one in [START_REF] Brachat | Symmetric tensor decomposition[END_REF] which applies only for symmetric tensors. The approach used in [START_REF] Bernardi | Computing symmetric rank for symmetric tensors[END_REF] for the rank of tensors in Ξ 2 (S d (E)) and in Ξ 3 (S d (E)) allows to avoid to loop again at step 4: if one doesn't get simple roots, then it is possible to use other techniques to compute the rank. Unfortunately the mathematical knowledge on the stratification by rank of secant varieties is nowadays not complete, hence the techniques developped in [START_REF] Bernardi | Computing symmetric rank for symmetric tensors[END_REF] cannot be used to improve algorithms for higher border ranks yet.

Algorithm 3.1: Decomposition algorithm

Input: a tensor T ∈ S δ1 (E 1 ) ⊗ • • • ⊗ S δ k (E k ).
Output: a minimal decomposition of T . Set r = 0;

(1) Determine if T * can be extended to Λ ∈ R * with rank H Λ = r;

(2) Find if there exists r distinct points P 1 , . . . ,

P s ∈ Ξ(S δ1 (E 1 ) ⊗ • • • ⊗ S δ k (E k )
) such that T ∈ P 1 , . . . , P s ≃ P s-1 ; Equivalently compute the roots of ker H Λ by generalized eigenvector computation [START_REF] Bernardi | Computing symmetric rank for symmetric tensors[END_REF] and check that the eigenspaces are simple; (3) If the answer to 2 is YES, then it means that

T ∈ Ξ o r (S δ1 (E 1 ) ⊗ • • • ⊗ S δ k (E k )) \ Ξ r-1 (S δ1 (E 1 ) ⊗ • • • ⊗ S δ k (E k ))
; therefore the rank of T is actually r and we are done; (4) If the answer to 2 is NO, then it means that

T ∈ Ξ o r (S δ1 (E 1 ) ⊗ • • • ⊗ S δ k (E k )
) hence its rank is bigger than r; Repeat this procedure from step 2 with r + 1.

We are going to characterize now under which conditions T * can be extended to Λ ∈ R * with H Λ of rank r (step 1).

We need the following technical property on the bases of A Λ , that we will consider: Definition 3.1. Let B be a subset of monomials in R. We say that B is connected to 1 if ∀m ∈ B either m = 1 or there exists i ∈ [1, n] and m ′ ∈ B such that m = x i m ′ .

Let B, B ′ ⊂ R δ1,...,δ k be a two sets of monomials connected to 1. We consider the formal Hankel matrix

H B,B ′ Λ = (h α+β ) α∈B ′ ,β∈B , with h α = T * (x α ) = c α if x α ∈ R δ1,...,
:= (H B,B ′ Λ ) -1 H B,B ′ x i,l ⋆Λ (h)
for every variable x i,l ∈ R. We use the following theorems which extend the results of [START_REF] Laurent | A Sparse Flat Extension Theorem for Moment Matrices[END_REF] to the cases of distinct sets of monomials indexing the rows and columns of the Hankel operators. They characterizes the cases where

K[x] = B ⊕ I Λ :
Theorem 3.2. Let B = {x β1 , . . . , x βr } and B ′ = {x β ′ 1 , . . ., x β ′ r } be two sets of monomials of in R δ1,...,δ k , connected to 1 and let Λ be a linear form that belongs to

( B ′ • B + δ1,...,δ k ) * . Let Λ(h) be the linear form of B ′ • B + * defined by Λ(h)(x α ) = Λ(x α ) if x α ∈ R δ1,...,δ k and h α ∈ K otherwise. Then, Λ(h) admits an extension Λ ∈ R * such that H Λ is of rank r with B and B ′ basis of A Λ iff (8) M B i,l (h) • M B j,q (h) -M B j,q (h) • M B i,l (h) = 0 (0 ≤ l, q ≤ k, 1 ≤ i ≤ n l , 1 ≤ j ≤ n q ) and det(H B ′ ,B Λ(h) ) = 0. Moreover, such a Λ is unique.
Proof. If there exists Λ ∈ R * which extends Λ(h), with H Λ of rank r and B and B ′ basis of A Λ then H B ′ ,B Λ(h) is invertible and the tables of multiplications by the variables x i,l :

M i,l := (H B ′ ,B Λ(h) ) -1 H B ′ ,B
x i,l ⋆Λ(h) (Proposition 2.15) commute. Conversely suppose that these matrices commute and consider them as linear operators on B . Then by [START_REF] Mourrain | A new criterion for normal form algorithms[END_REF], we have such a decomposition R = B ⊕ I where I is an ideal of R. As a matter of fact, using commutation relation and the fact that B is connected to 1, one can easily prove that the following morphism:

π : R -→ B p → p(M)(1)
is a projection on B whose kernel is an ideal I of R (note that for any p ∈ B , p(M )(1) = p).

We define Λ ∈ R * as follows: ∀p ∈ R, Λ(p) = Λ(p(M)(1)) where p(M) is the operator obtained by substitution of the variables x i,l by the commuting operators M i,l . Notice that p(M) is also the operator of multiplication by p modulo I.

Let us prove by induction on the degree of b ′ ∈ B ′ that for all b ∈ B :

(9) Λ(h)(b ′ b) = Λ(h)(b ′ (M)(b))
and thus by linearity that

(10) Λ(h)(b ′ p) = Λ(h)(b ′ (M)(p))
for all p ∈ B . The property is obviously true for b ′ = 1. Suppose now that b ′ ∈ B ′ is a monomial of degree strictly greater than zero. As B ′ is connected to 1, one has b ′ = x j,q b ′′ for some variable x j,q ∈ R and some element b ′′ ∈ B ′ of degree smaller than the degree of b ′ . By construction of the operators of multiplication (M i,l ), we have 9) is proved. Let us deduce now that Λ extends Λ(h) i.e that for all b + ∈ B + ans b ′ ∈ B ′ we have:

Λ(h)(b ′ b) = Λ(h)(b ′′ x j,q b) = Λ(h)(b ′′ M j,q (b)). Finally we have that Λ(h)(b ′ b) = Λ(h)(b ′ (M)(b)) and (
Λ(b ′ b + ) := Λ(h)((b ′ b + )(M)(1)) = Λ(h)(b ′ b + ).
Indeed, from [START_REF] Bro | Parafac, tutorial and applications[END_REF] we have: 1) belongs to B . Then, by definition of multiplication operators (M i,l ) we have

Λ(h)((b ′ b + )(M)(1)) = Λ(h)((b ′ (M) b + (M)(1)) = Λ(h)(b ′ [b + (M)(1)]) as b + (M)(
Λ(h)(b ′ [b + (M)(1)]) = Λ(h)(b ′ b + ).
Thus, we have [START_REF] Buczynski | Determinental equations for secant varieties and the Eisenbud-Koh-Stillman conjecture[END_REF] Λ(b

′ b + ) = Λ(h)(b ′ b + ) for all b + ∈ B + ans b ′ ∈ B ′ (i.e Λ extends Λ(h)).
We eventually need to prove that I Λ = I := Ker(π). By the definition of Λ we obviously have that I ⊂ I Λ. Let us prove that I Λ ⊂ I: assume p belongs to I Λ, then from ( 11) and thus coincides with ker H Λ. As Λ ′ coincides with Λ on B, the two elements of R * must be equal. This ends the proof of the theorem.

Λ(b ′ p(M)(1)) = Λ(h)(b ′ p(M)(1)) = 0 for all b ′ ∈ B ′ . As p(M)(1) ∈ B and det(H B ′ ,B Λ )(h) = 0,
The degree of these commutation relations is at most 2 in the coefficients of the multiplications matrices M i,l . A direct computation yields the following, for m ∈ B:

• If x i,l , m ∈ B, and

x j,q m ∈ B then (M B i,l • M B j,q -M B j,q • M B i,l )(m) ≡ 0 in K(h). • If x i,l m ∈ B, x j,q m ∈ B then (M B i,l • M B j,q -M B j,q • M B i,l )(m) is of degree 1 in the coefficients of M i,l , M j,q . • If x i,l m ∈ B, x j,q m ∈ B then (M B i,l • M B j,q -M B j,q • M B i,l )(m) is of degree 2 in the coefficients of M i,l , M j,q .
We are going to give an equivalent characterization of the extension property, based on rank conditions:

Theorem 3.3. Let B = {x β1 , . . . , x βr } and B ′ = {x β ′ 1 , . . ., x β ′
r } be two sets of monomials in R δ1,...,δ k , connected to 1. Let Λ be a linear form in ( B ′+ * B + δ1,...,δ k ) * and Λ(h) be the linear form of

B ′+ • B + * defined by Λ(h)(x α ) = Λ(x α ) if x α ∈ R δ1,...,δ k and h α ∈ K otherwise. Then, Λ(h) admits an extension Λ ∈ R * such that H Λ is of rank r with B and B ′ basis of A Λ iff all (r + 1) × (r + 1)-minors of H B ′+ ,B + Λ(h) vanish and det(H B ′ ,B Λ(h) ) = 0.
Proof. First, if such a Λ exists then H B ′ ,B Λ(h) is invertible and H Λ is of rank r. Thus all the (r + 1) × (r + 1)-minors of H B ′+ ,B + Λ(h) are equal to zero.

Reciprocally, assume all (r + 1) × (r + 1)-minors of H B ′+ ,B + Λ(h) vanish and det(H B ′ ,B Λ(h) ) = 0 then one can consider the same operators: ⋆Λ(h) . By definition of these operators one has that

M B,B ′ i,l (h) := (H B,B ′ Λ(h) ) -1 H B,B ′ x i,l
(12) Λ(h)(x i,l b b ′ ) = Λ(h)(M i,l (b) b ′ ) for all b ∈ B and b ′ ∈ B ′ . As the rank of H B ′+ ,B + Λ(h)
is equal to the rank of H B ′ ,B Λ(h) we easily deduce that ( 12) is also true for b ′ ∈ B ′+ . Thus we have

Λ(h)(x j,q x i,l b b ′ ) = Λ(h)(x i,l b b ′ x j,q ) = Λ(h)(M i,l (b) b ′ x j,q ) = = Λ(h)(x j,q M i,l (b) b ′ ) = Λ(h)([M j,q M i,l ](b) b ′ ) for all b ∈ B, b ′ in B ′ and x i,l , x j,q ∈ R. Then, Λ(h)([M j,q M i,l ](b) b ′ ) = Λ(h)([M i,l M j,q ](b) b ′ ) = Λ(h)(x j,q x i,l b b ′ ) for all b ∈ B and b ′ ∈ B ′ . As H B ′ ,B Λ(h) is invertible we deduce [M j,q M i,l ](b) = [M i,l M j,q ](b) for all b ∈ B. Thus M j,q M i,l = M i,l M j,q .
Finally, we conclude the proof by using Theorem 3.2 Proposition 3.4. Let B = {x β1 , . . . , x βr } and B ′ = {x β ′ 1 , . . ., x β ′ r } be two sets of monomials of in R δ1,...,δ k , connected to 1 and Λ be a linear form on B ′ + * B + . Then, Λ admits an extension Λ in R * such that H Λ is of rank r with B (resp. B') a basis of A Λ iff

(13) H + = H G ′ G t J , with H + = H B ′ + ,B + Λ H = H B ′ ,B Λ and (14)
G = H t W, G ′ = H W ′ , J = W t H W ′ . for some matrix W ∈ K B×∂B ′ , W ′ ∈ K B ′ ×∂B . Proof. According to Theorem 3.3, Λ ∈ B ′ + * B + * admits a (unique) extension Λ ∈ R * such that H Λ is of rank r with B (resp. B') a basis of A Λ, iff H B ′ + ,B + Λ = H B ′ + ,B + Λ is of rank r.
Let us decompose H + Λ as [START_REF] Cardoso | Blind signal separation: statistical principles[END_REF] with

H + = H B ′ + ,B + Λ , H = H B ′ ,B Λ . If we have G = H t W, G ′ = H W ′ , J = W t H W ′ , then H + = H H W ′ W t H W t H W ′ is clearly of rank ≤ rank H.
Conversely, suppose that rank H + = rank H. This implies that the image of G ′ is in the image of H. Thus, there exists

W ′ ∈ K B ′ ×∂B such that G ′ = H W ′ . Similarly, there exists W ∈ K B×∂B ′ such that G = H t W. Thus, the kernel of (H G ′ ) (resp. (H t G)) contains W ′ -I (resp. W -I ). As rank H = rank H + = r, the kernel of (H G ′ ) (resp. (H t G)) is the kernel of H + (resp H +t ). Thus we have J = G t W ′ = W t H W ′ .
Notice that if H is invertible, W ′ , W are uniquely determined.

Introducing new variables w, w ′ for the coefficients of the matrices W, W, solving the linear system G = H t W, G ′ = H W ′ , and reporting the solutions in the equation J = W t H W ′ , we obtain a new set of equations, bilinear in w, w ′ , which characterize the existence of an extension Λ on R * . This leads to the following system in the variables h and the coefficients w of matrix W. It characterizes the linear forms Λ ∈ R * δ1,...,δ k that admit an extension Λ ∈ R * such that H Λ is of rank r with B a basis of A Λ.

(15) [START_REF] Mourrain | Multivariate Polynomials, Duality, and Structured Matrices[END_REF], whose structure is imposed by equality (linear) constraints on its entries. If H is known (i.e. B × B ⊂ R δ1,...,δ k , the number of independent parameters in H B,B + Λ (h) or in W is the number of monomials in B × ∂B -R δ1,...,δ k . By Proposition 3.4, the rank condition is equivalent to the quadratic relations J -W t H t W = 0 in these unknowns.

H B,∂B Λ (h) -H B Λ (h) W(w) = 0, H ∂B,∂B Λ (h) -W t (w) H B Λ (h) W(w) = 0 with det(H B Λ (h)) = 0. The matrix H B + Λ is a quasi-Hankel matrix
If H is not completely known, the number of parameters in H is the number of monomials in B × B -R δ1,...,δ k . number of independent parameters in H B,∂B

Λ (h) or in W is then B × ∂B -R δ1,...,δ k .
The system ( 15) is composed of linear equations deduced from quasi-Hankel structure, quadratic relations for the entries in B × ∂B and cubic relations for the entries in B × ∂B in the unknown parameters h and w.

We are going to use explicitly these characterizations in the new algorithm we propose for minimal tensor decomposition.

Examples and applications

There exist numerous fields in which decomposing a tensor into a sum of rank-one terms is useful. These fields range from arithmetic complexity [START_REF] Bürgisser | of Grundlehren der Mathematischen Wissenschaften[END_REF] to chemistry [START_REF] Smilde | Multi-Way Analysis[END_REF]. One nice application is worth to be emphasized, namely wireless transmissions [START_REF] Sidiropoulos | Blind PARAFAC receivers for DS-CDMA systems[END_REF]: one or several signals are wished to be extracted form noisy measurements, received on an array of sensors and disturbed by interferences. The approach is deterministic, which makes the difference compared to approaches based on data statistics [START_REF] Comon | Handbook of Blind Source Separation, Independent Component Analysis and Applications[END_REF]. The array of sensors is composed of J subarrays, each containing I sensors. Subarrays do not need to be disjoint, but must be deduced from each other by a translation in space. If the transmission is narrow band and in the far field, then the measurements at time sample t recorded on sensor i of subarray j take the form:

T (i, j, t) = r p=1 A ip B jp C tp if r waves impinge on the array. Matrices A and B characterize the geometry of the array (subarray and translations), whereas matrix C contains the signals received on the array. An example with (I, J) = (4, 4) is given in Figure 1. Computing the decomposition of tensor T allows to extract signals of interest as well as interferences, all included in matrix C. Radiating sources can also be localized with the help of matrix A if the exact location of sensors of a subarray are known. Note that this framework applies in radar, sonar or telecommunications.

4.1. Best approximation of lower multilinear rank. By considering a kth order tensor as a linear map from one linear space onto the tensor product of the others, one can define the ith mode rank, which is nothing else but the rank of that linear operator. Since there are k distinct possibilities to build such a linear operator, one defines a k-uplet of ranks (r 1 , . . . r k ), called the multilinear rank of the kth order tensor. It is known that tensor rank is bounded below by all mode ranks r i :

(16) r ≥ r i , ∀1 ≤ i ≤ k
This inequality gives us an easily accessible lower bound. Let's turn now to an upper bound.

Proposition 4.1.

[3] The rank of a tensor of order 3 and dimensions n 1 × n 2 × n 3 , with n 1 ≤ n 2 , is bounded by

(17) r ≤ n 1 + n 2 ⌊n 3 /2⌋
This bound on maximal rank has not been proved to be always reached, and it is likely to be quite loose for large values of n i . Nevertheless, it is sufficient for our reasoning.

There are two issues to address. First, the algorithm we have proposed is not usable in large dimensions (e.g. significantly larger than 10). The idea is then to reduce dimensions n i down to r i before executing the algorithm, if necessary. Second, another problem in practice is the presence of measurement errors or modeling inaccuracies, which increase the tensor rank to its generic value. We do not know how to reduce tensor rank back to its exact value. The practical solution is then to compute the best approximate of lower multilinear rank (r 1 , . . . r k ), as explained in [START_REF] Comon | Tensor decompositions, alternating least squares and other tales[END_REF]. This best approximate always exists, and inequality [START_REF] Chiantini | Weakly defective varieties[END_REF] shows that reducing dimensions will indirectly reduce tensor rank. To compute it, it suffices to minimize ||T -(U (1) , U (2) , U (3) ) • C|| with respect to the three matrices U (i) , each of size n i × r i , under the constraint U (1)H U (1) = I. If properly initialized by a truncated HOSVD, a few iterations of any iterative algorithm will do it [START_REF] Lathauwer | On the best rank-1 and rank-(R1,R2,. . . RN) approximation of high-order tensors[END_REF]. The tensor of reduced dimensions is then given by C = (U (1)H , U (2)H , U (3)H ) • T . 4.2. Number of solutions. In the above mentioned applications, it is necessary to have either a unique solution, or a finite set of solutions from which the most realistic one can be picked up. For this reason, it is convenient to make sure that the tensor rank is not too large, as pointed out by the following propositions.

Proposition 4.2. [17]

A generic symmetric tensor of order k ≥ 3 and rank r admits a finite number of decompositions into a sum of rank one terms if r < r E (k, n), where:

(18) r E (k, n) = n+k-1 d n
Rank r E is usually referred to as the expected rank of order k and dimension n.

Note that this result is true for generic tensors of rank r, which means that there exists a set of exceptions, of null measure. This proposition has not yet been entirely extended to unconstrained tensors, which we are interested in. However, some partial results are available in the literature [START_REF] Abo | Induction for secant varieties of Segre varieties[END_REF], and the following conjecture is generally admitted Conjecture 4.3. A generic tensor of order k ≥ 3 and rank r admits a finite number of decompositions into a sum of rank one terms if r < r E (k, n), where:

(19) r E (k, n) = k i=1 n 1 1 -k + k i=1 n i where n i denote the dimensions, 1 ≤ i ≤ k
On the other hand, a sufficient condition for uniqueness has been proposed by Kruskal [START_REF]Three-way arrays: Rank and uniqueness of trilinear decompositions[END_REF], but the bound is more restrictive: Proposition 4.4. [START_REF]Three-way arrays: Rank and uniqueness of trilinear decompositions[END_REF] A tensor of order k ≥ 3 and rank r admits a finite number of decompositions into a sum of rank one terms if:

(20) r ≤ k i=1 κ i 2 -1
where κ i denote the so-called Kruskal's ranks of loading matrices, which generically equal the dimensions n i if the rank r is larger than the latter.

4.3.

Computer results. If we consider a 4 × 4 × 7 unconstrained tensor, it has an expected rank equal to 9, whereas Kruskal's bound generically equals 6. So it is interesting to consider a tensor with such dimensions but with rank 6 < r < 9. In such conditions, we expect that there are almost surely a finite number of solutions. This tensor would correspond to measurements received on the array depicted in Figure 1, if 7 time samples are recorded. In [START_REF] Chiantini | On generic identifiability of 3-tensors of small rank[END_REF] L. Chiantini and G. Ottaviani claim that a computer check shows that for a generic 4 × 4 × 7 tensor of rank 7, uniqueness of the decomposition holds. Our computer results have been obtained with 4 × 4 × 7 tensors of rank 7 randomly drawn according to a continuous probability distribution.

First we consider a 4 × 4 × 4 tensor whose affine representation is given by: T := 4 + 7 a 1 + 8 a 2 + 9 a

If we consider B

′ := (1, b 1 , b 2 , b 3 ) and B := (1, a 1 , a 2 , a 3 ), the corresponding matrix H B ′ ,B Λ is equal to H B ′ ,B Λ =     4
    t M B c3 =     0 3/2 0 -1/2 -2 -33/2 14 11/2 -2 -57/2 23 17/2 -2 3/2 2 -1/2    
whose eigenvalues are respectively (-1, 2, 4, 1), (-2, 4, 5, 1) and (-3, 2, 6, 1). The corresponding common eigenvectors are:

v 1 =     1 -1 -2 3     , v 2 =     1 2 2 2     , v 3 =     1 5 7 3     , v 4 =     1 1 1 1     .
We deduce that the coordinates (a 1 , a 2 , a 3 , b 1 , b 2 , b 3 , c 1 , c 2 , c 3 ) of the 4 points of evaluation are:

ζ 1 :=               -1 -2 3 ? ? ? -1 -2 -3      -1 -2 3 -1 -1 -1 -1 -2 -3       1 1 1 1 1 1 1 1 1              
.

Finally, we have to solve the following linear system in (γ 1 , γ 2 , γ 3 , γ 4 ):

T = γ 1 (1 + a 1 + a 2 + a 3 ) (1 + b 1 + b 2 + b 3 ) (1 + c 1 + c 2 + c 3 ) + γ 2 (1 -a 1 -2 a 2 + 3 a 3 ) (1 -b 1 -b 2 -b 3 ) (1 -c 1 -2 c 2 -3 c 3 ) + γ 3 (1 + 2 a 1 + 2 a 2 + 2 a 3 ) (1 + 2 b 1 + 2 b 2 + 3 b 3 ) (1 + 2 c 1 + 4 c 2 + 2 c 3 ) + γ 4 (1 + 5 a 1 + 7 a 2 + 3 a 3 ) (1 + 3 b 1 -4 b 2 + 8 b 3 ) (1 + 4 c 1 + 5 c 2 + 6 c 3 ), We get γ 1 = γ 2 = γ 3 = γ 4 = 1.
We consider such an example with 6 time samples, that is an element of R , several coefficients are unknown. Yet, as will see, they can be determined by exploiting the commutation relations, as follows.

4 ⊗ R 4 ⊗ R 6 : T := 1046 a 1 b 1 c 1 + 959 a 1 b 1 c 2 + 660 a 1 b 1 c 3 + 866 a 1 b 1 c 4 + 952 a 1 b 1 c 5 -1318 a 1 b 2 c 1 - 1222 a 1 b 2 c 2 -906 a 1 b 2 c 3 -1165 a 1 b 2 c 4 -1184 a 1 b 2 c 5 -153 a 1 b 3 c 1 +52 a 1 b 3 c 2 +353 a 1 b 3 c 3 + 354 a 1 b 3 c 4 + 585 a 1 b 3 c 5 + 852 a 2 b 1 c 1 + 833 a 2 b 1 c 2 + 718 a 2 b 1 c 3 + 903 a 2 b 1 c 4 + 828 a 2 b 1 c 5 - 1068 a 2 b 2 c 1 -1060 a 2 b 2 c 2 -992 a 2 b 2 c 3 -1224 a 2 b 2 c 4 -1026 a 2 b 2 c 5 +256 a 2 b 3 c 1 +468 a 2 b 3 c 2 + 668 a 2 b 3 c 3 +748 a 2 b 3 c 4 +1198 a 2 b 3 c 5 -614 a 3 b 1 c 1 -495 a 3 b 1 c 2 -276 a 3 b 1 c 3 -392 a 3 b 1 c 4 - 168 a 3 b 1 c 5 + 664 a 3 b 2 c 1 + 525 a 3 b 2 c 2 + 336 a 3 b 2 c 3 + 472 a 3 b 2 c 4 + 63 a 3 b 2 c 5 + 713 a 3 b 3 c 1 + 737 a 3 b 3 c 2 +791 a 3 b 3 c 3 +965 a 3 b 3 c 4 +674 a 3 b 3 c 5 -95 a 1 b 1 +88 a 1 b 2 +193 a 1 b 3 +320 a 1 c 1 + 285 a 1 c 2 +134 a 1 c 3 +188 a 1 c 4 +382 a 1 c 5 -29 a 2 b 1 -2 a 2 b 2 +198 a 2 b 3 +292 a 2 c 1 +269 a 2 c 2 + 138 a 2 c 3 +187 a 2 c 4 +406 a 2 c 5 +119 a 3 b 1 -139 a 3 b 2 +20 a 3 b 3 -222 a 3 c 1 -160 a 3 c 2 +32 a 3 c 3 + 9 a 3 c 4 -229 a 3 c 5 +122 b 1 c 1 +119 b 1 c 2 +112 b 1 c 3 +140 b 1 c 4 +108 b 1 c 5 -160 b 2 c 1 -163 b 2 c 2 - 176 b 2 c 3 -214 b 2 c 4 -
The columns H B ′ ,{m} are also known for m ∈ {b Using the first relation on b 3 , we can reduce a 1 b 3 , a 2 b 3 , a 3 b 3 and obtain 3 linear dependency relations between the monomials in B ∪ {a 2 1 , a 1 a 2 , a 1 a 3 , a 2 2 , a 2 a 3 , a 2 3 }. Using the commutation relations lcm(m1,m2) m1 N (m 1 ) -lcm(m1,m2) m2 N (m 2 ), for (m 1 , m 2 ) ∈ {(a 1 b 1 , a 2 b 1 ), (a 1 b 2 , a 2 b 2 ), (a 2 b 2 , a 3 b 2 )} where N (m i ) is the reduction of m i with respect to the prevision relations, we obtain 3 new linear dependency relations between the monomials in B ∪ {a 2 1 , a 1 a 2 , a 1 a 3 , a 2 2 , a 2 a 3 , a 2 3 }. From these 6 relations, we deduce the expression of the monomials in {a 2 1 , a 1 a 2 , a 1 a 3 , a 2 2 , a 2 a 3 , a In particular, we can compute the entries of H B,B Λ . Solving the system H B,B

  we also use the notation H B Λ and H B Λ . If B, B ′ are linearly independent, then H B,B ′ Λ is the matrix of H B,B ′ Λ in this basis {b 1 , . . . , b r } of B and the dual basis of B ′ in B ′ * . The catalecticant matrices of [36] correspond to the case where B and B ′ are respectively the set of monomials of degree ≤ k and ≤ d -k (k = 0, . . . , d).

  generated by ker H B + Λ . Proof. Let us first prove that {b 1 , . . . , b r } ∩ I Λ = {0}. Let p ∈ b 1 , . . . , b r ∩ I Λ . Then p = i p i b i with p i ∈ K and Λ(p b j ) = 0. The second equation implies that H B Λ • p = 0, where p = [p 1 , . . . , p r ] t ∈ K r . Since H B

  this yields a contradiction.Consequently, {b 1 ⋆ Λ, . . . , b r ⋆ Λ} span the image of H Λ . For any p ∈ R, it holds that p ⋆ Λ = r i=1 µ i (b i ⋆ Λ) for some µ 1 , . . . , µ r ∈ K. We deduce that p -r i=1 µ i b i ∈ I Λ . This yields the decomposition R = b 1 , . . . , b r ⊕ K, R = B ⊕ I Λ , and shows that b 1 , . . . , b r is a basis of A Λ . If 1 ∈ B , the ideal I Λ is generated by the relations x j b k -r i=1 µ j,k i b i ∈ I Λ . These are precisely in the kernel of H B + Λ . Proposition 2.12. If rank(H Λ ) = r < ∞, then A Λ is of dimension r over K and there exist ζ 1 , . . . , ζ d ∈ K n , where d ≤ r, and

Definition 2 . 13 .

 213 For T * ∈ R * δ1,...,δ k , we call generalized decomposition of T * a decomposition such that T * = d i=1 1 ζi • p i (∂) where the sum for i = 1, . . . , d of the dimensions of the vector spaces spanned by the inverse system generated by 1 ζi • p i (∂) is minimal. This minimal sum of dimensions is called the length of f . This definition extends the definition introduced in[START_REF] Iarrobino | Power sums, Gorenstein algebras, and determinantal loci[END_REF] for binary forms. The length of T * is the rank of the corresponding Hankel operator H Λ .Theorem 2.14. Let Λ ∈ R * such that Λ = r i=1 γ i 1 ζi with γ i = 0 and ζ i distinct points of K n , iff rank H Λ = r and I Λ is a radical ideal.Proof. If Λ = r i=1 γ i 1 ζi , with γ i = 0 and ζ i distinct points of K n . Let {e 1 , . . . , e r } be a family of interpolation polynomials at these points: e i (ζ j ) = 1 if i = j and 0 otherwise. Let I ζ be the ideal of polynomials which vanish at ζ 1 , . . . , ζ r . It is a radical ideal. We have clearly I ζ ⊂ I Λ . For any p ∈ I Λ , and i = 1, . . . , r, we have p ⋆ Λ(e i ) = Λ(p e i ) = p(ζ i ) = 0, which proves that I Λ = I ζ is a radical ideal. As the quotient A Λ is generated by the interpolation polynomials e 1 , . . . , e r , H Λ is of rank r.

  we deduce that p(M)(1) = 0 and that p belongs to I. Thus we have I Λ ⊂ I. Eventually, Λ extends Λ(h) with I Λ = I := Ker(π) and A Λ equal to R/I ≃ B which is a zero dimensional algebra of multiplicity r with basis B. If there exists another Λ ′ ∈ R * which extends Λ(h) ∈ B ′ • B + * with rank H Λ ′ = r, by proposition 2.11, ker H Λ ′ is generated by ker H B ′ ,B + Λ ′

Figure 1 .

 1 Figure 1. Array of 10 sensors decomposed into 4 subarrays of 4 sensors each.

  3 + 5 b 1 -2 b 2 + 11 b 3 + 6 c 1 + 8 c 2 + 6 c 3 + 21 a 1 b 1 + 28 a 2 b 1 + 11 a 3 b 1 -14 a 1 b 2 -21 a 2 b 2 -10 a 3 b 2 + 48 a 1 b 3 + 65 a 2 b 3 + 28 a 3 b 3 + 26 a 1 c 1 + 35 a 2 c 1 + 14 a 3 c 1 + 18 b 1 c 1 -10 b 2 c 1 + 40 b 3 c 1 + 36 a 1 c 2 + 48 a 2 c 2 + 18 a 3 c 2 + 26 b 1 c 2 -9 b 2 c 2 + 55 b 3 c 2 + 38 a 1 c 3 + 53 a 2 c 3 + 14 a 3 c 3 + 26 b 1 c 3 -16 b 2 c 3 + 58 b 3 c 3 + 68 a 1 b 1 c 1 + 91 a 2 b 1 c 1 + 48 a 3 b 1 c 1 -72 a 1 b 2 c 1 -105 a 2 b 2 c 1 -36 a 3 b 2 c 1 + 172 a 1 b 3 c 1 + 235 a 2 b 3 c 1 + 112 a 3 b 3 c 1 + 90 a 1 b 1 c 2 + 118 a 2 b 1 c 2 + 68 a 3 b 1 c 2 -85 a 1 b 2 c 2 -127 a 2 b 2 c 2 -37 a 3 b 2 c 2 + 223 a 1 b 3 c 2 + 301 a 2 b 3 c 2 + 151 a 3 b 3 c 2 + 96 a 1 b 1 c 3 + 129 a 2 b 1 c 3 + 72 a 3 b 1 c 3 -114 a 1 b 2 c 3 -165 a 2 b 2 c 3 -54 a 3 b 2 c 3 + 250 a 1 b 3 c 3 + 343 a 2 b 3 c 3 + 166 a 3 b 3 c 3 .

  117 b 2 c 5 + 31 b 3 c 1 + 57 b 3 c 2 + 65 b 3 c 3 + 73 b 3 c 4 + 196 b 3 c 5 -35 a 1 -21 a 2 + 54 a 3 -3 b 1 -3 b 2 + 24 b 3 + 50 c 1 + 46 c 2 + 20 c 3 + 29 c 4 + 63 c 5 -6. If we take B = {1, a 1 , a 2 , a 3 , b 1 , b 2 } and B ′ = {1, c 1 , c 2 , c 3 , c 4 , c 5 } we obtain the following known submatrix of H Λ :

  3 , a 1 b 1 , a 2 b 1 , a 3 b 1 , a 1 b 2 , a 2 b 2 , a 3 b 2 }.Thus we deduce the relations between these monomials and B by solving the systemH B ′ ,B Λ X = H B ′ ,{m} Λ .This yields the following relations inA Λ : b 3 ≡ -1. -0.02486 a 1 + 1.412 a 2 + 0.8530 a 3 -0.6116 b 1 + 0.3713 b 2 , a 1 b 1 ≡ -2. + 6.122 a 1 -3.304 a 2 +.6740 a 3 +.7901 b 1 -1.282 b 2 , a 2 b 1 ≡ -2.+4.298 a 1 -1.546 a 2 +1.364 a 3 +.5392 b 1 -1.655 b 2 , a 3 b 1 ≡ -2. -3.337 a 1 + 5.143 a 2 + 1.786 a 3 -2.291 b 1 + 1.699 b 2 , a 1 b 2 ≡ -2. + 0.03867 a 1 -0.1967 a 2 + 1.451 a 3 -2.049 b 1 + 3.756 b 2 , a 2 b 2 ≡ -2. + 3.652 a 1 -3.230 a 2 + .9425 a 3 -2.562 b 1 +4.198 b 2 , a 3 b 2 ≡ -2.+6.243 a 1 -7.808 a 2 -1.452 a 3 +5.980 b 1 +0.03646 b 2

.

  The eigenvectors of the transposed operator normalized so that the first coordinate is 1 are: They correspond to the vectors of evaluation of the monomial vector B at the roots of I Λ . Thus we known the coordinates a 1 , a 2 , a 3 , b 1 , b 2 of these roots. By expanding the polynomialγ 1 (1 + a 1 + a 2 + a 3 )) (1 + b 1 + b 2 + b 3 ) (1 + • • • ) + γ 2 (1 -a 1 -2 a 2 + 3 a 3 ) (1 -b 1 -b 2b 3 ) (1 + • • • ) + γ 3 (1 + 2 a 1 + 2 a 2 + 2 a 3 ) (1 + 2 b 1 + 2 b 2 + 3 b 3 ) (1 + • • • ) + γ 4 (1 + 5 a 1 + 7 a 2 + 3 a 3 ) (1 + 3 b 1 -4 b 2 + 8 b 3 ) (1 + • • • ) + γ 5 (1 + 8 a 1 + 6 a 2 -7 a 3 ) (1 + 4 b 1 -5 b 2 -3 b 3 ) (1 + • • • ) + γ 6 (1 + 3 a 1 + 4 a 2 -5 a 3 ) (1 -3 b 1 + 5 b 2 + 4 b 3 ) (1 + • • • )(where the • • • are terms linear in c i ) and identifying the coefficients of T which do not depend on c 1 , . . . , c 5 , we obtain a linear system in γ i , which unique solution is (2, -1, -2, 3, -5, -3). This allows us to compute the value Λ for any monomials in {a 1 , a 2 , a 3 , b 1 , b 2 , b 3 }.

  δ k and otherwise h α is a variable. The set of these new variables is denoted h.

	Suppose that H B,B ′ Λ ators	is invertible in K(h), then we define the formal multiplication oper-M B,B ′ i,l

  which is invertible. Thus, the rank is at least 6. Let us find if H Λ can be extended to a rank 6 Hankel matrix H Λ . If we look at H B ′+ ,B +

				
				
				
	46 285 269 -160 119 -163	  
	20 134 138	32	112 -176	 
				
	29 188 187	9	140 -214	 
	63 382 406 -229 108 -117
	Λ

  2 3 } as linear combinations of monomials in B: a 2 1 ≡ 12.08 a 1 -5.107 a 2 + .2232 a 3 -2.161 b 1 -2.038 b 2 -2., a 1 a 2 ≡ 8.972 a 1 -1.431 a 2 + 1.392 a 3 -3.680 b 1 -2.254 b 2 -2., a 1 a 3 ≡ -11.56 a 1 +9.209 a 2 +2.802 a 3 +1.737 b 1 +.8155 b 2 -2., a 2 2 ≡ -2. + 6.691 a 1 + 2.173 a 2 + 2.793 a 3 -5.811 b 1 -2.846 b 2 , a 2 a 3 ≡ -2. -11.87 a 1 + 9.468 a 2 +2.117 a 3 +3.262 b 1 +0.01989 b 2 , a 2 3 ≡ -2.+16.96 a 1 -8.603 a 2 +1.349 a 3 -6.351 b 1 -.3558 b 2 . Now, we are able to compute the matrix of multiplication by a 1 in B, which is obtained by reducing the monomials B •a 1 = {a 1 , a 2 1 , a 1 a 2 , a 1 a 3 , a 1 b 1 , a 1 b 2 } by the computed relations:

			0.0 -2.0	-2.0	-2.0	-2.0	-2.0
		 	1.0 12.08	8.972 -11.56 6.122	0.03867
	M a1 :=	    	0.0 -5.107 -1.431 9.209 -3.304 -0.1967 0.0 0.2232 1.392 2.802 0.6740 1.451
		
		 	0.0 -2.161 -3.680 1.737	0.7901 -2.049
			0.0 -2.038 -2.254 0.8155 -1.282	3.756


Stéphane Cordier
email: stephane.cordier@univ-orleans.fr

Brigitte Lucquin-Desreux
email: brigitte.lucquin@math.jussieu.fr

Simona Mancini
email: simona.mancini@math.jussieu.fr

A numerical test for smoothing effects of collision operators

Keywords: Kinetic model, focalization of a beam, collision operator (Boltzmann, Fokker-Planck, Lorentz), smoothing effects, propagation of singularity

come

INTRODUCTION

In this paper, we are interested in the evolution of a system of particles injected in a bounded region all with the same velocity and undergoing collisions with heavier particles present in the domain. When a particle (e.g. a photon) collides with a heavier particle (e.g. neutron or ions) we can consider that the collision is elastic: the velocity modulus does not vary (or equivalently, the kinetic energy is conserved); the velocity modulus can be treated as a parameter of the problem. On the other hand, the velocity direction is changed by the scattering event. This behavior is modeled by the so called Boltzmann-Lorentz collision operator (see [START_REF] Degond | The asymptotics of collision operators for two species of particles of disparate masses[END_REF]), which in the two-dimensional case reads:

Q BL f (θ) = S 1 B(θ -θ ′ ) (f (θ ′ ) -f (θ)) dθ ′ .
(1.1)

When collisions become grazing, (see [START_REF] Degond | The Fokker-Planck asymptotics of the Boltzmann collision operator in the Coulomb case[END_REF][START_REF] Desvillettes | On asymptotics of the Boltzmann equation when the collisions become grazing[END_REF][START_REF] Goudon | On Boltzmann equations and Fokker-Planck asymptotics: influence of grazing collisions[END_REF]) the Boltzmann-Lorentz operator converges to the Fokker-Planck-Lorentz (or Laplace-Beltrami) operator (see [START_REF] Buet | The grazing collision limit for the Boltzmann-Lorentz model[END_REF][START_REF] Lucquin-Desreux | Diffusion of electrons by multicharged ions[END_REF]):

Q F P f (θ) = ∂ 2 θθ f. (1.2)
The focalization of a beam of particles is a test used in photonics. It consists in studying the evolution of a beam of photons injected in a bounded region from one side of the boundary with velocity close to the speed of light and perpendicular to the boundary. Inside the region, the photons may collide with neutrons or ions and change their trajectory, or may freely move from one side to the other one a straight line without changing their velocity direction. In this paper, we investigate the rate of particles reaching the opposite side of the domain (in the one-dimensional space case a slab) with a velocity direction equal to the incoming one (in a slab, perpendicular to the planes). It seems clear that the number of particles reaching the opposite side of the domain depends on the amount of collisions that they undergo. We call F the focalization coefficient, i.e. the rate of particles which exit the domain with the same velocity direction they had when entering the region (i.e. in the slab with a perpendicular velocity direction).

We will show how this simple model may be useful to check whether or not the cut-off Boltzmann equation leads to some kind of smoothing effect as the Fokker-Planck-Landau equation. This paper is organized as follows. In section 2 we introduce the transport equation modeling the evolution of a beam of particles. This model is one-dimensional in space and two-dimensional in velocity. In section 3 we discuss the smoothing effect of the collision operators (Boltzmann-Lorentz and Fokker-Planck-Lorentz) on singular data. In section 4 we present the numerical result validating our analysis. Finally, in section 5, we summarize our results and discuss those that would be obtained by other numerical approximations.

THE MODEL

Let us consider a system of particles with velocities v of modulus |v| = 1 and which undergo elastic collisions. The evolution of this system of particles is described by the distribution function f = f (x, θ, t) representing the number of particles which, at time t > 0, are in a position x ∈ [0, 1], with velocity v = (cos θ, sin θ), θ ∈ /2π¡ . The distribution function satisfies an equation of the form:

∂ t f + cos θ ∂ x f = 1 τ Q(f ), (2.3)
where τ is the collision relaxation time, and Q is either the isotropic Boltzmann-Lorentz operator:

Q BL (f )(θ) = < f > -f , < f > = 1 2π 2π 0 f (x, θ, t) dθ, (2.4)
or the Fokker-Planck-Lorentz (or Laplace-Beltrami) operator (1.2). Note that this problem is one-dimensional in x and two-dimensional in v. The length of the domain where the evolution of the system of particles occurs is rescaled to one (this is related to a change of collision time as explained in [START_REF] Cordier | Numerical approximation of the Vlasov-Fokker-Planck-Lorentz model[END_REF]).

We complete equation (2.3) with the initial data f (t = 0) = 0 and the following incoming boundary conditions:

f (x = 0, θ, t) = g(θ, t) , for θ ∈ [-π/2, π/2] f (x = 1, θ, t) = 0 , for θ ∈ [π/2, 3π/2].
(2.5)

These boundary conditions model an entering beam of particles on the left hand side (x = 0) of the slab and no re-entering particle flux on the right hand side (x = 1). We assume that the beam is well focalized with a velocity perpendicular to the region where collisions occur, the entering profile g can be considered as a Dirac mass along the direction θ = 0

g(θ, t) = δ θ=0 , ∀t > 0. (2.6)
We note that if there were no collision (i.e. τ → ∞), every particle would be transmitted, i.e. it would cross the domain in a time T = 1 and exit with the same velocity direction θ = 0. In other words, all the particles in absence of collisions would travel along a straight line and their velocity would not vary (nor in modulus nor in direction).

On the other hand, if there are many collisions in the region (i.e. τ → 0), then the particles of the beam entering the domain are reflected, i.e. they exit the domain on the same side, x = 0. Moreover, a very small number of particles will reach the opposite boundary of the domain, i.e. x = 1, and an even smaller number of them will have the velocity direction θ = 0. In other words, the beam is mostly reflected and the number of particles reaching the opposite side of the domain with the right velocity direction is negligible.

Finally, for 0 < τ < ∞ , e.g. τ = O(1), some of the particles are transmitted, others are reflected. We say that a stationary solution is reached when the flux of outgoing particles (either transmitted or reflected) equals the flux of incoming particles.

Recently, this focalization problem has been studied in the 3-dimensional [START_REF] Cordier | Numerical approximation of the Vlasov-Fokker-Planck-Lorentz model[END_REF] and 2-dimensional [START_REF] Lucquin-Desreux | A finite element approximation of grazing collisions, to appear on Transp[END_REF] contexts. In both cases, the curve of the focalization coefficient is given as a function of the relaxation time τ and mainly concerns the case of the Fokker-Planck operator.

The dependence of the solution on singular data has been studied in the space homogeneous case, see for instance [START_REF] Desvillettes | On asymptotics of the Boltzmann equation when the collisions become grazing[END_REF][START_REF] Villani | Regularity estimates via the entropy dissipation for the spatially homogeneous Boltzmann equation without cut-off[END_REF]. Concerning the space inhomogeneous case, propagation of singularities has been proven under the cut-off assumption in [START_REF] Boudin | On the singularities of the global small solutions of the full Boltzmann equation[END_REF]; one can also find some regularity results for the non cut-off case for a particular Boltzmann type operator in [START_REF] Desvillettes | On a model Boltzmann equation without angular cutoff[END_REF].

The aim of this paper is to try to characterize at a numerical level, by use of an efficient and simple test case, these regularizing (or non regularizing) properties.

SMOOTHING EFFECTS

We expect the solution to keep in time its singularity in the Boltzmann case. On the contrary, in the Fokker-Planck case, the solution is instantaneously regularized. As a matter of fact, let us first consider the space-homogeneous problem ∂ t f = Qf , with f = f (θ, t) independent on x, and initial data given by: f (θ, 0) = δ θ=0

(by analogy with the boundary conditions (2.5),(2.6) at time t = 0). The solution in the Boltzmann-Lorentz case (i.e. Q = Q BL with τ = 1) is given by:

f (θ, t) = δ θ=0 exp(-t) + (1 -exp(-t))1/2π, (3.7)
where we have used the fact that ∂ t < f >= 0, so that < f >= 1. It is easy to see that this solution converges towards a constant equilibrium state, but a singular, measured value part remains for any finite time (with an exponential decay with respect to time).

In the Laplace-Beltrami case (i.e. Q = Q F P with τ = 1), the solution of the space-homogeneous problem is the elementary solution of the heat equation with periodic conditions:

f (θ, t) = 1 2π n∈ exp(-n 2 t) exp(inθ), t > 0.
In this case, the singularity of the Dirac at initial time disappears for any t > 0 and the solution is smooth with respect to θ, as seen on the exponential decay of the Fourier coefficients. On the other hand, it has been proven that the Fokker-Planck-Landau operator regularizes the solution as it is done for the heat equation (see [START_REF] Villani | Regularity estimates via the entropy dissipation for the spatially homogeneous Boltzmann equation without cut-off[END_REF]). In the non homogeneous case, the singularities are propagated along the characteristics in the cut-off Boltzmann case, i.e. when the cross section is integrable and the lost and gain terms can be separated (see [START_REF] Boudin | On the singularities of the global small solutions of the full Boltzmann equation[END_REF]); whereas, in the non cut-off case, some smoothing effect occurs for a particular Boltzmann type operator (see [START_REF] Desvillettes | On a model Boltzmann equation without angular cutoff[END_REF]).

Concerning problem (2.3), although no analytical solution in the nonhomogeneous case is known, it is believed that the solution will have the same properties than in the space-homogeneous case: the singularity will follow the characteristics in the Boltzmann-Lorentz case, whereas the solution will instantaneously become smooth in the Laplace-Beltrami case.

NUMERICAL METHOD AND RESULTS

From now on we assume that the relaxation time is τ = 1. Let us approximate the distribution function on a uniform grid both in space x = x i = i∆x, ∆x = 1/N x and i = 0...N x, and in velocity angle θ = θ j = -π + j∆θ, ∆θ = 2π/N θ and j = 0...N θ -1. The scheme is based on a time splitting:

• one time step for the transport equation using an explicit upwind scheme with time step restriction for preserving stability (and positiveness): ∆t = ∆x (we choose a CFL condition equal to 1 to avoid numerical dissipation)

• one time step for the collision part which is solved using either a classical second order quadrature formula in the Boltzmann case, or a second order finite difference scheme in the Fokker-Planck case. Note that one can also use the exact solution known for the Boltzmann case (see equation (3.7)). Moreover, this part is treated implicitly so that there is no time step condition on this part. We also refer to [START_REF] Cordier | Numerical approximation of the Vlasov-Fokker-Planck-Lorentz model[END_REF] and [START_REF] Lucquin-Desreux | A finite element approximation of grazing collisions, to appear on Transp[END_REF] for an implicit finite element approximation respectively in 2D and 3D (with respect to the velocity variable).

The time needed for the first particle entering the domain to exit (if no collision occurs) is T = L/|v| = 1. We remark that the solution reaches its stationary equilibrium for T = 10 (200 iterations). The number of time iterations to reach a given time is proportional to the number of points in the space discretization. We choose N x = 20 and we are interested with the limit N θ → ∞.

In Figure 1, we plot the outgoing distribution function against velocity. The initial Dirac mass is instantaneously smoothed by means of the Fokker-Planck operator, whereas in the case of the Boltzmann operator the distribution function is still picked. In order to differentiate the two cases (Boltzmann and Fokker-Planck), we refine the mesh with respect to θ. When doing this, we expect the focalization coefficient either to converge to 0 (in the Fokker-Planck case), or to a non-vanishing limit value (in the Boltzmann case). This is exactly what we observe in Figure 2, where we plot the focalization coefficient in terms of the number of iterations and for different numbers N θ of discretization points in θ. We observe that in the Boltzmann case (on the left) this coefficient tends to a finite non zero value when N θ → ∞ (θ = 0 being a grid point there are still particles in this direction even when refining the grid). Whereas it rapidly goes to zero in the Fokker-Planck case (on the right). This is due to the fact that, the beam diffusing on more and more neighboring grid points, the number of particles which velocity precisely remains in the direction θ = 0 decreases more and more and finally vanishes. In Figure 3, we also consider the case of Coulomb type operators, i.e. with a truncation around the value θ = 0. More precisely, we consider operators of the form (1.1), where the kernel B now depends on a small truncation parameter ε in the following way (χ denotes the characteristic function):

B(z) = B ε (z) = C α,ε χ {|z|>ε} |z| -α ,
for some α, α ≥ 0, where the constant C α,ε is given by: C α,ε = 1 for α ∈ [0, 3[, C α,ε = 1/ log ε for α = 3 and C α,ε = 1/ε 3-α for α > 3.

When ε goes to zero, the kernel of the Boltzmann-Lorentz operator has an integrable singularity for α < 1, which allows to separate in the integral the gain term and the loss term: this corresponds to the so-called Grad's assumption of "angular cut-off". On the other hand, for α ≥ 1, there is a non integrable singularity around θ = 0, but a Taylor expansion in the integral shows that there is a compensation between the gain and loss term (see [START_REF] Buet | The grazing collision limit for the Boltzmann-Lorentz model[END_REF] and [START_REF] Degond | The Fokker-Planck asymptotics of the Boltzmann collision operator in the Coulomb case[END_REF]): these values of the parameter α correspond to the "non cut-off" assumption. Moreover, for α ≥ 3, the Boltzmann-Lorentz operator Q ε BL associated with this kernel B ε converges (up to a fixed multiplicative constant) to Q F P when ε goes to zero (see [START_REF] Buet | The grazing collision limit for the Boltzmann-Lorentz model[END_REF]); the case α = 3 exactly corresponds to the Coulomb case.

We expect the focalization test to give the same results for the cutoff assumption (for example with α = 0.95 < 1) than in the (isotropic) Boltzmann case, i.e. that the focalization coefficient, which is precisely defined by:

F(t) = f (x = 1, θ = 0, t)
f (x = 0, θ = 0, t) converges to a fixed non zero value, when refining the grid with respect to the angle. On the contrary, we expect that this coefficient will go to zero in the "non cut-off" case (for example for α = 2 or α = 3), where an instantaneous regularization of the beam will occur, as for the Fokker-Plank-Lorentz case. These facts have been effectively confirmed by our numerical simulations, has shown in Figure 3, where in the last column we have quoted the ratio R between the focalization coefficient for N θ = 320, F 320 , and the same coefficient for N θ = 160, F 160 . We note that in the limit case α = 1 the focalization coefficient behaves like in the Boltzmann-Lorentz case; this is not quite surprising, since this corresponds to the limit case where the compensation between the gain and loss terms is the weaker one. We also remark a very sharp decrease of the ratio R around the "critical" value α = 1.

F for N θ = 320 N θ = 160 R = F 320 /
Let us finally remark that this focalization coefficient also naturally depends on the space (and time, since ∆t = ∆x) discretization. For example, by changing N x from 10 to 20, one obtains a result for N θ = 320 which differs from 2.24% in the isotropic (Boltzmann-Lorentz) case. But we also observe that the relative variation of the ratio R, which is in fact the "right" coefficient which allows to distinguish the regularizing case (i.e. the non cut-off case) from the Boltzmann one, is only of 0.04%. Also the relative error decreases when refining the mesh with respect to space (and time): by changing N x from 20 to 40, the relative error on the focalization coefficient with N θ = 320 is of 1.18% and on the ratio R it is of order 0.03%. These errors are a bit more important in the Coulomb case (α = 3): between the case N x = 20 and N x = 40, they are of the order 5.5% for the focalization coefficient with N θ = 320, and of order 3.08% for the ratio R, which is still quite reasonable.

CONCLUSION

Note that the smoothing effect can also be observed using other methods: Monte Carlo or spectral methods.

• Using stochastic or Monte Carlo methods, the distribution function is represented by pseudo-particles. The numerical method is based on a time splitting: one time step solving the transport part, one time step for collisions.

In the Boltzmann case (see [START_REF] Cercignani | The Boltzmann equation and its applications[END_REF][START_REF] Graham | Stochastic particle approximations for generalized Boltzmann models and convergence estimates[END_REF]), the method consists in performing a change of the particles velocity (choosing a post-collision velocity uniformly on the sphere) according to a random variable related to the collision time. It can be checked that the probability P of a particle exiting the domain without changing its velocity is not vanishing; more precisely, P = exp(-L/τ ), where L is the length of the domain, and τ is the collision time.

In the Fokker-Planck (or Laplace-Beltrami) case, the solution can be obtained either using the grazing collision limit, by making a large number of small deviations, or using that for small time steps. The solution of the collision part for a Dirac initial data behaves as a Gaussian with variance related to the time step. Thus, the velocities are chosen according to a Gaussian distribution. Moreover, the probability P of exiting the domain with the velocity of the beam is zero. This property can be used to separate the case of Boltzmann operator (P = 0) and of Fokker-Planck operator (P = 0).

• The same holds true in the spectral method (see [START_REF] Pareschi | Fast spectral methods for Boltzmann and Landau integral operators of gas and plasma kinetic theory[END_REF]). First, let us recall that the smoothness of a function is related to the decay of its Fourier coefficients (e.g. the Dirac measure corresponds to constant Fourier coefficients). When using a spectral method, the distribution function in velocity (or equivalently in angle) is represented by its Fourier coefficients. In the case of the Boltzmann operator, the distribution function at the right hand side has a singularity and thus the Fourier coefficients a k at x = L have a non zero limit when k goes to ∞. Concerning the Fokker-Planck (or Laplace-Beltrami) operator, the Fourier coefficients go to 0 as k → ∞ (exponentially fast). The difficulty with spectral methods is that the representation of Dirac masses requires a large number of modes.

The same differentiation should be observed on the classical non linear Boltzmann operator (for binary collision) -in which case the singularities follow the characteristics, see [START_REF] Boudin | On the singularities of the global small solutions of the full Boltzmann equation[END_REF] -and Fokker-Planck-Landau or non cutoff Boltzmann equation -in which case, regularity occurs, see [START_REF] Villani | Regularity estimates via the entropy dissipation for the spatially homogeneous Boltzmann equation without cut-off[END_REF] and [START_REF] Desvillettes | On a model Boltzmann equation without angular cutoff[END_REF]. Thus, this test serves to check smoothing effects, whatever the method is used for numerical simulations (discrete velocity, spectral or Monte Carlo methods).

Figure 1 :

 1 Figure 1: The outgoing profile f (x = 1, •, T = 10) versus the index j of the velocity angle θ j in the Boltzmann and Fokker-Planck case (N θ = 80)

Figure 2 :

 2 Figure 2: The focalization coefficient against the number of iterations in the Boltzmann-Lorentz case (left), and in the Fokker-Planck-Lorentz case (right)

Figure 3 :

 3 Figure 3: Focalization coefficient, N iter = 200

ACKNOWLEDGMENT

The authors thanks S. Méléard for helpful discussions, and acknowledge support from the TMR project "Asymptotic methods in kinetic theory" (TMR number: ERB FMRX CT97 0157), run by the European Community.


Hocine Menana
email: hocine.menana@univ-nantes.fr

Mouloud Féliachi
email: mouloud.feliachi@univ-nantes.fr

Modeling the response of a rotating eddy current sensor for the characterization of carbon fiber reinforced composites

Keywords: 82, 20, Wt Computational modeling; simulation -81, 70, Ex Nondestructive testing: electromagnetic testing, eddy-current testing

This work deals with the eddy current nondestructive testing of carbon fiber reinforced polymer composites (CFRPs), and the development of rapid numerical models for the simulation of the interaction between the electromagnetic field and complex structures. We show qualitatively the possibility of characterizing the CFRPs by using a rotating eddy current sensor. We make use of a numerical model which we have developed in a previous work for eddy current computation in CFRPs. This model is completed in this work to account for the impedance variation of the rotating eddy current sensor, used to determine the anisotropy factor and the fiber orientation in each ply of a layered CFRP plate.

Introduction

Several experimental works have highlighted the possibility of characterizing defects in CFRPs by using the eddy current nondestructive testing technique [1]- [2]. In this work which consists in modeling and simulations, we show qualitatively the possibility of characterizing the physical properties of layered CFRPs by using a rotating eddy current sensor. Modeling the electromagnetic field in layered CFRPs is a challenging problem due to their multi scale dimensions and their high anisotropy, associated to complex capacitive couplings due to the heterogeneity of their structures. The current calculation tools do not allow the modeling of the CFRPs on the scale of their structure; global anisotropic properties are thus assigned to them by using the homogenization techniques [3]. Semi analytical models based on the dyadic Green's functions have been proposed for the electromagnetic field modeling in anisotropic layered media [4]- [5]; these models are however limited to simple geometries. In this work, we make use of a numerical model which we have developed in a previous work for eddy current computation in CFRPs [6]. It is based on the T-φ formulation [7], where only the parts of the system carrying a source or an induced current are discretized; the inductor motion and the multiscale dimensions of the layered CFRPs are thus easily treated. This model is completed to account for the impedance variation of a rotating eddy current sensor, used to determine the anisotropy factor and the fiber orientation in each ply of a layered CFRP plate.

The modeled system

The modeled system is described in figure 1. It consists of a rectangular air cored coil of length l and width w, positioned at a distance h above a multilayered CFRP plate of surface L×W and thickness e. The coil is excited with a sinusoidal current source, and rotates around an axis perpendicular to the surface of the plate (figure 1). The structure of the considered CFRP is described in figure 2. These materials are made of several unidirectional plies of about 125µm thickness, stacked with different fiber orientations. Each unidirectional ply is constituted of carbon fibers embedded in an electrically non conductive matrix. In each ply, the electrical conductivity varies between 5×10 3 and 5×10 4 S/m in the direction of the fibers, and between 10 and 10 2 S/m in the transverse direction, depending on the type and the volume fraction of the fibers [8], [START_REF] Grimberg | Proceedings of the 8 th International Conference of the Slovenian Society for Non-Destructive Testing[END_REF]- [START_REF] Park | [END_REF]. When the plies are assembled, contacts between the fibers of adjacent plies give rise to a cross ply conductivity which is of lower value, being typically half the value of the transverse conductivity [8]. The homogenized conductivity tensor of a CFRP plate is given by (1), involving the longitudinal and transversal conductivities (σ L , σ T ) of the plies, the orientation (θ) of each ply, and the cross ply conductivity (σ cp ) [6,8].

                + - - + = cp T L T L T L T L σ θ σ θ σ θ σ σ θ σ σ θ σ θ σ σ 0 0 0 ) ( cos ) ( sin ) 2 sin( 2 0 ) 2 sin( 2 ) ( sin ) ( cos 2 2 2 2
(1)

The model equations

In the following, Ω S and Ω C denote respectively the volumes of the source (coil) and the CFRP plate regions. The latter is delimited by the frontier Γ. The source region is the support of a current density J s with a time harmonic variation. The electromagnetic problem is formulated in terms of the electric vector potential T and the scalar magnetic potential φ, governed by the equations ( 2)-( 3) in Ω C , involving the conductivity tensor σ of the CFRP plate, the free space permeability µ 0 and the angular frequency ω [6].

s H T T 0 0 1 ) ( ωµ ϕ ωµ j j σ - = ∇ - + × ∇ × ∇ - (2) 0 ) ( = ∇ - + ∇ ϕ T H . s
(3) The eddy currents density J in Ω C is given by

T J × ∇ = . ( 4
)
To cancel the normal component of the current density at the frontiers of Ω C , the following boundary condition has to be imposed

Γ = × / 0 T n . (5)
The source term H s in ( 2) is calculated by the Biot and Savart formulae given by ( 6), in which ||r-r'|| is the distance between a source point in Ω S and a point in Ω C where the magnetic field H s is evaluated.

Ω × = ∫ Ω d S 3 ) ( ) ( 4π 1 ) ( r' - r r' - r r' J r H s s (6)
Only the active parts of the system are discretized; the equations (2) (3) and ( 4) are discretized by using the finite difference method in Ω C . The source term H s is evaluated at each node of the finite difference mesh by mean of (6). The source region Ω S is discretized into Ns elementary volumes; the current source density J s is defined at the center of each elementary volume. The discrete form of the integral ( 6) is given as follows [6]

s s J H s M = . ( 7
)
Taking account of ( 7), the discrete forms of (2) (3) and ( 4) are given as follows [6] ( )

s J T s 0 0 0 M G I C R C ωµ ϕ ωµ ωµ j j j - = - + , (8)
( ) 0 G M D = - + ϕ T J s s , (9)
T J C = . ( 10
)
Details about the development of the matrices involved in the equations ( 7) to [START_REF] Park | [END_REF], as well as a validation of these equations are given in [6].

The impedance variation (∆Z) of the coil, due to the eddy currents flowing in the CFRP plate, is calculated by using equation (11) [4], where I s is the source current, and E i the electric field scattered by the CFRP plate, and calculated using the Biot and Savart formulae (12) once the eddy currents density J in Ω C is obtained.

∫ Ω Ω - = ∆ S d I Z s s J . E i 2 1 (11) ) ' ( , ' ) ' ( 4π ) ( 0 C C d j Ω ∈ Ω - - = ∫ Ω r r r r J r E ωµ i (12)
Introducing ( 12) in (11), we obtain

      Ω Ω - = ∆ Ω ∈ Ω ∈ ∫ ∫ Ω Ω C S C S s 2 s S C d d I j Z ' ' ) ' ( ) ( 4π 0 r r r r r J . r J ωµ . ( 13
)
The discrete form of ( 13) is given by ( 14), where vs n and vc m are elementary volumes of the source and the CFRP plate regions respectively, and ||r n -r m || is the distance between them as shown in figure 3. Nc is the number of the elementary volumes of the CFRP plate. Since the plies are very thin, the current density is assumed to be constant along the section of each ply; each ply is then represented by only one layer of finite difference mesh nodes. Each elementary volume vc m results from the surface defined by four adjacent finite difference mesh nodes (i,j,k; i,j+1,k; i+1,j,k; i+1,j+1,k) multiplied by the thickness (δp) of the ply number k. The current density J m in the element vc m is obtained by averaging the contribution of the four finite difference mesh nodes. ) )

∑∑

k j i k j i k j i k j i m , 1 , 1 , 1 , , , 1 , , 4 1 +
+ + + + + + = J J J J J (15) k j i k j i k j i k j i m p vc , , , 1 , , , , , 1 . r r
k j i k j i k j i k j i m , 1 , 1 , 1 , , , 1 , , 4 1 +
+ + + + + + = r r r r r (17)
The equation ( 14) can be written in a matrix form as follows

s J J Z M T = ∆ , ( 18
)
where

{ } { } { } { } { } { }           =           =           = z y x z s y s x s sc sc sc J J J J J J J J s , , M 0 0 0 M 0 0 0 M M . ( 19
)
The matrix sc M is of dimension (Ns×Nc); its elements are calculated as follows (20)

Numerical examples

In these numerical examples, we consider the system presented in figure 1. We show qualitatively the influences of the anisotropy factor, the fibers orientations (θ) in the CFRP plate and the coil design on the impedance variation of the latter. The numerical values of the modeled systems parameters are given in table1. Two CFRP plates are considered; the first one is constituted of one unidirectional ply of 125µm thickness and the second one is constituted of four unidirectional plies of 125µm thickness oriented in the directions (0°, 45°, 90°,-45°). The direction of the fibers of the upper ply is taken as reference (θ=0°). The rectangular coil, of 10 mm length, rotates above the CFRP plate with an angular step (∆α =3°); its impedance variation is recalculated at each position (α). First, we studied the influences of the plate anisotropy and the coil width on its impedance variation. We considered one unidirectional CFRP ply. Figure 4 represents the real and imaginary parts of the normalized impedance variation according to the coil rotation above the ply, for different values of the coil width (w). The impedance variation (∆Z) is in the form of lobes oriented in the direction of the fibers (θ=0°). We notice that increasing the coil width leads to a significant increase of ∆Z at (α =θ ±π) (coil oriented in the fibers direction), and a modest increase of ∆Z at (α =θ ±π/2) (coil oriented in the direction transverse to the fibers). The imaginary part is about only 2% of the total impedance variation. Figure 5 represents the real part of the normalized impedance variation according to the coil rotation, for different values of the transversal conductivity (σ T ) and a fixed value of the longitudinal conductivity (σ L ). We notice that increasing the transversal conductivity leads to a significant and non proportional increase of ∆Z at the position (α =θ ±π), and a modest proportional increase of ∆Z at the position (α =θ ±π/2). Figure 6 represents the real part of the normalized impedance variation according to the coil rotation, for different values of the longitudinal conductivity and a fixed value of the transversal conductivity. We notice that increasing the longitudinal conductivity leads to a significant and non proportional increase of ∆Z at (α =θ ±π), and a negligible increase of ∆Z at (α =θ ±π/2). Figures 5 and6 can be used to calibrate the sensor for the inverse problem which leads to the determination of the longitudinal and transversal conductivities from measurements of ∆Z at (α =θ ±π) and (α =θ ±π/2).

In the second step, we studied the influence of the coil width on the determination of the fiber orientations in a multi-layer CFRP plate. Figure 7 shows the normalized impedance variation according to the coil rotation above a CFRP plate constituted of four plies oriented in the directions (0°, 45°, 90°, -45°), for different widths (w) of the coil. When the coil width is too large (w=2mm), the influence of the fiber orientation in the lower plies on the impedance variation is hidden by the influence of the upper plies. On the other hand, a small width (w=1mm) leads to a weak interaction between the coil and the lower plies. The coil width must then satisfy a compromise between the sensitivity and the distinguishability. Figure 8 shows the eddy currents repartition in the four plies plate when the rectangular coil is oriented in the direction of the fibers of the first ply (α=θ=0°). The effect of the anisotropy on the eddy current repartition in the plies is noticeable; the eddy currents are stretched in the direction of the high conductivity.

Finally, to give an idea about the computing time, the latter is of about 30 min for 120 positions of the coil on the four plies plate discretized into (50×50×4) nodes. A PC of 3.4 GHz processor is used. The matrices in equations ( 7) to (18) are calculated only once; this makes the calculation fast.

Conclusion

We have shown the possibility of characterizing CFRPs by using a rotating eddy current sensor of a rectangular shape. Although the results in the treated examples are only qualitative, interesting phenomena have been observed, giving rise to new ideas which will be developed in future works. In the short run, we envisage the development of an inverse model associated to the rotating eddy current sensor for the electrical conductivity characterization in CFRPs, and applications to defects characterization such as delaminations and fiber breakings. An experimental validation is also necessary to determine the frequency limit of validity of the model equations. In the long run, we project to take account of the capacitive effects in the modelling.

Figures and tables captions

Fig. 1 .Fig. 2 .Fig. 3 .Fig. 4 .Fig. 5 .Fig. 6 .Fig. 7 .Fig. 8 .

 12345678 Fig. 1. The modeled system

  Figures

  Fig. 2

Fig. 4

 4 Fig. 4

Fig. 5

 5 Fig. 5

Fig. 6

 6 Fig. 6

Fig. 7

 7 Fig. 7

Fig. 8

 8 Fig. 8

Table 1 .

 1 The numerical values of the modeled system parameters

Acknowledgments

Work supported by Nantes Métropole.

Tables