Jean-Pierre Gazeau

Geometric study of the beta-integers for a Perron number and mathematical quasicrystals

R. Nous étudions géométriquement les ensembles de points de R obtenus par la betanumération que sont les β-entiers Z β ⊂ Z[β] o ù β est un nombre de Perron. Nous montrons qu'il existe deux schémas de coupe-et-projection canoniques associés à la β-numération, o ù les β-entiers se relèvent en certains points du réseau Z m (m = degré de β) , situés autour du sous-espace propre dominant de la matrice compagnon de β . Lorsque β est en particulier un nombre de Pisot, nous redonnons une preuve du fait que Z β est un ensemble de Meyer. Dans les espaces internes les fenêtres d'acceptation canoniques sont des fractals dont l'une est le fractal de Rauzy (à quasi-homothétie près). Nous le montrons sur un exemple. Nous montrons que Z β ∩ R + est de type fini sur N, faisons le lien avec la classification de Lagarias des ensembles de Delaunay et donnons une borne supérieure effective de l'entier q dans la relation :

A. We investigate in a geometrical way the point sets of R obtained by the β-numeration that are the β-integers Z β ⊂ Z[β] where β is a Perron number. We show that there exist two canonical cut-and-project schemes associated with the β-numeration, allowing to lift up the βintegers to some points of the lattice Z m (m = degree of β) lying about the dominant eigenspace of the companion matrix of β . When β is in particular a Pisot number, this framework gives another proof of the fact that Z β is a Meyer set. In the internal spaces, the canonical acceptance windows are fractals and one of them is the Rauzy fractal (up to quasi-dilation). We show it on an example. We show that Z β ∩R + is finitely generated over N and make a link with the classification of Delone sets proposed by Lagarias. Finally we give an effective upper bound for the integer q taking place in the relation: x, y ∈ Z β =⇒ x + y (respectively xy ) ∈ β -q Z β if x + y (respectively xy ) has a finite Rényi β-expansion.

Introduction

Gazeau [Gaz], Burdik et al [Bu] have shown how to construct a discrete set Z β ⊂ Z[β] ⊂ R which is a Delone set [Mo], called set of β-integers (or beta-integers), when β > 1 is a Pisot number of degree greater than 2. A beta-integer has by definition no fractional part in its Rényi β-expansion [Re] [Pa]. As basic feature, this Delone set is self-similar, namely βZ β ⊂ Z β .

2000 Mathematics Subject classification : 11A63, 11R06, 52C23. Keywords : quasicrystal, Perron number, Pisot number, Delone set, beta-numeration.

Since the general notion of β-expansion of real numbers (see section 2 for definitions) was created by Rényi for any real number β > 1, the set of beta-integers Z β , defined as the set of real numbers equal to the integer part of their β-development, is defined without ambiguity in full generality and is self-similar by construction: βZ β ⊂ Z β . The main questions we may address are the following: (Q1) For which β > 1 is Z β a Delone set ? or equivalently (Q1') for which β > 1 is Z β a uniformly discrete set ? since the sets Z β of beta-integers are always relatively dense by construction. Now Delone sets are classified into several types (see the definitions in the Appendix) so that the following question is also fundamental: (Q2) For which class of β > 1 is Z β a Delone set of a given type ?

The uniform discretness property of Z β is a crucial property which is not obtained for all real number β, but very few general results are known nowadays. Thurston has shown that it is the case when β is a Pisot number [Th]. It is conjectured that it is also the case when β is a Perron number. Apart from the Pisot case, many open questions remain (Bertrand-Matthis [START_REF] B -M | Questions diverses relatives aux systèmes codés: applications au θ-shift[END_REF], Blanchard [Bl]) and are expressed in terms of the β -shift. Schmeling [Sc] has proved that the class C 3 of real numbers β > 1 such that the Renyi-expansion d β (1) of 1 in base β contains bounded strings of zeros, but is not eventually periodic, has Hausdorff dimension 1. For all β in this class C 3 , the β-shift is specified [Bl]. It is obvious that the specification of the β-shift is equivalent to the fact that Z β is uniformly discrete. So that the class C 3 would contain all Perron numbers. The idea of exploring relationships between the β-shift and the algebraic properties of β in number theory is due to A. Bertrand-Matthis [START_REF] B -M | Nombres de Perron et questions de rationnalité[END_REF]. In this direction, some results are known (Akiyama [Ak] [Ak1]). Parry [Pa] has proved that the β-shift is sofic when β is a Pisot number. Lind [Li] conversely has shown that β is a Perron number if the β-shift is sofic.

In section 2 we will recall some basic facts about the β-numeration and the beta-integers.

In section 3, we will establish the geometrical framework which is attached to the algebraic construction of the set of the beta-integers when β is a Perron number in general (of degree m

2). Namely, by geometric framework, we mean that we will show the existence of two cut-and-project schemes (see the definitions in the Appendix) embedded in a canonical way in the Jordan real decomposition of R m where this decomposition is obtained by the action of the companion matrix of β, respectively of its adjoint, the second cut-and-project scheme being the dual of the first one. This will be done without invoking any substitution system on a finite alphabet [AI] or the theory of Perron-Frobenius [Mi]. These cut-and-project schemes will consist of an internal space which will be an hyperplane of R m complementary to a onedimensional line on which the set of β-integers will be set up in a natural way, together with the usual lattice Z m in R m . The constituting irreducible subspaces of the internal spaces will appear by construction as asymptotic linear invariants. This will allow us to deduce several results when β is a Pisot number: a minimal acceptance window in the internal space closely related to the Rauzy fractal, a geometrical proof that Z β is a Meyer set, the fact that Z β is finitely generated over N. We will make a link on an example with the Rauzy fractal when the beta-integers arise from substitution systems of Pisot type (for instance Rauzy [Ra], Arnoux and Ito [AI], Messaoudi [Me] [Me1], Ito and Sano [IS], Chap. 7 in Pytheas Fogg [PF]). At this point, we should outline that the main difference with the substitutive approach is that the matrices involved may have negative coefficients (compare with the general approach of Akiyama [Ak] [Ak1]).

The additive properties of Z β will be studied in section 4 by means of the canonical cutand-project schemes when β is a Pisot number: in A), we shall show that the elements of Z β ∩ R + can be generated over N by elements of Z β of small norm, in finite number, using truncated cones whose axis of revolution is the dominant eigenspace of the companion matrix of β and a Lemma of Lind [Li] on semigroups; in B), we will provide a geometrical interpretation of the maximal preperiod of the β-expansion of some real numbers coming from the addition of two beta-integers, of the finite sets T and T in the relations [Bu]

Z + β + Z + β ⊂ Z + β + T , Z + β -Z + β ⊂ Z β +
T and an upper bound of the integer q taking place in the relation x, y ∈ Z + β =⇒ x ± y ∈ β -q Z β when x + y and xy have finite β-expansions.

Beta-numeration and beta-integers

Let β ∈ (1, +∞) \ N. We will refer in the following to Rényi [Re], Parry [Pa] and Frougny [Fro] [Fro1] [Bu]. For all x ∈ R we will denote by x , resp. {x} = xx , the usual integer part of x, resp. its fractional part. Let us denote by T (x) = {βx} the ergodic transformation sending

[0, 1] into itself. For all x ∈ [0, 1] , the iterates T n (x) := T (T n-1 (x)), n 1, with T 0 := I d by convention, provide the sequence (x -i ) i 1 of digits, with x -i := βT i-1 (x) , in the finite alphabet A = {0, 1, • • • , β }. The element x is then equal to its Rényi β-expansion +∞ j=1 x -j β -j also denoted by 0.x -1 x -2 x -3 . . .. The Rényi β-expansion of 1 will be denoted by d β (1). The operator T on [0, 1] induces the shift σ : (x -1 , x -2 , . . .) → (x -2 , x -3 , . . .) on
the compact set A N (with the usual product topology). The closure of the subset of A N invariant under σ takes the name of β-shift. The knowledge of d β (1) suffices to exhaust all the elements in the β-shift (Parry [Pa]). For this let us define the following sequence (c i ) i 1 in

A N : c 1 c 2 c 3 • • • =    t 1 t 2 t 3 • • • if the Rényi β-expansion d β (1) = 0.t 1 t 2 • • • is infinite, t 1 t 2 • • • t r-1 (t r -1) ω if d β (1) is finite and equal to 0.t 1 t 2 • • • t r ,
where ( ) ω means that the word within ( ) is indefinitely repeated. Then the sequence (y -i ) i 1 in A N is exactly the sequence of digits provided by the iterates of y = +∞ i=1 y -i β -i by T n if and only if the following inequalities are satisfied: (y -n , y -(n+1) , . . .) < (c 1 , c 2 , c 3 , . . .) for all n 1 where " < " means lexicographical smaller. These inequalities will be called conditions of Parry. We will now use finite subsets of the β-shift.

D 2.1. -Let Z + β = {x k β k +x k-1 β k-1 +• • •+x 1 β+x 0 | x i ∈ A, k 0, and (x j , x j-1 , . . . , x 1 , x 0 , 0, 0, • • • ) < (c 1 , c 2 , • • • ) for all j, 0 j k } be the discrete subset of R + of the real numbers equal to the integer part of their Rényi β-expansion. The set Z β = Z + β ∪ -Z + β is called the set of β-integers. For all x ∈ R + , if x = p i=-∞ x i β i with p
0, is obtained by the greedy algorithm, then (x i ) i p will satisfy the conditions of Parry. We will denote by int(x) = p i=0 x i β i the integer part of its Rényi β-expansion, respectively by frac(

x) = -1 i=-∞ x i β i its fractional part. The element 1 = β 0 belongs to Z + β .
Let us now turn to the case where β is a positive real algebraic integer. Then there ex-

ists an irreducible polynomial P(X ) = X m -m-1 i=0 a i X i , a i ∈ Z with m = degree(β) such that P(β) = 0. Then β = m-1 i=0 a m-1-i β -i . If a j
0 for all j and (a n , a n+1 , . . .) < (a m-1 , a m-2 , . . . , a 0 , 0, 0, . . .) for all n m -2, then the Rényi β-expansion of β would be

m-1 i=0 a m-1-i β -i from which we would deduce d β (1) = m-1 i=0 a m-1-i β -i-1
as well. But the coefficients a i do not obey the conditions of Parry in general. More considerations on the relations between β-expansions and algebraicity can be found in [Be] [Sch]. Bertrand-Matthis [Be] and Schmidt [Sch] have proved that, when β is a Pisot number, x ∈ Q(β) if and only if the Rényi β-expansion of x is eventually periodic; in particular the Rényi β-expansion of any Pisot number is eventually periodic.

[Be1] [Be2] [Be3] [Fro1] [Ak] [Ak1]
Let us recall that a Perron number β, resp. a Lind number, resp. a Salem number, will be a real algebraic integer β > 1 whose conjugates β (i) are of modulus strictly less than β, resp. of modulus less than β with at least one conjugate of modulus β [La], resp. of modulus less than 1 with at least one conjugate of modulus one. A Pisot number β will be a real algebraic integer β > 1 for which all the conjugates are in the open unit disc in the complex plane.

Canonical cut-and-project schemes over

Z β Assume that β > 1 is a Perron number of degree m 2, dominant root of the irreducible polynomial P(X ) = X m -a m-1 X m-1 -a m-2 X m-2 -• • • -a 1 X -a 0 , a i ∈ Z, a 0 ≠ 0.
All the elements rβ k with k 1, r ∈ {1, 2, . . . , β } are obviously in Z β . We are looking for asymptotic linear invariants associated with them, hence, by linearity, associated with the powers β k , k 1, of β, when k tends to infinity. By linearity, they will be also associated to the beta-integers. Let us set up the general situation. For all k 0 , write

β k = z m-1,k β m-1 + z m-2,k β m-2 + • • • + z 1,k β + z 0,k , where all the integers z 0,k , z 1,k , • • • , z m-1,k be- long to Z. Denote Z k = t z 0,k z 1,k z 2,k . . . z m-1,k , B = B (0) = t 1 β β 2 . . . β m-1 , B ( j) = t 1 β ( j) β ( j) 2 . . . β ( j) m-1
, where t means transposition and the elements

β ( j) , j ∈ {1, 2, • • • , m -1}, are the conjugate roots of β = β (0) in the minimal polynomial of β . Set B k =             β k β (1) k β (2) k . . . β (m-1) k             and Q =            0 1 0 • • • 0 0 0 1 • • • 0 . . . . . . . . . 0 0 0 • • • 1 a 0 a 1 • • • a m-1           
the m × m matrix with coefficients in Z. The transposed matrix of Q is denoted by t Q. It is the companion matrix of P(X ) (and of β). For all p, k ∈ {0, 1, • • • , m -1}, we have:

z p,k = δ p,k
the Kronecker symbol. It is obvious that, for all k 0, we have

Z k+1 = t Q Z k . Denote C =         1 β β 2 • • • β m-1 1 β (1) β (1) 2 • • • β (1) m-1 . . . . . . . . . . . . 1 β (m-1) β (m-1) 2 • • • β (m-1) m-1        
the Vandermonde matrix of order m. We obtain C Z k = B k by the real and complex embeddings of Q[β] since all the coefficients z j,k , j ∈ {0, 1, • • • , m -1}, are integers and remain invariant under the conjugation operation.

T 3.1. -If V 1 denotes the vector defined by the first column of C -1 , then the limit

lim k→+∞ Z k -1 Z k exists and is equal to the unit vector u := V 1 -1 V 1 .
Moreover, all the components of V 1 are real and belong to the Z -module

Z[β]
β m-1 P (β) .

Proof. -Since P(X ) is minimal, all the roots of P(X ) are distinct. Hence, the determinant

of C is i< j (β (i) -β ( j) ) and is not zero. Let C -1 = (ξ i j ). Then C • C -1 = I , that is ξ 1i +ξ 2i β ( j) +ξ 3i β ( j) 2 +• • •+ξ mi β ( j) m-1 = δ i, j+1 , i = 1, 2, . . . , m, j = 0, 1, . . . , m -1 (1)
On the other hand, the Lagrange interpolating polynomials associated with {β, β (1) , β (2) , . . . , β (m-1) } are given by

L s (X ) = m-1 j=0 j≠s X -β ( j) β (s) -β ( j) s = 0, 1, . . . , m -1.
For m arbitrary complex numbers

y 1 , y 2 , • • • , y m , let us denote by σ r = σ r (y 1 , y 2 , • • • , y m ) = 1 i 1 i 2 ••• ir m r j=1 y i j the r-th elementary symmetric function of the m numbers y 1 , y 2 , • • • , y m . The degree of L s (X ) is m -1 and L s (X ) can be expressed as L s (X ) = m-1 r=0 (-1) r σ (s) r X m-r-1 / m-1 r=0 r≠s (β (s) -β (r ) )
where σ (s)

r = σ r (β, β (1) , • • • , β (s-1) , β (s+1) , • • • , β (m-1)
) denotes the r-th elementary symmetric function of the m -1 numbers β, β (1)

, • • • , β (s-1) , β (s+1) , • • • , β (m-1) where β (s) is missing. Since these polynomials satisfy L s (β (k) ) = δ s,k for all s, k = 0, 1, • • • , m -1, comparing with
(1), we obtain, by identification of the coefficients

ξ j i = (-1) m-j σ (i-1) m-j m-1 r=0 r≠i-1 (β (i-1) -β (r ) ) = (-1) m-j σ (i-1) m-j P (β (i-1) ) for all i, j = 1, 2, • • • , m. We have: L s (X ) = m j=1 ξ j,s+1 X j-1 , s = 0, 1, • • • , m -1. Now C • Z k = B k for all k 0, hence Z k = C -1 • B k . Each component z i,k , 0 i m -1, k 0 of Z k can be expressed as z i,k = m j=1 ξ i+1, j β ( j-1) k . Since β is a Perron number, we have |β ( j) | < β for all j, 1 j m -1. Hence, for all j, 1 j m -1, lim k→+∞ β ( j) β k
= 0, and therefore lim k→+∞ z i,k

β k = ξ i+1,1 , i = 0, 1, • • • , m -1. Moreover, lim k→+∞    m-1 i=0 |z i,k | 2    1/2 β k = lim k→+∞ Z k β k = m-1 i=0 |ξ i+1,1 | 2 = V 1
hence the result. The fact that all the components of V 1 are real and belong to the Z -module Z[β]/(β m-1 P (β)) comes from the following more precise Proposition.

P 3.2. -The components (ξ j,1 ) j=1,...,m of V 1 are given by the following explicit functions of the coefficients a i of P(X ):

ξ j,1 = a j-1 β j-1 +a j-2 β j-2 +•••+a 1 β+a 0 β j P (β)
. In particular,

ξ m,1 = 1 P (β) . Proof. -We have L 0 (X ) = m j=1 ξ j,1 X j-1 and P(X ) = m-1 j=0 (X -β ( j) ) = L 0 (X )(X - β)P (β). All the coefficients of L 0 (X ) satisfy the following relations: -βP (β)ξ 1,1 = -a 0 , -βP (β)ξ 2,1 +ξ 1,1 P (β) = -a 1 , -βP (β)ξ 3,1 +ξ 2,1 P (β) = -a 2 , • • • , -βP (β)ξ m,1 +ξ m-1,1 P (β) = -a m-1 , ξ m,1 P (β) = 1. Hence the result recursively from ξ 1,1 noting that P (β) ∈ R -{0}. T 3.3. -Let u B := B/ B . Then: (i) u • u B = B -1 V 1 -1 > 0, (ii) the limit lim k→+∞ Z k+1 Z k exists and is equal to β , (iii) u is an eigenvector of t Q of eigenvalue β and the eigenspace of R m associated with the eigenvalue β of t Q is Ru, (iv) u B is an eigenvector of the adjoint matrix ( t Q) * = Q associated with the eigenvalue β and for all x ∈ C m : lim k→+∞ β -k ( t Q) k (x) = (x • B) V 1 . Proof. -(i) and (ii): From the relation C • C -1 = I d we deduce the equality V 1 • B = 1. Hence u • B = V 1 -1 > 0 . Then, for all k 0, t Z k • B = β k = Z k t ( Z k Z k -u + u) • B > 0 which tends to infinity when k tends to + ∞. Since u -Z k / Z k tends to zero when k goes to infinity, Z k behaves at infinity like β k / (u • B), hence the limit; (iii): for all k 0, t Q(u) = t Q(u -Z k Z k + Z k Z k ) = t Q(u -Z k Z k ) + Z k+1 Z k Z k+1 Z k+1
. The first term is converging to zero and the second one to βu when k goes to infinity, from Theorem 3.1. Hence, the result since all the roots of P(X ) are distinct and the (real) eigenspace associated with -1) . Let us specify their respective actions on R m .

β is 1 -dimensional; (iv): it is clear that B is an eigenvector of the adjoint matrix Q . If h 0 , h 1 , • • • , h m-1 ∈ C , x = m-1 j=0 h j Z j , where Z 0 , Z 1 , • • • , Z m-1 is the canonical basis of C m , we have: β -k ( t Q) k (x) = m-1 j=0 h j β -k Z k+ j = m-1 j=0 h j β j Z k+ j β k+ j , but, from the proof of Theorem 3.1, lim k→+∞ Z k+ j β k+ j = V 1 and m-1 j=0 h j β j = x • B . We deduce the claim. Let us denote by t Q C the automorphism of C m which is the complexification operator of t Q . Its adjoint Q C obviously admits {B, B (1) , B (2) , • • • , B (m-1) } as a basis of eigenvectors of respective eigenvalues β, β (1) , β (2) , • • • , β (m
Let s 1 , resp. t , be the number of real, resp. complex (up to conjugation), embeddings of the number field Q(β) . We have m = s + 2t . Assume that the conjugates of β are

β, β (1) , • • • , β (s-1) , β (s) , β (s+1) , • • • , β (m-2) = β (s+2t -2) , β (m-1) = β (s+2t -1) where β (q) is real if q s -1 and β (s+2 j) = β (s+2 j+1) = |β (s+2 j) |e iθ j , j = 0, 1, • • • , t -1, is complex with non- zero imaginary part. Let us recall that V 1 denotes the vector defined by the first column of C -1 (Theorem 3.1). C 3.4. -(i) A basis of eigenvectors of t Q C is given by the m column vectors {W k } k=1,2,••• ,m of respective components ξ j,k = a j-1 β (k-1) j-1 +a j-2 β (k-1) j-2 +•••+a 1 β (k-1) +a 0 β (k-1) j P (β (k-1) ) with j = 1, 2, • • • , m; in particular, ξ m,k = 1 P (β (k-1) ) ; (ii) a real Jordan form for t Q is given by the diagonal matrix Diag(β, β (1) , • • • , β (s-1) , D 0 , D 1 , • • • , D t -1 ) in the basis of eigenvectors {V j } j=1,••• ,m with V 2 = W 2 , • • • , V s = W s , V s+2 j+1 = Im(W s+2 j+1 ), V s+2 j+2 = Re(W s+2 j+1 ), j = 0, 1, • • • , t -1 and where the 2 × 2 real Jordan blocks D j are   |β (s+2 j) | cos θ j -|β (s+2 j) | sin θ j |β (s+2 j) | sin θ j |β (s+2 j) | cos θ j   ;
(iii) a real Jordan form of the adjoint operator ( t Q) * = Q is given by the same diagonal matrix

Diag(β, β (1) , • • • , β (s-1) , D 0 , D 1 , • • • , D t -1 ) in the basis of eigenvectors {X j } j=1,••• ,m with X 1 = B, X 2 = B (1) , X 3 = B (2) , • • • , X s = B (s-1) , X s+2 j+1 = Im(B (s+2 j) ), X s+2 j+2 = Re(B (s+2 j) ), j = 0, 1, • • • , t -1. The t planes RX s+2 j+1 + RX s+2 j +2 , j = 0, 1, . . . , t -1 are all orthogonal to V 1 ,
and thus also to u.

Proof. -(i): We apply, componentwise in the equation

( t Q)V 1 = βV 1 , the Q -automor-
phisms of C which are the real and complex embeddings of the number field Q(β). Since t Q has rational entries and V 1 has its components in the Z-module

β 1-m (P (β)) -1 Z[β],
we deduce the claim:

( t Q)W j = β ( j-1) W j with j = 1, 2, • • • , m and where W 1 = V 1 ; (ii): the restric- tions of t Q C to the (real) t Q -invariant subspaces of R m
have no nilpotent parts since all the roots of P(X ) are distinct. Hence, a real Jordan form of t Q is the one proposed with Jordan blocks which are 2×2 on the diagonal [HS]. (iii): in a similar way the equation QB = βB implies

QB ( j) = β ( j) B ( j) with j = 0, 1, • • • , m -1.
Obviously Q C and t Q C have the same eigenvalues and Q and t Q the same real 2 × 2 Jordan blocks on the diagonal. The corresponding basis of eigenvectors is given by the vectors X i [HS]. The orthogonality between V 1 and the vector X s+2 j+1 , resp. X s+2 j+2 , j = 0, 1, . . . , t -1, arises from the relation C • C -1 = I d . We deduce the claim for the planes.

The linear invariants associated with the powers of β are the invariant subspaces given by Corollary 3.4. Let us turn to the beta-integers. Beta-integers are particular Z-linear combinations of powers of β. We will show how to construct the set Z β using the above linear invariants, namely, the set Z β will appear in a natural way on the line Ru B as image of a point set close to the expanding line Ru.

R . -The conditions of Parry, used here in the context of matrices t Q without any condition on the signs of the entries, give the same results as those obtained with the Perron-Frobenius theory (Minc [Mi]), when this one is applicable, that is when t Q has non-negative entries: first, the dimensionality one for the dominant eigenspace of t Q ; second, the equality

lim k→+∞ β -k ( t Q) k (x) = (x • B) V 1 , for x ∈ C m , in Theorem 3.3 (compare with Ruelle [Ru]
p136 when t Q has non-negative entries), and its consequences.

T 3.5. -Let π B be the orthogonal projection mapping of R m onto RB and define

L = {x k Z k + x k-1 Z k-1 + • • • + x 1 Z 1 + x 0 Z 0 | x i ∈ A, k 0 , and (x j , x j-1 , • • • , x 1 , x 0 , 0, 0, • • • ) < (c 1 , c 2 , • • • ) for all j, 0 j k } the t Q-invariant subset of Z m . Then: (i) the mapping k j=0 x j β j → k j=0 x j Z j : Z + β → L (with the same coefficients x j ) is a bijection, (ii) the mapping π B | Z m is one-to-one onto its image Z[β] B -1 u B : for any k 0, a 0 , a 1 , • • • , a k ∈ Z, we have π B k i=0 a i Z i = k i=0 a i β i B -1 u B and
(L ) = Z + β B -1 u B .
Proof. -(i): this mapping Z + β → L is obviously surjective. Let us show that it is injective.

Assume there exists a non-zero element

k j=0 x j β j in Z + β such that k j=0 x j Z j = 0 . Since t k j=0 x j Z j B = 0 = k j=0
x j β j , this would mean that zero could be represented by a non-zero element. This is impossible by construction; (ii): for all k 0 , we have

π B (Z k ) = β k B -1 u B , hence the result by linearity. The injectivity of π B | Z m comes from the assertion (i). P 3.6. -Let u B,i = X i -1 X i if i = 1, 2, • • • , s, u B,i = ( X i 2 + X i+1 2 ) 1/2 ( Re( B (i-1) -2 )X i + I m( B (i-1) -2 )X i+1 ) if i = s + 1, • • • , m with i -(s + 1) even, and u B,i = ( X i-1 2 + X i 2 ) 1/2 (-I m( B (i-1) -2 )X i +Re( B (i-1) -2 )X i+1 ) if i = s +1, • • • , m with i -(s+1) odd. Denote by π B,i : R m → Ru B,i , i = 1, 2, • • • , s the orthogonal projection mappings to the 1-dimensional eigenspaces of Q , resp. π B,i : R m → Ru B,i +Ru B,i+1 , i = s +1, • • • , m with i -(s + 1)
even, the orthogonal projection mappings to the irreducible 2-dimensional eigenspaces of Q .

Then, for all k

0, a 0 , • • • , a k ∈ Z, we have π B,i ( k j=0 a j Z j ) = k j=0 a j β (i-1) j X i -1 u B,i , i = 1, 2, • • • , s and, for all i = s + 1, • • • , m with i -(s + 1) even, π B,i ( k j=0 a j Z j ) = 1 m-1 k=0 |β (i-1) | 2k 1/2   Re( k j=0 a j β (i-1) j ) I m( k j=0 a j β (i-1) j ) -I m( k j=0 a j β (i-1) j ) Re( k j=0 a j β (i-1) j )     u B,i u B,i+1   .
Proof. -It suffices to apply the real and complex embeddings of Q(β) to the relation

π B ( k j=0 a j Z j ) = (( k j=0 a j Z j ) • B) B -2 B = k j=0 a j β j B -2 B: for complex embed- dings, X i 2 + X i+1 2 = m-1 k=0 |β (i-1) | 2k and B (i-1) -2 B (i-1) means:   Re( B (i-1) -2 ) I m( B (i-1) -2 ) -I m( B (i-1) -2 ) Re( B (i-1) -2 )     X i X i+1   = 1 ( X i 2 + X i+1 2 ) 1/2   u B,i u B,i+1   .
The explicit expressions given above will allow us below to compare the "geometric" Rauzy fractals deduced from the present study and the "algebraic" Rauzy fractal. Before stating the main theorem about the existence of canonical cut-and-projection schemes associated with the beta-integers when β is a general (non-integer) Perron number, let us first consider the case of equality u = u B and show that it is rarely occuring.

P 3.7. -The equality u = u B holds if and only if β is a Pisot number, root > 1 of the polynomial X 2 -aX -1 , with a 1 . Proof. -The condition u = u B is equivalent to V 1 colinear to B , that is ξ j,1 β -j+1 = a non-zero constant, for all j = 1, 2, • • • , m. The condition is sufficient: if β is such a Pisot
number, such equalities hold. Conversely, if such equalities hold, this implies in particular that

ξ 1,1 β -1+1 = ξ m,1 β -m+1 . Thus we obtain a 0 β m-2 = 1 , that is necessarily m = 2 and a 0 = 1 .
The Perron number β is then a Pisot number of negative conjugate -β -1 which satisfies β 2a 1 β -1 = 0 , where a 1 = ββ -1 is an integer greater than or equal to 1. This is the only possibility of quadratic Pisot number of norm -1 ([Fro1], Lemma 3).

T 3.8. -Denote by E the line Ru B in R m . There exist two canonical cut-and-

project schemes E p 1 ←-(E × D R m , Z m ) p 2 -→ D associated with Z β ⊂ E (see the definitions
in the Appendix). They are given by, in case (i): the orthogonal projection mapping π B as p 1 , ⊕ F F as internal space D, p 2 = ⊕ F π F , where the sums are over all irreducible t Q-invariant subspaces F of R m except Ru and where π F is the projection mapping to F along its t Qinvariant complementary space, in case (ii): as p 1 the orthogonal projection mapping π B , ⊕ F F as internal space D where the sum is over all irreducible Q-invariant subspaces F of R m except E , as p 2 the sum ⊕ i≠1 π B,i of all the orthogonal projection mappings except π B,1 = π B ; in the case (ii), the internal space D is orthogonal to the line Ru.

Proof. -In both cases, the fact that p 2 (Z m ) is dense in D arises from Kronecker's theorem (Appendix B in [Mey]): since β is an algebraic integer of degree m, the m real num-

bers 1 = β 0 , β 1 , • • • , β m-1
are linearly independent over Q . Hence, for all > 0 and all m-tuple of real numbers x 0 , x 1 , • • • , x m-1 such that the vector (say) x = t (x 0 x 1 . . . x m-1 ) belongs to D , there exist a real number w and m rational integers

u 0 , u 1 , • • • , u m-1 such that | x j -β j w -u j | / √
m for all j = 0, 1, . . . , m -1. In other terms, there exists a point u = t (u 0 , u 1 , • • • , u m-1 ) ∈ Z m such that its image p 1 (u) is wB ∈ Ru B and its image p 2 (u) is close to x up to . Hence the result. As for the restriction of the projection mapping

p 1 = π B = π B,1 : R m → E to the lattice Z m ,
it is injective after Theorem 3.5. The orthogonality between D and u comes directly from Corollary 3.4.

The mapping p

1 (Z m ) → D : x → x * = p 2 • ( p 1 | Z m ) -1 (x)
will be denoted by the same symbol (.) * in the cases (i) and (ii), the context making the difference. P 3.9. -Let β be a Pisot number, root > 1 of the polynomial

X 2 -aX -1 , with a 1. Put c a = (1+aβ) β √ 2+aβ(β-1)
. Then the two canonical cut-and-project schemes given by (i) and (ii) in Theorem 3.8 are identical and the inclusion of Z β B -1 u B in the following model set holds:

Z β B -1 u B = π B (L ∪ (-L )) ⊂ {v ∈ π B (Z 2 ) | v * ∈ [-c a u B,2 , +c a u B,2 ]} where u B,2 = t (-β 1) B -1 .
Proof. -The two cut-and-project schemes are identical: by Proposition 3.7 the equality u = u B holds and the line Ru B,2 , which is obviously orthogonal to the line Ru B , is t Qinvariant. Now, if g denotes an arbitrary element of L , it can be written g

= x k ( t Q) k Z 0 + x k-1 ( t Q) k-1 Z 0 + • • • + x 1 ( t Q)Z 0 + x 0 Z 0 for a certain integer k 0 with x i ∈ A and (x j , x j-1 , • • • , x 1 , x 0 , 0, 0, • • • ) < (c 1 , c 2 , • • • ) for all j, 0 j k . We have Z 0 = su + s ⊥ u B,2 with s = B -1 and s ⊥ = -β B -1 . Then g = k j=0 x j ( t Q) j Z 0 = k j=0 x j sβ j u + s ⊥ (-1) j β -j u B,2 . Thus p 2 (g) = s ⊥ k j=0 x j (-1) j β -j u B,2 and p 1 (g) * = p 2 (g) |s ⊥ | β +∞ j=0 β -j = |s ⊥ | β 1 1-β -1 which is equal to c a since B = √ 2 + aβ .
This constant is independent of k . Hence we have

p 1 (g) ⊂ {v ∈ π B (Z 2 ) | v * ∈ [-c a u B,2 , +c a u B,2 ]} and the claim. Let C = { m-1 j=0 α j Z j | α j ∈ [0; 1] for all j = 0, 1, • • • , m -1 } be the m-cube at the origin. For all irreducible t Q-invariant subspace F of R m , put δ F = max x∈C π F (x) , λ F the absolute value of the eigenvalue of t Q on F and c F = β δ F 1-λ m F . Denote by Ω F the closed interval centred at 0 in F of length 2c F if dim F = 1, resp. the closed disc centred at 0 in F of radius c F if dim F = 2.
T 3.10. -Let β be a Pisot number of degree m 2 and Ω = ⊕ F Ω F where the sum is over all irreducible t Q-invariant subspace F of R m except Ru. Then the inclusion of

Z β B -1 u B in the following model set defined by Ω holds: Z β B -1 u B = p 1 (L ∪ (-L )) ⊂ {v ∈ p 1 (Z m ) | v * ∈ Ω }
in the cut-and-project scheme given by the case (i) in Theorem 3.8.

Proof.

-

If g = k j=0 x j Z j ∈ L with k = d m -1 , and d 1 an integer, then g = d-1 q=0 m-1 l =0 x qm+l ( t Q) qm Z l = d-1 q=0 ( t Q) qm    m-1 l =0 x qm+l Z l    .
Hence p 1 (g

) * = p 2 (g) = ⊕ F π F (g) = F d-1 q=0   ( t Q | F ) qm π F    m-1 l =0 x qm+l Z l       , with: π F (g) d-1 q=0 β λ qm F π F    m-1 l =0 Z l    β δ F +∞ q=0 λ qm F = β δ F 1 -λ m F = c F .
This constant is independent of d, hence of k = d m -1. It is easy to check that it is an upper bound for p 2 (g) if k ≡ -1(mod m) and also for all g ∈ -L . We deduce the claim.

C 3.11. -If β is a Pisot number of degree m 2 , then Z β is a Meyer set.

Proof. -If β a Pisot number, the set Z β , viewed as the set of vertices of an aperiodic tiling, is obtained by concatenation of prototiles on the line, which are in finite number by Thurston [Th]. And it is relatively dense by construction. Now, by Theorem 3.10 it is included in a model set. This proves the claim (see the Appendix).

In both cases of cut-and-project scheme, as given by Theorem 3.8 where the duality between the matrices Q and t Q clearly appears, the internal space represents the contracting hyperplane, whereas the line Ru is the expanding direction, when β is a Pisot number. The duality between both cut-and-project schemes is connected to the substitutive approach by the following (Arnoux and Ito [AI], Chap. 7 in Pytheas Fogg [PF]): the abelianized Z k of the iterates of the substitution satisfy Z k+1 = QZ k , and gather now about the line RB. If one takes the projection on RB of the new set L (defined similarly as L ) along the other eigenspaces, one recovers Z β (up to a scalar factor). A striking feature of the internal spaces is that the numeration in base β ( j) (conjugates of β) appears as canonical ingredient to control the distance between a point of L and its orthogonal projection to the expanding line Ru, in particular at infinity. D 3.12. -Let β be a Pisot number of degree m 2 . The closure

Z + β B -1 u B * of
the set p 2 (L ) is called the canonical acceptance window associated with the set of beta-integers

Z +
β in both cases (case (i) or (ii) in Theorem 3.8) of cut-and-project scheme: in the case (i) it will be denoted by R i and in the case (ii) by R .

The notations R and R i (R i ⊂ Ω by Theorem 3.10) with an " R " like Rauzy are used to recall the close similarity between these sets and the Rauzy fractal (Rauzy [Ra], Arnoux and Ito [AI], Messaoudi [Me], Ito and Sano [IS], Chap. 7 in Pytheas Fogg [PF]). The fact is that the set R is exactly the Rauzy fractal up to the multiplication by a non-zero scalar factor on each irreducible Q-invariant subspace (by definition we will speak of quasi-dilation). Let us show it on an example.

"Tribonacci" case [Me]: let us consider the irreducible polynomial P(X ) = X 3 -X 2 -X -1.

Its dominant root is denoted by β, and α and α are the two other complex conjugates roots of P(X ). In this case, the Rauzy fractal is "algebraically" defined by E := { ∞ i=3 i α i | i ∈ {0, 1} and i i+1 i+2 = 0 for all integer i 3}. The condition imposed on the sequence ( i ) i 3 is exactly that given by the conditions of Parry. Indeed [START_REF] F | Finite beta-expansions[END_REF] and section 2), d β (1) = 0.111 and the lexicographical maximal sequence is

c 1 c 2 c 3 • • • = (110) ω . Now (Proposition 3.6) B (1) = t (1 α α 2 ) and X 2 2 + X 3 2 = 1 + αα + α 2 α 2 = β. We deduce that R = β -1/2 E = π B,2 (L ) with the following (metric) identification of C : φ(C) = Ru B,2 + Ru B,3
where φ is the isometry which sends the vector

  1 0   , resp.   0 1   , to   u B,2 u B,3   , resp. to   u B,3 -u B,2   .
P 3.13. -The canonical acceptance window R (relative to the case (ii) of cutand-project scheme in Theorem 3.8) is compact and connected. Its interior int (R ) is simply connected, contains the origin. The set R is such that:

(i) int (R ) = R ; (ii) it induces a
tiling of the internal space D modulo the lattice φ(Z + Zα):

D = z∈φ(Z+Zα) (R + z); (iii) (R + z) ∩ int (R + z ) = ∅ for all z, z ∈ φ(Z + Zα), z ≠ z .
Proof. -Since R = β -1/2 E , we deduce the properties of R from those of E already established in Rauzy [Ra], Messaoudi [Me] and [START_REF] M | Frontière du fractal de Rauzy et système de numérotation complexe[END_REF].

P 3.14. -The boundary of R is a fractal Jordan curve. A point z belongs to the boundary of R if and only if it admits at least 2 distinct Rényi α-expansions. A point belonging to the boundary of R admits 2 or 3 distinct Rényi α-expansions, never more.

Proof. -The properties of the boundary of E are given in Ito and Kimura [IK] and Messaoudi [START_REF] M | Frontière du fractal de Rauzy et système de numérotation complexe[END_REF]). Hence the claim.

The properties of R i follow from the equality: p 2 (R ) = R i , where p 2 refers to the case (i) of cut-and-project scheme in Theorem 3.8, and from Proposition 3.13 and 3.14: in particular, it has also a fractal boundary. We will speak of "geometrical" Rauzy fractals for R and R i and of "algebraic" Rauzy fractal for E . They are similar objects as far as they concentrate all the information about the beta-integers and the completions of their real and complex embeddings (Rauzy [Ra]). The respective canonical acceptance windows associated with Z β are R ∪ (-R ) and R i ∪ (-R i ) in the two cut-and-project schemes.

Additive properties of Z β

In this section, β will be a Pisot number of degree m 2.

A) Cones, generators and semi-groups.-We will show that any element of L is generated by a finite number of elements of L of small norm, over N. By projection to E by π B , the ambiant 1-dimensional space of the beta-integers (Theorem 3.8), this will imply the same property for Z β . This finiteness property, stated in Corollary 4.6, constitutes a refinement of Theorem 4.16 (i) (Lagarias) for the Meyer sets Z β .

First let us fix the notations and simplify them somehow. Let π : R m → Ru be the projection mapping along its t Q-invariant complementary space (instead of denoting it by π Ru ), and p 2 the projection mapping of the cut-and-project scheme (i) in Theorem 3.8. Let π : R m → Ru be the orthogonal projection mapping and π ⊥ = I dπ (π ⊥ is the mapping ⊕ i≠1 π b,i in the case (ii) of cut-and-project scheme in Theorem 3.8). The basic ingredient will be the construction of semi-groups of finite type associated with cones whose axis of revolution is the expanding line Ru, following an idea of Lind [Li] in another context.

Truncating them in a suitable way at a certain distance of the origin will be the key for finding generators of L over N. In the first Lemma we will consider the possible angular openings of these cones around the expanding line Ru for catching the points of L . For θ > 0 , define the cone

K θ := {x ∈ R m | θ p 2 (x) π(x) , 0 π(x) • u }. For r, w > 0 , define K θ (r) := {x ∈ K θ | π(x) r }, K θ (r, w) := {x ∈ K θ | r π(x) w }. If A is an arbitrary subset of R m , denote by sg(A ) := { finite m i x i | m i ∈ N, x i ∈ A } the
semigroup generated by A . Let ρ be the covering radius of the subset L ∪ (-L ) with respect to the band R i × Ru: ρ is the smallest positive real number such that for any z ∈ R m such that p 2 (z) ∈ R i the closed ball B(z, ρ) contains at least one element of L ∪ (-L ). A lower bound of ρ is given by the covering radius √ m/2 of the lattice Z m . Referring R m to the basis {B, V 2 , V 3 , • • • , V s+2t } and using Corollary 3.4 and Theorem 3.10 we easily deduce the following upper bound of ρ: 1 2 B -1 + F c F , where the sum F means there and in the following everywhere it will be used "the sum over all irreducible t Q-invariant subspaces F of R m except Ru". The notation diam(•) will be put for the diameter of the set (•) in the following. P 4.1. -(i) For all θ > 0 , there exists an integer j 0 = j 0 (θ) 0 such that

Z j ∈ K θ for all j j 0 ; (ii) if t Q is nonnegative, and min{ξ j,1 | j = 1, 2, • • • , m } > 2 V 1 (diam(R i ))
, then the following equality j 0 (θ) = 0 holds for all 0 < θ < θ min , where

θ min := -2 + (diam(R i )) -1 V 1 -1 min{ξ j,1 | j = 1, 2, • • • , m }. Proof. -(i) Let θ > 0 .
We have just to prove that π(Z j ) • u tends to +∞ and not to -∞ when j goes to +∞ . Let j 0.

Write Z j = π(Z j ) + p 2 (Z j ) = π (Z j ) + π ⊥ (Z j ); hence π(Z j ) -π (Z j ) = π ⊥ (Z j ) -p 2 (Z j ) π ⊥ (Z j ) + p 2 (Z j ) = π ⊥ (Z j -π(Z j )) + p 2 (Z j ) 2 p 2 (Z j ) 2diam(R i ). On the other hand π B (π (Z j )) -π B (Z j ) = (Z j • u)(u B •u)u B -B -1 β j u B = | (Z j •u)(u B •u)-B -1 β j | = π B (π ⊥ (Z j )) π ⊥ (Z j ) diam(R ). Hence, since u B • u > 0 (Theorem 3.3), | Z j • u -(u B • u) -1 B -1 β j | (u B • u) -1 diam(R ). Consequently π(Z j ) -(u B •u) -1 B -1 β j u = π(Z j ) -π (Z j ) + π (Z j ) - (u B •u) -1 B -1 β j u π(Z j )-π (Z j ) + π (Z j )-(u B •u) -1 B -1 β j u 2diam(R i )+ (u B • u) -1 )diam(R ). The quantity π(Z j ) • u tends to +∞ as (u B • u) -1 B -1 β j when j → +∞ . Then there exists j 0 such that Z j • u 2diam(R i ) + ((u B • u) -1 + θ)diam(R i ), for all j j 0 . As a consequence π(Z j ) • u Z j • u -2diam(R i ) ((u B • u) -1 + θ)diam(R i ) > 0 for all j j 0 . We claim that Z j ∈ K θ for all j j 0 . Indeed, since p 2 (Z j ) diam(R i ) , the inequalities hold: θ p 2 (Z j ) (θ + (u B • u) -1 ) p 2 (Z j ) (θ + (u B • u) -1 )diam(R i ) π(Z j ) • u = π(Z j ) for all j j 0 .
(ii) If t Q is nonnegative the coefficients a i in P(X ) are nonnegative with a 0 ≠ 0 and at least one of the coefficients a k , k 1, is non-zero since β is assumed to be a Pisot number and not a Salem number. Hence (Proposition 3.2), since

P (β) > 0, we have V 1 -1 ξ 1,1 = π (Z 0 ) = Z 0 • u = a 0 βP (β) > 0 and V 1 -1 ξ j+1,1 = π (Z j ) = Z j • u, for all j = 1, 2, • • • m - 1 with V 1 -1 min{ξ j,1 | j = 1, 2, • • • m -1} a 0 β m P (β) > 0. Because {Z 0 , Z 1 , • • • , Z m-1 } is the canonical basis of R m , any Z j , j m
, can be written as a combination of the elements of this basis with positive coefficients. Hence,

Z j • u V 1 -1 min{ ξ l ,1 | l = 1, 2, • • • , m -1 } for all j 0. But the relation Z j = π(Z j )+ p 2 (Z j ) = π (Z j )+π ⊥ (Z j ) implies that π(Z j )-π (Z j ) = π ⊥ (Z j ) -p 2 (Z j ). Hence, |π(Z j ) • u -Z j • u| π ⊥ (Z j ) + p 2 (Z j ) 2 p 2 (Z j ) 2(diam(R i )) for all j 0. Therefore π(Z j ) • u V 1 -1 min{ ξ l ,1 | l = 1, 2, • • • , m -1 } - 2diam(R i )
which is > 0 by assumption for all j 0. Hence, by definition of

θ min , π(Z j )• u = π(Z j ) θ min (diam(R i )) θ min p 2 (Z j ) θ p 2 (Z j
) for all j 0 and 0 < θ θ min .

We deduce that Z j ∈ K θ for all j 0 and 0 < θ θ min . Let us observe that the conditions of the present assertion are generally not fulfilled.

We now turn to the question of generating the elements of L by a finite number of them over N. The idea we will follow is simple: let us consider the set of the semi-groups generated by a finite number of (arbitrary) elements of L ∩ K θ for all θ > 0; in this set, we will show the existence of semi-groups (θ > 0 fixed) containing K 2θ ∩ L , that is containing L except a finite number of elements of L close to the origin. Then we will minimize this finite number of excluded elements. For this we will consider the maximal possible values of θ. In final this will provide a suitable value of θ and a control of the norms of the generating elements of the semi-group which will contain L .

L 4.2. -(Lind [Li])

Let θ > 0 . If δ = (2θ + 2) -1 and x ∈ K 2θ with π(x) = π(x) • u > 4, then [x -K θ (1, 3)] ∩ K 2θ contains a ball of radius δ.
Proof. - [Li] Take y = 2u + 3(π(x) • u) -1 p 2 (x) . We will show that the ball centred at xy and of radius δ satisfies our claim. Suppose z < δ . Then xy + z ∈ K 2θ . Indeed,

2θ p 2 (x -y + z) 2θ 1 -3(π(x) • u) -1 p 2 (x) + δ 1 -3(π(x) • u) -1 (π(x) • u) + 2θ(2θ + 2) -1 = (π(x) • u) -2 -2(2θ + 2) -1 but p 2 (y) = 2 . We deduce 2θ p 2 (x -y + z) π(x -y + z) • u. Let us show that y -z ∈ K θ . We have 2θ p 2 (y) = 6θ(π(x) • u) -1 p 2 (x) 3(π(x) • u) -1 (π(x) • u) = 3. Therefore θ p 2 (y -z) θ( p 2 (y) + δ) 3 2 + θ(2θ + 2) -1 = 2 -(2θ + 2) -1 π(y -z) • u. Now, since δ < 1 , we have the inequalities 1 π(y -z) • u 3, establishing the result. T 4.3. -Let θ > 0 . If r is such that r > ρ(2θ + 2), then K 2θ ∩ L ⊂ sg(K θ (r) ∩ L ).
Proof. -Lemma 4.2 implies the following assertion:

if x ∈ K 2θ is such that π(x) • u > 4r with r > ρ(2θ + 2) , then [x -K θ (r, 3r)] ∩ K 2θ contains a ball of radius rδ > ρ . But
ρ is by definition the covering radius of L ∪ (-L ), hence this ball intersects L . Now, let A = K θ (4r) ∩ L be the finite point set of L and let us show that K 2θ ∩ L ⊂ sg(A ). First the inclusion K 2θ (4r) ∩ L ⊂ sg(A ) holds. We now proceed inductively. Suppose K 2θ (r ) ∩ L ⊂ sg(A ) for some r 4r . We will show that this implies K 2θ (r + r) ∩ L ⊂ sg(A ) , which will suffice by induction. For this, let us take g ∈ L ∩ [K 2θ (r + r) K 2θ (r)]. From Lemma 4.2 and the above, there exists an element, say y , in L , contained in [g -K θ (r, 3r)] ∩ K 2θ (r ). By assumption, y ∈ sg(A ) and y = gx for some x ∈ K θ (r, 3r) ∩ L L ⊂ sg(A ). Therefore g = x + y ∈ sg(A ) + sg(A ) ⊂ sg(A ). This concludes the induction. L 4.4. -For all θ > 0 , the following set: -Any element g ∈ L \ L (θ f ) can be expressed as a finite combination over N of elements of the finite point set

L (θ) := { x ∈ L | p 2 (x) ∈ R i , x ∈ K θ (ρ(2θ + 2)), x ∈ K 2θ } is finite. Proof. -The proof is clear since all g ∈ L such that π(g) • u > 2ρ(2θ + 2) belongs to K 2θ . Define θ f := max{θ > 0 | #(L (θ)) is minimal } (where #(•) denotes the number of elements of the set (•)). If t Q is nonnegative
K θ f (ρ(2θ f + 2)) ∩ L . Proof. -It is a consequence of Theorem 4.3 with θ = θ f and r = ρ(2θ f + 2).
C 4.6. -There exist two disjoint finite subsets

F = B π B (g) • u B | g ∈ L (θ f )) and F = {g 1 , g 2 , • • • , g η } ⊂ B π B (g) • u B | g ∈ K θ f (ρ(2θ f + 2)) ∩ L of Z + β such that Z + β ⊂ F ∪ N[g 1 , g 2 , • • • , g η ].
(1)

The generating elements g i ∈ F satisfy: 1) and F is empty, then the degree m of β divides η.

g i ρ(2θ f + 2) B -1 V 1 -1 + diam(R i ). If the couple (F , F ) is such that η = #F is minimal for the inclusion relation (
Proof. -To obtain the inclusion (1) it suffices to project L by π B and to apply Theorem 4.3 and 4.5 and Lemma 4.4. Let us show the upper bound on the norms of the elements of F .

If g ∈ K θ f (ρ(2θ f + 2)) ∩ L is decomposed as g = π(g) + t , where t ∈ R i , then, by Theorem 3.3 (i), we have: π B (g) π(g) B -1 V 1 -1 + diam(R i ). But π(g) ρ(2θ f + 2).
We deduce the claim. Now if

Z + β ⊂ N[g 1 , g 2 , • • • , g η ] the group Z[g 1 , g 2 , • • • , g η ] contains Z β and the equality Z[Z β ] = Z[g 1 , g 2 , • • • , g η ] necessarily holds. By Theorem 4.16 we deduce that m divides η since the rank of Z[Z β ] = Z[N[Z β ]
] is η when η is the smallest integer such that the set inclusion (1) holds and that F is empty. B) Preperiods in the addition of beta-integers.-The Delone set Z β endowed with the usual addition and multiplication cannot have a ring structure otherwise it would contain Z but it is obvious that Z β contains no subset of the type λZ, λ > 0. This absence of ring structure on Z β for the usual laws can be partially overcome by controlling the fractional parts of the Rényi β-expansions of x + y and xy when x, y ∈ Z β . This is the aim of this paragraph to focus on the geometrical meaning of the sets T and T as stated in Theorem 4.8 and of the exponent q in its Corollary 4.9.

The projection mappings will be the ones redefined (in a simpler way) at the beginning of the subsection A). Let R > 0 and I be an interval of R having compact closure. Let us extend the m-cube C for reasons which will appear below. Let

C = { m-1 j=0 α j Z j | α j ∈ [-1; 1] for all j = 0, 1, • • • , m -1 }. For all irreducible t Q-invariant subspace F of R m , put δ F = max x∈C π F (x) , λ F the absolute value of the eigenvalue of t Q on F and c F = β δ F 1-λ m F . Denote by Ω F the closed interval centred at 0 in F of length 2c F if dim F = 1, resp. the closed disc centred at 0 in F of radius c F if dim F = 2. Let Ω = ⊕ F Ω F . We will denote by T I ,R := {x ∈ R m | p 2 (x) ∈ β -1 R Ω , π B (x) • u B ∈ B -1 I
} the slice of the band defined by β -1 R Ω in the internal space, extended by symmetrization with respect to β -1 R Ω (compare the definitions of Ω and Ω in Theorem 3.10), of axis the expanding line

Ru . Let F R := {frac(z) | z = a k β k + a k-1 β k-1 + • • • + a 1 β + a 0 , a i ∈ Z, |a i | R} ⊂ [0, 1). L 4.7. -The set { B π B (g) • u B | g ∈ T [0,1),R+ β ∩ Z m } is a finite subset of Z[β] ∩ [0, 1) and the following inclusion holds: F R ⊂ { B π B (g) • u B | g ∈ T [0,1),R+ β ∩ Z m }.
Proof. -The finiteness of the set is obvious: it is a discrete set in a subset of R m having compact closure. The inclusion relation is a reformulation of Lemma 2.1 in [Bu]. Let us briefly

recall the proof. Let z = k j=0 a j β j with a i ∈ Z, |a i | R . We have also z = k j=-∞ x j β j as β -expansion of z . Therefore z-int(z) = k i=0 a j β k -k j=0 x j β j . Since 0 x j β and |a i | R , frac(z) ∈ [0, 1
) is a polynomial in β , the coefficients of which have their absolute values bounded by R + β . Here the coefficients may be negative or positive. This is why we have introduced C instead of C . We deduce the result in a similar way as in the proof of Theorem 3.10 for the computation of the upper bound c F , except that now it is with Ω , c F and the fact that the absolute value of the digits is less than R + β ; this obliges to multiply Ω by the factor (R + β )/ β . The set F R is finite (Lemma 6.6 in [So]), and (Proposition 3.5) is in one-to-one correspondence with a subset of the finite point set T [0,1),R+ β . We deduce the claim.

Let L I ,R := min [ln(β (i-1) -1 )] -1 ln m-1 k=0 (β (i-1) ) 2k 1/2 ψ I ,R+ β
where the minimum is taken over the real positive embeddings of Q(β) ( i = 1, 2, . . . , s and β (i-1) > 0)

and where ψ I ,R := max{ y | y ∈ T I ,R }. Let us consider an element z ∈ F R . Its β -expansion:

+∞ j=1 z -j β -j is eventually periodic [Be] [Sch] and therefore can be written

k 0 (z) j=1 z -j β -j + +∞ k=0 k 0 (z)+(k+1)r (z)
j=k 0 (z)+k r(z)+1 z -j β -j where the integers k 0 (z), r(z) 1 are minimal. We will denote by J R = max{k 0 (z) | z ∈ F R } the maximal preperiod of the β -expansions of the elements of F R . An upper bound of J R will be computed below.

Z + β + Z + β ⊂ Z + β + T, Z + β -Z + β ⊂ Z β + T , where T = { B π B (g) • u B | g ∈ T [0,+1),3 β ∩ Z m } and T = { B π B (g) • u B | g ∈ T (-1,+1),2 β ∩ Z m }. Proof. -(i) Let x = x k β k + • • • + x 0 and y = y l β l + • • • + y 0 denote two elements of Z + β . Then z = x + y is of the form z = a j β j + • • • + a 0 with 0 a j 2 β .
Write now the β -expansion of z as z = +∞ j=1 z -j β -j + e j=0 z j β j and assume it is finite. Then it admits only a β-expansion up till the term indexed by its preperiod k 0 (z) and the period has necessarily the form given above with r(z) = 1 and z -j = 0 as soon as j > k 0 (z).

Then k 0 (z) j=1 z -j β -j = (a j β j + • • • + a 0 ) -( e i=0 z i β i ). This means that the fractional part k 0 (z) j=1 z -j β -j is a polynomial of the type f i=0 b i β i with -β b i 2 β hence with |b i | 2 β . The set F 2 β is finite (Lemma 4.7
) and the set of all possible fractional parts of elements of Z + β is exactly in one-to-one correspondence with a subset of the finite point set

T [0,1),3 β ∩ Z m of Z m . Therefore, there exists a unique g z = f i=0 b i Z i ∈ T [0,1),3 β ∩ Z m such that B π B (g z ) •u B = f i=0 b i β i = k 0 (z)
j=1 z -j β -j = frac(z) . Let us apply the real and complex embeddings of the number field Q(β). It gives:

k 0 (z) j=1 z -j (β (i-1) ) -j = f j=0 b j (β (i-1)
) j for all i = 2, 3, . . . , m. For the real embeddings in particular this implies (Proposition 3.6):

π B,i (g z ) = π B,i ( f j=0 b j Z j ) = f j=0 b j (β (i-1) ) j X i u B,i = k 0 (z) j=1 z -j (β (i-1) ) -j X i u B,i for all i = 1, 2, • • • , s with all z -j 0.
The case of real embeddings will provide a direct computation of the first upper bound L [0;1),2 β of the preperiod and merits to be isolated. Indeed, since in this case 0 < β (i-1) < 1 for all i ∈ {2, 3, • • • , s}, with s assumed 2, and that all the digits z -j are positive, we necessarily have:

X i -1 β (i-1) -j ψ [0,1),3 β as soon as j is large enough. Recall that X i = m-1 k=0 (β (i-1) ) 2k 1/2
. With the definition of L [0;1),2 β , this implies that the sum of the positive terms k 0 (z) j=1 z -j (β (i-1) ) -j cannot contain any term indexed byj with j > L [0;1),2 β . Hence, k 0 (z) L [0;1),2 β . As for the negative real embeddings and the complex embeddings they will provide the second upper bound of the preperiod by the computation of J 2 β : indeed, its calculation gives an upper bound of the number of terms k 0 (z) in the fractional part of z, hence, after reducing frac(z) to the same denominator, which will be β k 0 (z) , we immediately get the result; (ii) (This is reformulation of Theorem 2.4 in [Bu]) First, we have

F β ⊂ F 2 β , second Z + β + Z + β ⊂ Z + β + F 2 β , Z + β -Z + β ⊂ Z β + (F β ∪ -F β ). Since Ω is invariant by inversion and that F β ∪ -F β ⊂ { B π B (g) • u B | g ∈ T [0,+1),2 β ∩ Z m } ∪ { B π B (g) • u B | g ∈ T (-1,0],2 β ∩ Z m } = { B π B (g) • u B | g ∈ T (-1,+1),2 β ∩ Z m } (Lemma 4.7), we deduce the claim. C 4.9. -Let q = min{L (-1,+1),2 β , J 2 β }.
Then, for all x, y ∈ Z β such that x + y and xy have finite β-expansions, the following relations hold:

x + y (resp. x -y) ∈ β -q Z β . Proof. -Indeed, T ⊂ T ∪ (-T ) . Hence Z β + Z β ⊂ Z β + (T ∪ (-T )) . Since T ∪ (-T ) = { B π B (g) • u B | g ∈ T (-1,+1),3 β ∩ Z m } ,
we deduce the exponent q from the definition of L I ,R and from Theorem 4.8.

Computation of an upper bound of the maximal preperiod J R .-We will use the case (ii) of cut-and-project scheme in Theorem 3.8. Let {Z -j } j 0 be the sequence of vectors defined by Z 0 = ( t Q) j Z -j . We denote as usual the algebraic norm of β by N

(β) = N Q(β)/Q (β) = m-1 i=0 β (i) . Recall that a 0 = (-1) m-1 N (β) .
L 4.10. -(i) The following limit holds: lim j→+∞ Z -j = +∞ ; (ii) for all j ∈ N, Z -j ∈ N (β) -j Z m . In particular, if β is a unit of the number field Q(β) , then all the elements Z -j belong to Z m .

Proof. -(i) Since |β (i) | -1 > 1 for all i = 1, 2, . . . , m -1, the inverse operator ( t Q) -1 acts as a dilation by a factor of modulus strictly greater than one on each t Q-invariant subspace F in R m except Ru: all the non-zero components of the vector Z -j (which never belongs to Ru) in the system {V i } i=2,3,••• ,m diverge when j tends to infinity, hence the claim. (ii) Solving the equation where ( k (α)) k 0 is the sequence of digits of the Rényi β-expansion of α and (r (n) 1 , r

Z 0 = ( t Q) 1 Z -1 shows that Z -1 can be written Z -1 = -a -1 0 (a 1 Z 0 + a 2 Z 1 + • • • + a m-1 Z m-2 -Z m-1 ) ∈ N (β) -1 Z m . Since by construction we have Z j = ( t Q) -1 (Z j+1 ) for all j ∈ Z , applying ( t Q) -1 to the last equality clearly gives Z -2 ∈ N (β) -2 Z m and, by induction Z -h ∈ N (β) -h Z m for all h 0 . Now it is classical that β is a unit of Q(β) if and only if N (β) = ±1 establishing the result. T 4.11. -Denote by B R =      x ∈ R m | π B,i (x) ψ [0;1),R+ β 1 -|β (i-1) | m + β m-1 k=0 |β (i-1) | 2k 1/2 1 -|β (i-1) | , i = 2, 3, • • • , m      the cylinder (band) of axis the expanding line Ru and V R = {x ∈ B R | B π B (x) • u B ∈ [0, 1)} the slice of the band B R . Then this slice is such that J R #(V R ∩ N (β) -m Z m ).
(n) 2 , • • • , r (n) m ) ∈ Z m .
Recall that 0 (α) = α = 0 . The real and complex embeddings of the number field Q(β) applied to T n (α) provide the m equalities, with j = 1, 2, • • • , m:

β ( j-1) n •    m-1 i=0 p i β ( j-1) i - n k=1 k (α) β ( j-1) -k    = m k=1 r (n) k β ( j-1) -k . We deduce that m-1 k=0 |β (i-1) | 2k 1/2 π B,i ( m k=1 r (n) k (β ( j-1) ) -k ) = m k=1 r (n) k (β ( j-1) ) -k m-1 i=0 |p i ||β ( j-1) | n+i + β n k=0 |β ( j-1) | k

Appendix.-Classification of Delone sets

We will say that a subset Λ of R n is (i) uniformly discrete if there exists r > 0 such that xy r for all x, y ∈ Λ, x ≠ y, (ii) relatively dense if there exists R > 0 such that, for all z ∈ R n , there exists λ ∈ Λ such that the ball B(z, R) contains λ, (iii) a Delone set if it is relatively dense and uniformly discrete. Delone sets are basic objects for mathematical quasicrystals [START_REF] L | Mathematical Quasicrystals and the problem of diffraction[END_REF] [MVG]. D 4.12. -A cut-and-project scheme consists of a direct product E × D, where E and D are Euclidean spaces of finite dimension, and a lattice L in E × D so that, with respect to the natural projections p 1 : E × D → E , p 2 : E × D → D : (i) p 1 restricted to L is one-to-one onto its image p 1 (L), (ii) p 2 (L) is dense in D. We will denote by * the following operation: Meyer sets were introduced in [Mey]. By definition, we will say that Λ, assumed to be a relatively dense subset of R n , is a Meyer set of R n if it is a subset of a model set. Other equivalent definitions can be found in [Mo] or [Mey] The class of finitely generated Delone sets of R n is strictly larger than the class of Delone set of finite type of R n , which is itself larger than the class of Meyer sets of R n [La] [La1]. T 4.16. -( Lagarias [La]) Let Λ be a Delone set in R n such that ηΛ ⊂ Λ for a real number η > 1. The following assertions hold: (i) If Λ is finitely generated, then η is an algebraic integer. If the rank of Z[Λ] is s, then the degree of η divides s, (ii) If Λ is a Delone set of finite type, then η is a Perron number or is a Lind number.

* := p 2 • ( p 1 | L ) -1 : p 1 (L) → D.
Although Z β is associated with two canonical cut-and-project schemes when β is a noninteger Perron number, the converse of the assertion (ii) of Theorem 4.16 seems to be an open problem. It is at least already related to the question Q1' of the introduction and to various arithmetical and dynamical problems [ABEI].

  and the condition (ii) in Proposition 4.1 satisfied, then the equality #(L (θ)) = 0 holds for θ < θ min and therefore θ f θ min /2 . T 4.5. -(Minimal decomposition).

  T 4.8. -(i) For all x, y ∈ Z + β such that x + y has a finite β-expansion the following relation holds: x + y ∈ β -L Z + β where L := min{L [0,1),2 β , J 2 β }; (ii) the following inclusions hold:

Proof. -

 - Each element α ∈ F R can be written α = m-1 i=0 p i β i with p i ∈ Z and m-1 i=0 p i Z i ∈ T [0,1),R+ β ∩Z m (Lemma 4.7). Thus, |p i | ψ [0,1),R+ β for all i = 0, 1, • • • , m-1. Now ([Sch] and section 2), the following equality holds for all n 0:T n (α) = β n •

11

  -|β ( j-1) | ψ [0,1),R+ β (1 -|β ( j-1) | m ) k ∈ N (β) -m Z m . Its projections by the projection mappings π B,i , i = 2, 3, • • • , m to the Qinvariant subspaces of R mare bounded by constants which are independant of n . The restriction of the lifting of the operator T to V R ∩ N (β) -m Z m has self-avoiding orbits (to have a preperiod) whose length is necessarily smaller than the number of available points in the volume V R . We deduce the upper bound #(V R ∩ N (β) -m Z m ) of J R .

  D 4.13. -A subset Λ of a finite dimensional Euclidean space E is a model set (also called a cut-and-project set) if there exist a cut-and-project scheme (E × D, L) and a subset Ω of D with nonempty interior and compact closure such thatΛ = Λ(Ω) = {p 1 (l ) | l ∈ L, p 2 (l ) ∈ Ω }, equivalently = { v ∈ p 1 (L) | v * ∈ Ω }.The set Ω is called acceptance window.

  . For instance, Λ is a Meyer set if and only if it is a Delone set and there exists a finite set F such that Λ -Λ ⊂ Λ + F ; or if and only if it is relatively dense and Λ -Λ is uniformly discrete. The above definition shows that the class of Meyer sets of R n contains the class of model sets of R n . T 4.14. -( Meyer [Mey]) Let Λ be a Delone set in R n such that ηΛ ⊂ Λ for a real number η > 1. If Λ is a Meyer set, then η is a Pisot or a Salem number. D 4.15. -A Delone set Λ is said to be finitely generated if Z[Λ -Λ] is finitely generated. A Delone set Λ is said to be of finite type if Λ -Λ is such that its intersection with any closed ball of R n is a finite set.

Acknowledgements.-We are indebted to Christiane Frougny, Valérie Berthé, Pierre Arnoux and a referee for very useful and valuable comments and discussions.


Luc Brun

Benoit Gaüzère

Sébastien Bougleux

Florian Yger

A new Sinkhorn algorithm with Deletion and Insertion operations

Introduction

This report is devoted to the continuous estimation of an -assignment(Definition 1). Roughly speaking, an -assignment between two sets V 1 and V 2 may be understood as a bijective mapping between a sub part of V 1 and a sub part of V 2 . The remaining elements of V 1 (not included in this mapping) are mapped onto an pseudo element of V 2 . We say that such elements are deleted. Conversely, the remaining elements of V 2 correspond to the image of the pseudo element of V 1 (Figure 1). We say that these elements are inserted.

Let us note that if V 1 and V 2 have the same size, the bijective mapping induced by an -assignment may involve all elements of V 1 , each element being mapped onto a single element of V 2 . In this sense, an -assignment is more general than a bijective mapping. Moreover, the main advantage of an -assignment is that it provides us the freedom to not map any element which is then assigned to the element of V 2 or belong to the image of the element of V 1 . This last property allows us to reject some mappings if for example, these mappings are associated to a large cost.

An -assignment function may be associated to an -assignment matrix (Figure 1(b)) just like any bijective mapping is associated to a permutation matrix. Given two sets, V 1 and V 2 of respective sizes n and m, an -assignment matrix is encoded by a (n+1)×(m+1) matrix, where n + 1 and m + 1 play respectively the roles of the element of V 1 and the one of V 2 . The last column of index m + 1 of such a matrix encodes the deletions while the last line encodes the insertions. By construction, there is a single 1 in each of the first n rows and m columns, the remaining elements being set to 0. Given V 1 and V 2 , one can define a (n + 1) × (m + 1) cost matrix encoding the cost of the mapping of any element of V 1 onto an element of V 2 as well as the cost of deleting each element of V 1 and inserting each element of V 2 . Finding an -assignment minimizing the sum of mappings, deletions and insertions costs is a direct extension of the Linear Sum Assignment Problem (LSAP) called the Linear Sum Assignment Problem with Edition [START_REF] Bougleux | Linear Sum Assignment with Edition[END_REF] (LSAPE). Given an -assignment matrix X and a cost matrix C, this cost may be formulated as:

min X n+1 i=1 m+1 j=1 c i,j x i,j
where X is taken over all -assignment matrices.

We define in previous works [START_REF] Bougleux | A Hungarian Algorithm for Error-Correcting Graph Matching[END_REF][START_REF] Bougleux | Linear Sum Assignment with Edition[END_REF], an adaptation of the Hungarian algorithm which allows to find an optimal solution to the above problem in O(min(n, m) 2 max(n, m)). However, while providing an optimal solution, this algorithm does not readily allow the computation of the gradient of the associated operation. This last drawback, does not allow to easily insert such an algorithm into a deep learning pipeline. On the other hand, the Sinkhorn algorithm [START_REF] Dennis | concerning nonnegative matrices and doubly stochastic matrices[END_REF], is based on a continuous relaxation of the problem where permutation matrices are replaced by bi-stochastic matrices with an entropic regularization. This algorithm is the workhorse of computational optimal transport [START_REF] Peyré | Computational optimal transport: With applications to data science[END_REF] and is based on iterative matrix multiplications hereby allowing the backpropagation of the gradient [START_REF] Genevay | Learning generative models with sinkhorn divergences[END_REF]. The aim of this technical report is to transpose the results of the Sinkhorn algorithm to assignment matrices. Just like the Sinkhorn algorithm which does not provide a permutation matrix but rather a bi-stochastic matrix, our algorithm will provide an bi-stochastic matrix (Definition 3). This last point may be of advantage within the Neural Network framework where the hard decisions corresponding toassignment matrices may not allow a proper propagation of the gradient.

More formally, given a similarity matrix S (which may be easily deduced from a cost matrix), we aim at finding two diagonal matrices D 1 and D 2 such that D 1 SD 2 is a bi-stochastic matrix. Section 2 provides the main definitions and notations used in the remaining part of this report. The existence and uniqueness of a solution is demonstrated in Section 3 while Section 4 provides a constructive algorithm which convergence is demonstrated. Let us note that while Section 3 is a simple adaption of the original proof [START_REF] Dennis | concerning nonnegative matrices and doubly stochastic matrices[END_REF], Section 4 is significantly different from [START_REF] Dennis | concerning nonnegative matrices and doubly stochastic matrices[END_REF] since the arguments used for bistochastic matrices in the original proof do not hold for bi-stochastic matrices.

Definitions and notations

Definition 1 ( -assignment).

Let n and m be two strictly positive integers. An -assignment is a mapping ϕ : {1, . . . , n + 1} → P({1, . . . , m + 1}) satisfying the following constraints:

   ∀i ∈ {1, . . . , n}, |ϕ(i)| = 1 ∀j ∈ {1, . . . , m}, |ϕ -1 (j)| = 1 m + 1 ∈ ϕ(n + 1)
where P({1, . . . , m + 1}) is the power set of {1, . . . , m + 1}.

Each element of i ∈ {1, . . . , n} is thus mapped onto a set composed of a single element of {1, . . . , m + 1} (|ϕ(i)| = 1) and in the same way the set of antecedents of each j ∈ {1, . . . , m} is reduced to one element (|ϕ -1 (j)| = 1). Hence the only element of {1, . . . , n + 1} which can be mapped onto a set composed of several elements is n + 1. In the same way, m + 1 is the only element which may have several antecedents. The constraint m + 1 ∈ ϕ(n + 1) ensures that n + 1 is mapped to at least one element and that m + 1 has at least an antecedent.

In the example of Figure 1 we have n = 3 and m = 2. Elements 1, 2, 3 are respectively mapped onto {b}, {a}, {3}. Where the last mapping corresponds to a deletion of 3 (which is mapped onto m+1 = 3). Consequently m+1 = 3 has two antecedents 3 and 4 = n+1.

Definition 2 ( -row/column stochastic matrix).

A non negative (n + 1) × (m + 1) matrix X is called an -row stochastic matrix iff:

∀i ∈ {1, . . . , n} m+1 j=1 X i,j = 1 X is called an -column stochastic matrix iff: ∀j ∈ {1, . . . , m} n+1 i=1 X i,j = 1

Definition 3 ( -bi-stochastic matrix).

A non negative (n + 1) × (m + 1) matrix X is called an -bi-stochastic matrix iff: m+1) , X is called an -assignment matrix and there is a one-to-one mapping between -assignments and -assignment matrices.

           m+1 j=1 X i,j = 1 ∀i ∈ {1, . . . , n} n+1 i=1 X i,j = 1 ∀j ∈ {1, . . . , m} x n+1,m+1 = 1 If X ∈ {0, 1} (n+1)×(
Let us note that any -bi-stochastic matrix is a bi-stochastic matrix on which the bi-stochastic constraints are relaxed on the last line and last column. So any squared bi-stochastic matrix is also an -bi-stochastic matrix (the reverse being obviously false).

Definition 4. -diagonal

If A is a (n+1)×(m+1) matrix and ϕ an -assignment then the set A 1,ϕ(1) , . . . , A n+1,ϕ(n+1) is called an -diagonal of A corresponding to ϕ. if A is squared and ϕ is the identity, the diagonal is called the main diagonal.

Note that a 1,ϕ(1) , . . . , a n,ϕ(n) is a sequence (as ϕ(i) is unique for i ∈ {1, . . . , n} while a n+1,ϕ(n+1) is a set). The above definition is a straightforward of the usual notion of diagonal where ϕ is required to be a permutation. In the following we will only consider -diagonals of matrices which will be simply called diagonal.

Definition 5. total support

If A is a nonnegative matrix, A is said to have total support if A = 0 and if every positive element of A lies on a positive -diagonal. A nonnegative matrix that contains a positive diagonal is said to have a support.

If µ and ν define set of indices respectively contained in {1, . . . , n+1} and {1, . . . , m+ 1} then :

• A[µ, ν] denotes the sub matrix of A restricted to indices µ and ν,

• A(µ, ν] denotes the sub matrix of A restricted to indices not contained in µ, i.e. {1, . . . , n + 1} -µ, and to the indices contained in ν,

• A[µ, ν) denotes the sub matrix of A restricted to indices contained in µ and not contained in ν, i.e. {1, . . . , m + 1} -ν.

• A(µ, ν) denotes the sub matrix of A restricted to the indices not contained in µ and ν.

Definition 6. Secable rectangular matrix

A rectangular n × m non negative matrix A is said to be secable if one can find :

• a partition of {1, . . . , n} into two sets X and Y and

• a partition of {1, . . . , m} into two sets Z and T such that:

A[X, T ] = 0 and A[Y, Z] = 0 Z T X A[X,Z] 0 Y 0 A[Y,T]
Let us note that this notion of secable matrix is quite close from the one of block diagonal matrix. However, A[X, T ] is not required to be squared.

Existence and uniqueness

= diag(x 1 , . . . , x n , 1), D 2 = diag(y 1 , . . . , y m , 1) , D 1 = diag(x 1 , . . . , x n , 1) and D 2 = diag(y 1 , . . . , y m , 1). If p i = x i xi and q i = y i yi : n+1 i=1
x i a ij y j = 1, ∀j = 1, . . . , m;

(1)

m+1 j=1 x i a ij y j = 1, ∀i = 1, . . . , n (2)
n+1 i=1 p i x i a ij q j y j = 1, ∀j = 1, . . . , m; (3)
m+1 j=1 p i x i a ij q j y j = 1, ∀i = 1, . . . , n (4)
Let

E j = {i|a ij > 0}, F i = {j|a ij > 0} and put α = {i ∈ {1, . . . , n}|p i = min i p i = p}, β = {j ∈ {1, . . . , m}|q j = max j q j = q}
Let us note that since x n+1 = y m+1 = x n+1 = y m+1 = 1 we have p n+1 = q m+1 = 1, p ≤ 1 and q ≥ 1. Moreover, we have by hypothesis E j ∩ {1, . . . , n} = ∅ and F i ∩ {1, . . . , m} = ∅ for all (i, j) ∈ {1, . . . , n} × {1, . . . , m}.

Let us fist show that α = ∅ ⇐⇒ β = ∅. If α = ∅, then For all i ∈ {1, . . . , n} p i > 1. Then using equation (3), we have for any j ∈ {1, . . . , m}:

1 = n+1 i=1 p i x i a ij q j y j > q j n+1 i=1
x i a ij y j = q j Hence q j < 1 for all j ∈ {1, . . . , m} and thus β = ∅.

Conversely, if β = ∅, we have q j < 1 for all j ∈ {1, . . . , m} . Using equation ( 4) for i ∈ {1, . . . , n}:

1 = m+1 j=1 p i x i a ij q j y j < p i m+1 j=1
x i a ij y j = p i Hence, p i > 1 for all i ∈ {1, . . . , n} and α = ∅.

In this case we consider the alternative definitions for α and β: α = {i ∈ {1, . . . , n}|p i = max i p i = p}, β = {j ∈ {1, . . . , m}|q j = min j q j = q} Since p i > 1 for all i ∈ {1, . . . , n} and q j < 1 for all j ∈ {1, . . . , m}, α and β are non empty.

Using initial definitions for α and β, let us thus consider i 0 ∈ α and j 0 ∈ β. Then using (3):

q j0 = 1 n+1 i=1 p i x i a ij0 y j0 ≤ 1 p i0 n+1 i=1 x i a ij0 y j0 = p -1 i0
where the last equality comes from (1). Similarly, using (4):

p i0 = 1 m+1 j=1 x i a ij q j y j ≥ 1 q j0 m+1 j=1 x i a ij y j = q -1 j0
where the last equality comes from (2). Whence q j0 = p -1 i0 = p -1 . But in this case, we have using (3):

n+1 i=1 p i x i a ij0 q j0 y j0 = i∈Ej 0 p i x i a ij0 q j0 y j0 = i∈Ej 0 p i p x i a ij0 y j0 = 1
This last equality is compatible with (1) only if p i = p for all i ∈ E j0 . Dropping sub indices, we have for all j ∈ β and all i ∈ E j p i = p. Thus

j∈β E j ⊆ α ∪ {n + 1}. Hence if j ∈ β and i ∈ α ∪ {n + 1}, i ∈ j∈β E j . So a ij = 0.
More concisely, we have:

A(α ∪ {n + 1}, β] = 0.
In the same way, p i0 = q -1 j0 = q -1 implies using (4):

m+1 j=1 p i0 x i0 a i0j q j y j = j∈Fi 0 q j q x i0 a i0j y j = 1
which is compatible with (2) only if q j = q for all j ∈ F i0 . Thus for all i ∈ α and for all j ∈ F i we have q j = q. Thus i∈α

F i ⊂ β ∪ {m + 1} and A[α, β ∪ {m + 1}) = 0
On α × β we have p i q j = pq = 1. Thus:

a ij = b ij x i y j = b ij x i y j ⇒ p i q j b ij = b ij ⇒ b ij = b ij Hence B[α, β] = B [α, β]
Moreover, for any j ∈ β we have:

n+1 i=1 b i,j = n+1 i=1 x i a i,j y j = i∈Ej ⊂α∪{n+1} x i a i,j y j = i∈α∪{n+1} x i a i,j y j = i∈α∪{n+1} b i,j = 1
In the same way, we have for any j ∈ β, i∈α∪{n+1} b i,j = 1. But in this case using

B[α, β] = B [α, β],
we have for j ∈ β:

1 = i∈α∪{n+1} b i,j = i∈α b i,j + b n+1,j = i∈α b ij + b n+1,j
Thus 1 -b n+1,j + b n+1,j = 1 which induces b n+1,j = b n+1,j which imposes q = 1. Indeed, since j ∈ β, we have q j = q and q j b n+1,j = qb n+1,j = b n+1,j . In the same way for i ∈ α:

m+1 j=1 b ij = m+1 j=1 x i a i,j y j = j∈Fi⊂β∪{m+1} x i a i,j y j = j∈β∪{m+1} x i a i,j y j = j∈β∪{m+1} b i,j = 1
and we have in the same way for i ∈ α, j∈β∪{m+1} b i,j = 1. We then obtain for i ∈ α:

j∈β∪{m+1} b i,j = j∈β b i,j + b i,m+1 = j∈β b i,j + b i,m+1
Using the previous equality 1

-b i,m+1 + b i,m+1 = 1, hence b i,m+1 = b i,m+1 and p = 1 (since pb i,n+1 = b i,n+1 ). Thus B[α ∪ {n + 1}, β ∪ {m + 1}] = B [α ∪ {n + 1}, β ∪ {m + 1}
] is an bi-stochastic matrix (where n + 1 and m + 1 plays the role of the last row and column respectively).

Let us briefly show that {1, . . . , n} -α and {1, . . . , m} -β are simultaneously empty or non empty. Let us fist suppose that α = {1, . . . , n} and let us consider j ∈ β ∪{m+1}. Since E j ∩ {1, . . . , n} = ∅ by hypothesis, it exists i ∈ {1, . . . , n} = α such that j ∈ F i . But since i ∈ α, we have F i ⊂ β ∪ {m + 1} and thus a contradiction. In the same way, if β = {1, . . . , m}, let us consider i ∈ α ∪ {n + 1}. Since F i ∩ {1, . . . , m} = ∅, it exists j ∈ {1, . . . , m} = β such that i ∈ E j ⊂ α ∪ {n + 1}. Again a contradiction.

If A is non secable the configuration where both {1, . . . , n} -α and {1, . . . , m} -β are non empty correspond to a partition of {1, . . . , n + 1} into α ∪ {n + 1} and its complementary and a partition of {1, . . . , m + 1} into β ∪ {m + 1} and its complementary with no connections between α ∪ {n + 1} and the complementary of β ∪ {m + 1} nor any connection between β ∪ {m + 1} and the complementary of α ∪ {n + 1} (Figure 2). Such a decomposition being refused, we have α = {1, . . . , n} and β = {1, . . . , m} Hence

A[α ∪ {n + 1}, β ∪ {m + 1}] = A and D 1 AD 2 = D 1 AD 2 and D 1 and D 2 are unique (p = q = p n+1 = q m+1 = 1).
If the non secable property of A does not hold and A(α ∪ {n + 1}, β] and A[α, β ∪ {m + 1}) exist, B(α, β) and B (α, β) exist, include the row n + 1 and the column m + 1, are bi-stochastic matrices and have a size lower than the one of A. Furthermore,

B(α, β) = D 1 A(α, β)D 2 and B (α, β) = D 1 A(α, β)D 2
where D i and D j have like D i and D j (from which they are derived) a positive main diagonal with a 1 at last position. The argument may be repeated on these submatrices until

D 1 AD 2 = D 1 AD 2 is established. Given that B(α, β) = B (α, β), we already know that B[α ∪ {n + 1}, β ∪ {m+1}] = B [α∪{n+1}, β ∪{m+1}
] and that A and hence B and B are zeros elsewhere (Figure 2). Hence B is equal to B . Note however, that since α∪{n+1} = {1, . . . , n+1} (and the same for β) D 1 and D 2 are no longer unique.

β m+1 α 0 0 B(α, β) n + 1

A constructive algorithm

For any i ∈ {1, . . . , n} and any j ∈ {1, . . . , m} let us consider the series (x i,p ) p∈N and (y j,p ) p∈N defined as follows:

∀i ∈ {1, . . . , n} x i,0 = m+1 j=1 a ij -1 x i,p+1 = χ -1
i,p x i,p ∀j ∈ {1, . . . , m} y j,0 = 1 y j,p+1 = γ -1 j,p y j,p Moreover we also define:

   ∀i ∈ {1, . . . , n} χ i,p = m+1 j=1 x i,p a i,j y j,p = m+1 j=1 γ -1 j,p-1 x i,p a i,j y j,p-1 ∀j ∈ {1, . . . , m} γ j,p = n+1 i=1 χ -1
i,p x i,p a i,j y j,p with for all p:

x n+1,p = y m+1,p = 1 χ n+1,p = γ m+1,p = 1

Let us denote by L p the (n+1)×(m+1) matrix whose entries in n×(m+1) are equal to (x i,p a i,j y j,p-1 ) and whose last row is filled with zeros but a 1 at position (n + 1, m + 1).

In the same way let us denote by C p the (n + 1) × (m + 1) matrix whose entries in (n + 1) × m are equal to x i,p a i,j y j,p and whose last column is filled with zeros except a 1 at position (n + 1) × (m + 1).

If χ p ∈ R n+1 and γ p ∈ R m+1 denote the vectors encoding respectively (χ i,p ) i∈{1,...,n+1} and (γ j,p ) j∈{1,...,m+1} we have for p ≥ 1:

χ p = L p γ -1 p-1 γ p = C T p χ -1 p (5)
L p is row stochastic. Indeed, for any i ∈ {1, . . . , n} and p ≥ 1:

m+1 j=1 (L p ) i,j = m+1 j=1 x i,p a i,j y j,p-1 = m+1 j=1 χ -1 i,p-1 x i,p-1 a i,j y j,p-1 = χ -1 i,p-1 m+1 j=1 x i,p-1 a i,j y j,p-1 = 1
and the last line of L p contains a single entry equal to 1.

Moreover, C p is column stochastic for p ≥ 1. Indeed for each j ∈ {1, . . . , m}:

n+1 i=1 (C p ) i,j = n+1 i=1 x i,p a i,j y j,p = n+1 i=1 χ -1 i,p-1 γ -1 j,p-1 x i,p-1 a i,j y j,p-1 = γ -1 j,p-1 n+1 i=1 χ -1 i,p-1 x i,p-1 a i,j y j,p-1 = 1 L p =        x 1,p a 1,1 y 1,p-1 x 1,p a 1,2 y 2,p-1 . . . x 1,p a 1,m y m,p-1 x 1,p a 1,m+1 x 2,p a 2,1 y 1,p-1 x 2,p a 2,2 y 2,p-1 . . . x 2,p a 2,m y m,p-1 x 2,p a 2,m+1
. . . Note that C 0 is not column stochastic. One noticeable effect of this negative property is that γ 0 = 1 while χ 0 = 1. Moreover the last column contains a single positive entry equal to 1. Hence C T p is row stochastic and equation 5 involves two row stochastic matrices.

x n,p a n,1 y 1,p-1 x n,p a n,2 y 2,p-1 . . . x n,p a n,m y m,p-1 x n,p a n,m+1 0 0 . . . 0 1        C p =        x 1,p a 1,
       C T p =        x 1,p a 1,
Combining both equations of 5 we have:

χ p = L p C T p-1 χ -1 p-1 -1
where the inverse notation applied to a vector denotes the element-wise inverse operation. Since L p is row stochastic, we have:

χ p -1 = L p C T p-1 χ -1 p-1 -1 -1 Using 1 x -1 = 1-x
x we obtain:

χ p -1 = L p 1 -C T p-1 χ -1 p-1 C T p-1 χ -1 p-1 -1
where is the element-wise product also known as Hadamard product. Since, for p ≥ 2, C T p-1 is row stochastic we have C T p-1 1 = 1 and thus:

χ p -1 = L p C T p-1 1 -χ -1 p-1 C T p-1 χ -1 p-1 -1 = L p diag C T p-1 χ -1 p-1 -1 C T p-1 1 -χ -1 p-1 Using 1 -1 x = x-1
x we obtain:

χ p -1 = L p diag C T p-1 χ -1 p-1 -1 C T p-1 (χ p-1 -1) χ -1 p-1 = L p diag C T p-1 χ -1 p-1 -1 C T p-1 diag(χ p-1 ) -1 (χ p-1 -1) = L p diag (γ p-1 ) -1 C T p-1 diag(χ p-1 ) -1 (χ p-1 -1)
Moreover the left and right multiplications of C T p-1 by diagonal matrices is equivalent to a multiplications of its lines by γ -1 p-1 and its columns by χ -1 p-1 . More precisely we have:

diag (γ p-1 ) -1 C T p-1 diag(χ p-1 ) -1 =        χ -1 1,p-1 γ -1 1,p-1 x 1,p-1 a 1,1 y 1,p-1 . . . χ -1 n,p-1 γ -1 1,p-1 x n,p-1 a n,1 y 1,p-1 γ -1 1,p-1 a n+1,1 y 1,p χ -1 1,p-1 γ -1 2,p-1 x 1,p-1 a 1,2 y 2,p-1 . . . χ -1 n,p-1 γ -1 2,p-1 x n,p-1 a n,2 y 2,p-1 γ -1 2,p-1 a n+1,2 y 2,p . . . χ -1 1,p-1 γ -1 m,p-1 x 1,p-1 a 1,m y m,p-1 . . . χ -1 n,p-1 γ -1 m,p-1 x n,p-1 a n,m y m,p-1 γ -1 m,p-1 a n+1,m y m,p 0 . . . 0 1        =        x 1,p a 1,
       = C T p
Hence we have:

∀p ≥ 2 χ p -1 = L p C T p (χ p-1 -1) (6)
As L p and C T p are row stochastic matrices, so is L p C T p . Moreover, L p C T p is a square (n + 1) × (n + 1) matrix. Let us note that due do this row stochastic property the last equality is equivalent to: ∀p ≥ 2 χ p = L p C T p χ p-1 Examining more precisely the matrix L p C T p we have:

       (L p C T p ) i,j = x i,p x j,p m k=1 a i,k a j,k y k,p y k,p-1 ∀(i, j) ∈ {1, . . . , n} 2 (L p C T p ) n+1,j = 0 ∀j ∈ {1, . . . , n} (L p C T p ) i,n+1 = x i,p m k=1 a i,k a n+1,k y k,p y k,p-1 + x i,p a i,m+1 ∀i ∈ {1, . . . , n} (L p C T p ) n+1,n+1 = 1
Let us note that i and j may be interchanged in the first equation above. Hence:

∀(i, j) ∈ {1, . . . , n} 2 (L p C T p ) i,j = (L p C T p ) j,i
Moreover we have for all (i, j) ∈ {1, . . . , n} × {1, . . . , m}:

   χ i,0 = x i,0 m+1 j=1 a i,j y j,0 = m+1 j=1 ai,j m+1 j=1 ai,j = 1 γ j,0 = n+1 i=1 χ -1 i,0 x i,0 a i,j y j,0 = n+1 i=1 ai,j m+1 j=1 ai,j
Since the sum of a line (or a column) of A cannot be equal to zeros, it exists two positive numbers γ, γ such that:

0 < γ1 ≤ γ 0 ≤ γ1
We have thus:

0 < γ1 ≤ γ 0 ≤ γ1 γ -1 1 ≤ γ -1 0 ≤ γ -1 1 γ -1 1 ≤ L 1 γ -1 0 ≤ γ -1 1 0 < γ -1 1 ≤ χ 1 ≤ γ -1 1
where the last inequality is deduced from the fact that χ 1 = L 1 γ -1 0 , and the fact that all the entries of L 1 are non negative.

Hence using ( 6) and a basic recursion we have :

∀p ≥ 1, 0 < γ -1 1 ≤ χ p ≤ γ -1 1
Using γ p = C T p χ -1 p which is row stochastic we obtain:

∀p ≥ 1, 0 < γ1 ≤ γ p ≤ γ1
From now on, let us suppose that all entries of the last column and the last line of A are positive.

Let us suppose that for some i ∈ {1, . . . , n} we have lim p→+∞ x i,p = +∞. Since we have:

χ i,p = x i,p   m j=1
a i,j y j,p + a i,m+1   we also have lim p→+∞ χ i,p = +∞ which is impossible since χ is bounded. It exists thus an upper bound M x such that:

∀i ∈ {1, . . . , n + 1}∀p ≥ 2 x i,p ≤ M x
In this case we have:

1 = y j,p+1 n+1 i=1 a i,j x i,p+1 ≤ y j,p+1 M x n+1 i=1 a i,j ≤ y j,p+1 M x A 1
Hence we have:

y j,p+1 ≥ 1 M x A 1
In the same way, let us suppose that for some j ∈ {1, . . . , m} we have lim p→+∞ y j,p = +∞. Since we have:

γ j,p = y j,p n i=1 χ -1 i,p a i,j x i,p + a n+1,j
we also have lim p→+∞ β j,p = +∞ which is forbidden since γ is also upper bounded. So, it exists an upper bound M y such that: ∀j ∈ {1, . . . , m + 1}, ∀p ≥ 2 y j,p ≤ M y

In this case we have:

1 = x i,p+1 m+1 j=1 a i,j y j,p ≤ x i,p+1 M y m+1 j=1 a i,j ≤ x i,p+1 M y A ∞ Hence: x i,p+1 ≥ 1 M y A ∞ 1.
The series (x i,p ) p and (y j,p ) p being both lower bounded, all non zeros entries of L p C T p are lower bounded by a positive value. Moreover since L p C T p is row stochastic, all its entries are bounded by 1. We say that the entries of L p C T p are uniformly positive.

Moreover:

∀i ∈ {1, . . . , n} (L p C p ) i,i = x 2 i,p m k=1 (a i,k ) 2 y k,p y k,p-1
Since both series (x i,p ) p and (y j,p ) p are lower bounded and that a line of A cannot be equal to zeros, we have (L p C T p ) i,i > 0. Let us additionally note that (L p C T p ) n+1,n+1 = 1 > 0.

3. Let us consider the Graph G(L p C T p ) = (V, E) where V = {1, . . . , n + 1} and an edge (i, j) connects node i to node j iff (L p C T p ) j,i > 0. In this case since for any i ∈ {1, . . . , n}, (L p C T p ) i,n+1 ≥ x i,p a i,m+1 > 0. The node n + 1 is adjacent to all nodes in {1, . . . , n}. Conversely, since (LC T p ) n+1,j = 0 for all j ∈ {1, . . . , n}, no node j is incident to n + 1. The node n + 1 defines a source component. Moreover, for any pair (i, j) ∈ {1, . . . , n} 2 , since (L p C T p ) i,j = (L p C T p ) j,i , node i is adjacent to node j and vice versa. Hence n + 1 is the only source component of G(L p C T p ) which is thus quasi-strongly connected.

Using [START_REF] Proskurnikov | Recurrent averaging inequalities in multi-agent control and social dynamics modeling[END_REF], Lemma 3 with T=1 we can conclude that χ converges towards a consensus of the form c1. Since the last entry of χ is a constant equal to 1 we have c = 1 and:

lim p→+∞ χ p = 1 Since C T
p is row stochastic we have:

lim p→+∞ γ p = lim p→+∞ C T p χ -1 p = 1 Since χ i,p = xi,p
xi,p+1 converges towards 1, (x i,p ) p is a Cauchy serie in a complete space(R + ). Hence (x i,p ) p converges. The same argument holds for (y j,p ) p .

Let us consider the two diagonal matrices D 1,p = diag(x 1,p , . . . , x n+1,p ) and D 2,p = diag(y 1,p , . . . , y m+1,p ) together with the matrix S p = D 1,p AD 2,p . By construction we have for any i ∈ {1, . . . , n}:

(S p 1) i = m+1 j=1
x i,p a i,j y j,p = χ i,p Thus: lim p→∞ (S p 1) i = lim p→∞ χ i,p = 1. The matrix S p converges thus towards an row stochastic matrix.

Moreover, for any j ∈ {1, . . . , m}:

(S T p 1) j = n+1 i=1
x i,p a i,j y j,p = n+1 i=1

x i,p+1 a i,j y j,p + (x i,p -x i,p+1 )a i,j y j,p Since x i,p converges for any i ∈ {1, . . . , n + 1}, y j,p is bounded, and the above sums are finite, it exists for any η > 0 a value p 0 such that for any p ≥ p 0 , we have :

n+1 i=1 x i,p a i,j y j,p - n+1 i=1
x i,p+1 a i,j y j,p < η Both sums converge (or diverge) thus toward a same value. Moreover:

n+1 i=1
x i,p+1 a i,j y j,p = γ j,p

We have thus: lim p→∞ (S T p 1) j = lim p→∞ γ j,p = 1. The matrix S p converges thus toward an column stochastic matrix and hence an bi-stochastic matrix.

Two iterative algorithms

The code (in python) corresponding to the construction of matrices D 1,p and D 2,p is provided in Figure 4. You may note the fact that we set x n+1 and y m+1 to 1 respectively on line 9 and 15. This point together with the use of rectangular matrices is the main Figure 5: A version of our computation of an assignment matrix which updates directly the matrix S p difference between this algorithm and the classical Sinkhorn algorithm. The convergence criterion which allows to avoid to loop up to the maximum number of iterations is based on the fact that both χ p and γ p converge toward a vector of 1.

An equivalent code computing directly the matrix S p is provided in Figure 5. In this case the setting of x n+1 and y m+1 to 1 is performed on line 9 and 19. The stopping criterion is based on the computation of the distance of the current matrix to the set of assignment matrices. To do so, we compute the distance between the vector of 1 and S T p 1 after each row normalization(line 11). In the same way, we compute the distance between a vector of 1 and S p 1 after each column normalization (line 15). After convergence, we apply a last row normalization before setting the last line of our assignment matrix to the complement to 1 of each column.

From similarity to cost matrices and vice versa

The Sinkhorn algorithm is well known for providing an approximation of the Linear Sum Assignment Problem (Section 7) which can be formulated as:

max X n i=1 n j=1 s i,j x i,j
where S = (s i,j ) is our similarity matrix and X = (x i,j ) is taken over all bi stochastic matrices. The optimal solution being a permutation matrix, hence a binary matrix.

This maximization problem may be translated into a minimization problem by considering the matrix c1 n×n -S, where 1 n×n is a n × n matrix filled of 1 and c is a positive constant greater than all values of S. We have indeed:

n i=1 n j=1 (c -s i,j )x i,j = n i=1 n j=1 cx i,j - n i=1 n j=1 s i,j x i,j = c n i=1 n j=1 x i,j - n i=1 n j=1 s i,j x i,j = c n i=1 1 - n i=1 n j=1 s i,j x i,j = cn - n i=1 n j=1 s i,j x i,j
Hence c and n being constant, minimize n i=1 n j=1 (c-s i,j )x i,j is equivalent to maximize n i=1 n j=1 s i,j x i,j . The matrix c1 n×n -S is usually interpreted as a cost matrix. This last point is important if one wants to compare the Sinkhorn algorithm to an optimal Hungarian algorithm which performs a minimization of costs instead of a maximization of similarities.

As stated in Section 1, our algorithms being an extension of the Sinkhorn algorithm we expect them to converge to :

max x n+1 i=1 m+1 j=1 s i,j x i,j
where X = (x i,j ) is taken over all -bi stochastic matrices. However, the transformation of this maximization of similarities into a minimization of costs, is slightly more complex in the case of assignment matrices. To do so, let us consider a (n + 1) × (m + 1) matrix C = (c i,j ) with:

c i,j =        2c if i ≤ n ∧ j ≤ n c lr if i = n + 1 ∧ j ≤ m c lc if i ≤ n ∧ j = m + 1 0 if i = n + 1 ∧ j = m + 1
where c, c lr , c lc are three positive constants. Considering the cost matrix C -S we have:

n+1 i=1 m+1 j=1 (c i,j -s i,j )x i,j = n+1 i=1 m+1 j=1 c i,j x i,j - n+1 i=1 m+1 j=1 s i,j x i,j
If c lc = 2c and c lr = 0: We have:

n+1 i=1 m+1 j=1 c i,j x i,j = n i=1 m+1 j=1 2cx i,j = 2c n i=1 1 = 2cn
If c lr = 2c and c lc = 0: We have:

n+1 i=1 m+1 j=1 c i,j x i,j = m j=1 n+1 i=1 2cx i,j = 2c m j=1 1 = 2cm
If c lr = c and c lc = c: We have:

n+1 i=1 m+1 j=1 c i,j x i,j = 2c n i=1 m j=1 x i,j + c n i=1 x i,m+1 + c m j=1 x n+1,j = c n i=1 m+1 j=1 x i,j + c m j=1 n+1 i=1 x i,j = cn + cm
In all cases we have thus:

n+1 i=1 m+1 j=1 (c i,j -s i,j )x i,j = Q - n+1 i=1 m+1 j=1 s i,j x i,j (7)
where Q = cn if c lc = 2c and c lr = 0, Q = cm if c lr = 2c and c lc = 0 and finally Q = cn + cm if c lr = c lc = c. The minimization of the left part of equation 7 (minimization of costs) is thus equivalent to a maximization of the similarities. Let us note that the trivial solution consisting to take C = c1 (n+1)×(m+1) does not provide an equivalence between both problems since additional terms related either to the last column or the last row forbid to state that one problem is equal to a constant minus the other problem. Moreover, the last solution ( c lr = c and c lc = c) is the only one allowing to ensure that all coefficients of the similarity matrix are positive when transforming a cost matrix into a similarity matrix.

Experiments

We proposed in Section 4 two algorithms converging toward an unique solution if the conditions defined in Section 3 are satisfied. The aim of this section is to measure experimentally the convergence of our algorithms toward a solution maximizing : n+1 i=1 m+1 j=1 s i,j x i,j over all bi stochastic matrices X. Where S is the input matrix. Such a problem is called a Linear Sum Assignment Problem with Edition (LSAPE).

Deviation of the Sinkhorn algorithm from the optimal solution

Sinkhorn algorithm provides an approximate solution to the well known Linear Sum Assignment Problem (LSAP): where X is taken over the set of bi stochastic matrices. From a certain point of view, LSAP may be considered as a restriction of the LSAPE with squared matrices and no deletions/insertions. Let us first evaluate the error induced by the use of the Sinkhorn algorithm. To do this, we define matrices filled by random number in the interval [START_REF] Bougleux | Linear Sum Assignment with Edition[END_REF][START_REF] Bougleux | A Hungarian Algorithm for Error-Correcting Graph Matching[END_REF]. For each matrix size we compute 100 matrices and compute for each matrix both the solution produced by the Sinkhorn algorithm and the optimal one produced by an Hungarian algorithm. The results of this experiment are displayed in Figure 6 for n ∈ {10, . . . , 200}.

max X n i=1 m j=1 s i,j x i,j
Considering Figure 6, the error is approximately constant for all sizes of matrices lies between 20% and 23%.

Deviation of our algorithms from the optimal solution

In order to test our algorithm we use the same kind of random matrices but with a specific procedure for the last row and column which encode respectively the affinity of each element toward insertions and deletions:

S i,j := rand()+1 if i<n ∧ j<m 0 if i=n ∧ j=m h rand() else (8)
For h ≤ 0.5 we can insure that for any (i, j) ∈ {1, . . . , n} × {1, . . . , m}, s i,j ≥ s n+1,j +s j,m+1 . In other terms we always get a greater sum by substituting i onto j than by deleting i and then inserting j. Conversely, if s i,j < s n+1,j + s j,m+1 the substitution of i onto j will never be part of an optimal -assignment since this operation can be replaced, with a greater value of the sum, by the removal of i and the insertion of j.

Let us first focus on values of h lower than .5. Figure 7 shows on the first two columns the relative error of both sinkhorn Sp and sinkhorn D1D2 according to an optimal LSAPE algorithm [START_REF] Bougleux | A Hungarian Algorithm for Error-Correcting Graph Matching[END_REF] for increasing sizes of the matrix. Let us note that [START_REF] Bougleux | A Hungarian Algorithm for Error-Correcting Graph Matching[END_REF] minimizes a sum of costs. We compare both algorithms using the results of Section 6. For each matrix size 100 random matrices are generated and the results are averaged for all three algorithm (sinkhorn Sp, sinkhorn D1D2 and the optimal one). Our algorithms provide an approximation of the LSAPE which is slightly above 20%, hence comparable with the one provided by the Sinkhorn algorithm for the LSAP problem. Interestingly, our algorithm provides better approximations for small matrix sizes while Figure 6 suggest an opposite behavior for the Sinkhorn algorithm.

The last column of Figure 7 allows to compare more precisely sinkhorn Sp and sinkhorn D1D2. Both algorithms seems to be equivalent since the relative error between both methods do not exceed 1%.

For h greater than .5 we observe in Figure 8 that we keep an error of about 20% for h = 1.0 and h = 2.0. In these cases the values of the last line and the last column remain comparable with the inner values of the random matrix. However, for larger values of h, namely h ∈ {4.0, 6.0, 8.0} we observe a large increase of the relative error especially in the case of squared matrices. We can conclude from these experiments that our algorithms do not converge to the expected value when the values of the last column/line are very large compared to the inner values. More precisely, when we have: s i,j s n+1,i + s i,m+1 for (i, j) ∈ {1, . . . , n} × {1, . . . , m} A simple solution to fix this problem, consists in simplifying the similarity matrix by removing (setting to a low value) any entry (i, j), i ∈ {1, . . . , n} and j ∈ {1, . . . , m} such that (i, j) cannot belong to any optimal solution. Given the similarity matrix S, such entries are characterized by s i,j < s n+1,j + s i,m+1 . In such cases, the substitution of i onto j may be advantageously replaced by the removal of i and the insertion of j. This last point forbids the assignment of i onto j in any optimal -assignment. Figure 9 represents the relative errors according to the optimal solution performed by Sinkhorh Sp and sinkhorn D1D2 using this simplification of the similarity matrix. In this experiment the low value replacing any entry of the matrix S which can not be included in any optimal solution has been fixed to 10 -4 . One can first observe that we get the same behavior for the squared (n × n) case and the rectangular one (n × 2n). We can further observe that all errors remain below 20% for all sizes. Moreover, the relative error appear to be decreasing as a function of h for both algorithms. This may be explained by the fact that as h get higher, the simplified matrix becomes more and more trivial. Indeed for largest values of h, simplified similarity matrices correspond to trivial matrices with a constant value (equal to 10 -4 in this experiment) for all entries (i, j) in {1, . . . , n}×{1, . . . , m} and a last row and column which remains unchanged and greater than 10 -4 by several orders of magnitude. In such cases our algorithms converge immediately to the optimal solution which correspond to the removal of all elements in {1, . . . , n} and the insertion of all elements in {1, . . . , m}.

Execution times of our algorithms

The execution times of both sinkhorn Sp and sinkhorn D1D2 computed either on a graphic card (Nvidia Quadro P2000) or on a CPU (Intel Core i5 650@3.2GHz) are displayed in Figures 10 for n×n squared matrices and in Figure 11 for n×2n rectangular matrices. In both figures and for each size, both algorithms have been run 100 times and the execution times have been averaged.

Considering the squared case (Figure 10) we can see that on graphics cards both sinkhorn Sp and sinkhorn D1D2 take roughly the same amount of time. This is clearly not the case on CPU where the updates of the matrix S p instead of the two diagonal matrices D 1 and D 2 induce a large difference, by a factor greater than 10, between the execution times of sinkhorn Sp and sinkhorn D1D2. We can further observe that if we do not observe major differences between the cuda and cpu execution times for sinkhorn D1D2 this is clearly not the case for sinkhorn Sp which is more than 10 times accelerated by the use of the graphic card. Finally, we can note that the execution times are decreasing as a function of h. As previously, this last point is due to the fact that as h get higher the simplified similarity matrices become more and more trivial and our iterative algorithms need less and less iterations to converge.

Concerning the rectangular case (Figure 11), we observe the same trends than in the squared case. However the ratio between the execution times of sinkhorn Sp and sinkhorn D1D2 on CPU is, in this case, about 6. Moreover, we observe for such matrices a factor approximately equal to 2 between the execution times of sinkhorn D1D2 on GPU and CPU. As in the squared case the CUDA implementation provides a large speedup for sinkhorn Sp. Finally, the execution times of both algorithms decrease as h increase. This phenomenon already encountered in the squared case is due to the same reasons.

Conclusion

We have presented in this report (Section 5) two algorithms. The proof of their convergence is provided in Section 4 while conditions of the existence and uniqueness of the limits are provided in Section 3. Section 6 provides simple methods to transform the sum maximization problem addressed by our algorithms into sum minimization (minimization of a sum of costs).

As shown in Section 7, these algorithms provide an approximate solution to the Linear Sum Assignment Problem with Edition (LSAPE). The relative error of these algorithms compared to the optimal solutions is similar, and even much lower in some cases, to the relative error between the classical Sinkhorn and the optimal solutions to the Linear Sum Assignment Problem (LSAP). The main difference between these algorithms and the Hungarian based algorithms providing the optimal solution is that our algorithms are iterative and differentiable and may thus be easily inserted within a backpropagation based learning framework such as artificial neural networks.

Compared to the LSAP, the LSAPE problem allows to manage assignments between sets of different sizes by allowing the possibility to reject some elements from the matching through the use of insertion/deletion operations. Let us note that these insertion/deletion operations are integrated to the matching algorithm. The LSAPE thus avoids any artificial preprocessing step consisting in the selection of the more promising objects of both sets for matching. Our algorithms output an bi-stochastic matrix which mainly differs from the results provided by [START_REF] Peyré | Computational optimal transport: With applications to data science[END_REF][START_REF] Genevay | Learning generative models with sinkhorn divergences[END_REF] through the explicit notions of insertion and deletion.

Let us finally note that we have required in Section 5 that all the coefficients of the last line and column of the similarity matrix should be positive. This requirement is certainly too strong and weaker conditions may certainly be established by future works. An alternative future research direction, consists in addressing the LSAPE problem using the recent advances in Graph Neural Networks (GNN).

Figure 1 :

 1 Figure 1: (a) An example of -assignment function. 1 is mapped onto b, 2 onto a, 3 is deleted. (b) its associated assignment matrix

Figure 2 :

 2 Figure 2: Decomposition of matrix B.

1 dFigure 4 :18

 14 Figure 4: The code corresponding to the construction scheme described in Section 4

Figure 6 :

 6 Figure 6: Relative Error between the results provided by the Sinkhorn Algorithm and the optimal solution provided by an Hungarian method.

Figure 7 :

 7 Figure 7: Relative errors of sinkhorn Sp and sinkhorn D1D2 according to the optimal solution.

Figure 8 :Figure 9 :Figure 10 :

 8910 Figure 8: Relative errors of sinkhorn Sp and sinkhorn D1D2 according to the optimal solution.

Figure 11 :

 11 Figure 11: Execution times of sinkhorn Sp and sinkhorn D1D2 on n × 2n rectangular matrices using both graphic card and cpu computation.

Theorem 3.1. Let

  A be a nonnegative (n+1)×(m+1) matrix such that A[{1, . . . , n}, {1, . . . , m}] does not contain any line or column filled with 0. A necessary and sufficient condition that there exists an bi-stochastic matrix B of the form D 1 AD 2 where D 1 and D 2 are diagonal matrices with positive main diagonals and a last entry equal to 1 is that A has total support. If B exists then it is unique. Also D 1 and D 2 are unique if and only if A is non secable.Proof. Let us suppose that B = D 1 AD 2 and B = D 1 AD 2 are -bi-stochastic matrices where D 1

  1 y 1,p x 1,p a 1,2 y 2,p . . . x 1,p a 1,m y m,p 0 x 2,p a 2,1 y 1,p x 2,p a 2,2 y 2,p . . . x 2,p a 2,m y m,p 0 . . . x n,p a n,1 y 1,p x n,p a n,2 y 2,p . . . x n,p a n,m y m,p 0 a n+1,1 y 1,p a n+1,2 y 2,p . . . a n+1,m y m,p 1

  1 y 1,p x 2,p a 2,1 y 1,p . . . x n,p a n,1 y 1,p a n+1,1 y 1,p x 1,p a 1,2 y 2,p x 2,p a 2,2 y 2,p . . . x n,p a n,2 y 2,p a n+1,2 y 2,p . . . x 1,p a 1,m y m,p x 2,p a 2,m y m,p . . . x n,p a n,m y m,p a n+1,m y m,p

					
					     
	0	0	. . .	0	1
		Figure 3: L

p and C p matrices.

  1 y 1,p x 2,p a 2,1 y 1,p . . . x n,p a n,1 y 1,p a n+1,1 y 1,p x 1,p a 1,2 y 2,p x 2,p a 2,2 y 2,p . . . x n,p a n,2 y 2,p a n+1,2 y 2,p . . . x 1,p a 1,m y m,p x 2,p a 2,m y m,p . . . x n,p a n,m y m,p a n+1,m y m,p

	0	0	. . .	0	1