Jérémy Berthomieu
email: jeremy.berthomieu@lip6.fr

Jean-Charles Faugère
email: jean-charles.faugere@inria.fr

Guessing Linear Recurrence Relations of Sequence Tuples and P-recursive Sequences with Linear Algebra

Keywords: CCS Concepts, Comput. method. → Symbolic calculus algorithms; Gröbner bases, P-recursive sequences, FGLM algorithm

Given several n-dimensional sequences, we first present an algorithm for computing the Gröbner basis of their module of linear recurrence relations.

A P-recursive sequence (u i ) i∈N n satisfies linear recurrence relations with polynomial coefficients in i, as defined by Stanley in 1980. Calling directly the aforementioned algorithm on the tuple of sequences (i j u i ) i∈N n j for retrieving the relations yields redundant relations. Since the module of relations of a P-recursive sequence also has an extra structure of a 0-dimensional right ideal of an Ore algebra, we design a more efficient algorithm that takes advantage of this extra structure for computing the relations.

Finally, we show how to incorporate Gröbner bases computations in an Ore algebra K t1, . . . , tn, x1, . . . , xn , with commutators

into the algorithm designed for P-recursive sequences. This allows us to compute faster the Gröbner basis of the ideal spanned by the first relations, such as in 2D/3D-space walks examples.

INTRODUCTION

Computing linear recurrence relations of multi-dimensional sequences is a fundamental problem in Computer Science. As this is the core of the paper, we start with an example. The binomial sequence b = ( i 1 ! /(i 2 ! (i 1 -i 2 )!)) (i 1 ,i 2 )∈N 2 satisfies Pascal's rule, a linear recurrence relation:

∀ (i1, i2) ∈ N 2 , bi 1 +1,i 2 +1 = bi 1 ,i 2 +1 + bi 1 ,i 2 . ( 1
)
It has constant coefficients. But this sequence also satisfies relations with polynomial coefficients in (i1, i2): Thanks to equations (2a) and (2b), we can compute any term of b starting with b0,0 = 1. This makes b P-recursive. Deciding if a sequence is P-recursive, as defined by Stanley [START_REF] Stanley | Differentiably finite power series[END_REF] finds applications for instance in Combinatorics and in Computer algebra [START_REF] Salvy | D-finiteness: Algorithms and applications[END_REF]. We refer the reader to the series of papers on planar [START_REF] Bostan | Non-D-finite excursions in the quarter plane[END_REF][START_REF] Bousquet-Mélou | Walks with small steps in the quarter plane[END_REF][START_REF] Bousquet-Mélou | Walks confined in a quadrant are not always d-finite[END_REF] and 3D-space walks [START_REF] Bostan | On 3-dimensional lattice walks confined to the positive octant[END_REF][START_REF] Bostan | Explicit formula for the generating series of diagonal 3D rook paths[END_REF].

∀ (i1, i2) ∈ N 2 , (i1 -i2 + 1) bi 1 +1,i 2 = (i1 + 1) bi 1 ,i 2 , ( 2a
) (i2 + 1) bi 1 ,i 2 +1 = (i1 -i2) bi
The generating series i 1 ,...,in∈N ui 1 ,...,in x i 1 1 • • • x in n of a P-recursive sequence (ui 1 ,...,in ) (i 1 ,...,in)∈N n is a D-finite function in x1, . . . , xn, see [START_REF] Lipshitz | D-finite power series[END_REF]Theorem 3.7]. A D-finite function is such that the dimension of the vector space spanned by all its derivatives is finite over K(x1, . . . , xn), this allows one to characterize such functions with initial conditions and differential equations. For instance, the Dynamical Dictionary of Mathematical Functions [START_REF] Benoit | The Dynamic Dictionary of Mathematical Functions (DDMF)[END_REF] generates dynamically and automatically web pages on elementary and special functions thanks to the equations they satisfy. Likewise, a P-recursive sequence is uniquely defined by the linear recurrence relations with polynomial coefficients it satisfies and a finite number of initial terms.

Related work

For one-dimensional sequences, C-relations (linear recurrence relations with constant coefficients) can be guessed with Berlekamp -Massey algorithm [START_REF] Berlekamp | Nonbinary BCH decoding[END_REF][START_REF] Massey | Shift-register synthesis and BCH decoding[END_REF], while for multidimensional sequences, they can be guessed with Berlekamp -Massey -Sakata algorithm [START_REF] Sakata | Finding a minimal set of linear recurring relations capable of generating a given finite two-dimensional array[END_REF][START_REF] Sakata | The BMS algorithm[END_REF] or Scalar-FGLM [START_REF] Berthomieu | Linear Algebra for Computing Gröbner Bases of Linear Recursive Multidimensional Sequences[END_REF][START_REF] Berthomieu | Linear Algebra for Computing Gröbner Bases of Linear Recursive Multidimensional Sequences[END_REF]. All these algorithms only guess and check that the relations are valid for a finite number of terms of the sequence. Then, usually, one need to prove that they are satisfied for all the terms of the sequence. In this paper, a table shall denote a finite subset of terms of a sequence: it is one of the input parameter of the algorithms since one cannot handle an infinite sequence.

Likewise, for one-dimensional sequences, P-relations (linear recurrence relations with polynomial coefficients) can be guessed with Beckermann -Labahn algorithm [START_REF] Beckermann | A uniform approach for the fast computation of matrix-type pade approximants[END_REF] for computing Hermite-Padé approximants.

When computing parameterized sums of multivariate sequences, or integrals of multivariate functions, in general, it is hopeless to obtain a closed form. Indeed, even for elementary univariate function we may not have closed forms when computing their integral on some domain. This motivated Zeilberger in 1990 [START_REF] Zeilberger | A holonomic systems approach to special functions identities[END_REF] to introduce the so-called creative telescoping method. It is now one of the most efficient ways to do so and has been a tremendously active research topic. We refer the reader to the up-to-date nice survey on the state of the art in creative telescoping in [START_REF] Chyzak | The ABC of Creative Telescoping -Algorithms, Bounds, Complexity[END_REF] and the references therein.

Although sequences can satisfy some linear recurrence re-lations with constant or polynomial coefficients, they still may be neither C-recursive nor P-recursive. Given a subset S of steps in {0, ±1} 2 , one can define the sequence u = (un,i,j) (n,i,j)∈N 3 of the number of planar walks of length n ending at (i, j) with steps in S. Determining if u is Precursive has applications in lattice path enumeration [START_REF] Banderier | Basic analytic combinatorics of directed lattice paths[END_REF][START_REF] Bousquet-Mélou | Walks with small steps in the quarter plane[END_REF][START_REF] Bousquet-Mélou | Walks confined in a quadrant are not always d-finite[END_REF]. Likewise, now 3D-space walks are investigated. In [START_REF] Bostan | On 3-dimensional lattice walks confined to the positive octant[END_REF], the authors use several different methods to prove that somes sequences of walks are P-recursive among which: the algebraic kernel method, the reduction to smaller dimensions and the guessing and proving method.

Contributions

We start by recalling the definitions of C-recursive (linear recursive with constant coefficients) sequences and define Crecursive tuples of sequences. This framework allows us to transform a P-recursive sequence into a tuple of sequences.

In Section 3, we extends Scalar-FGLM [START_REF] Berthomieu | Linear Algebra for Computing Gröbner Bases of Linear Recursive Multidimensional Sequences[END_REF][START_REF] Berthomieu | Linear Algebra for Computing Gröbner Bases of Linear Recursive Multidimensional Sequences[END_REF] for guessing the module of C-relations satisfied by a tuple of tables up to a degree bound d. The algorithm returns a (d+1)-truncated Gröbner basis of their module of relations. The main idea of this algorithm is the same as Scalar-FGLM's: extracting a full-rank submatrix of maximal rank.

P-relations with coefficients of degree at most δ satisfied by a sequence u = (ui 1 ,...,in ) (i 1 ,...,in)∈N n are C-relations satisfied by the tuple of sequences (i

j 1 1 • • • i jn n ui 1 ,...,in ) (i 1 ,...,in)∈N n , with j1, . . . , jn ≥ 0 and j1 + • • • + jn ≤ δ.
Hence, a first idea is to compute them with the extended Scalar-FGLM for instance. However, this strategy is not optimal because the sequences are not independent: the module of relations has an additional structure, see also [START_REF] Quadrat | A constructive study of the module structure of rings of partial differential operators[END_REF]. This is why, in Section 4, we design algorithm P-recursive-FGLM to take advantage of this particular structure. This algorithm input is a P-recursive table, a degree bound on the order of the relations and a degree bound on the polynomial coefficients. It returns a Gröbner basis of the right ideal of relations in a skew-polynomial ring K t1, . . . , tn, x1 . . . , xn with commutation rules x k x -x x k = t k t -t t k = t k xx t k = 0 for k = and t k x k -x k t k = x k . P-recursive-FGLM uses an extra criterion à la FGLM: if relation (2b) is only known to be valid for a given pair (i1, i2), then we still know that for any j1, j2 ≥ 0,

i j 1 1 i j 2 2 (i2 + 1) bi 1 ,i 2 +1 = i j 1 1 i j 2 2 (i1-i2) bi 1 ,i 2 .
This makes the extracted matrix smaller. Sometimes an estimate of the Gröbner basis staircase size is known. To optimize the number of table queries, we present an output sensitive probabilistic algorithm in Section 4.2.

The relations computed by P-recursive-FGLM can be output one after the other. When several relations are found, it seems profitable not to look for any relation in the ideal they span. In Section 5, we use a mixed approach by computing truncated Gröbner bases in the aforementioned skewpolynomial ring to determine "needed redundant relations". On the one hand, they are redundant for they are deduced from previously found relations. On the other hand, they are computed faster since this approach does not require any extra queries to the table. Furthermore, these relations can help us to close the staircase in order to guess prematurely the P-recursiveness of the table for a given degree. We illustrate this with Examples 6 and 7

In Section 6, we show that the matrices handled by Multi-Scalar-FGLM and P-recursive-FGLM are heavily structured. When the monomial ordering is LEX, we can make use of the quasi-Hankel structure to improve the complex-ity of the linear system solving [START_REF] Bostan | Solving Toeplitzand Vandermonde-like Linear Systems with Large Displacement Rank[END_REF]. We also investigate the number of table queries performed by the algorithms.

DEFINITIONS AND NOTATION

We recall the classical definitions of a C-recursive sequence and of a P-recursive sequence. We shall also see how to deal with C-relations between terms of multiple sequences and see that from a P-recursive sequence, we can create a tuple of sequences satisfying C-relations.

In all the paper, we shall use standard notation with bolded letters corresponding to vectors and sequences: for n ∈ N * , i = (i1, . . . , in) and x = (x1, . . . , xn). As usual,

x i denotes x i 1 1 • • • x in n and |i| = i1 + • • • + in. Finally, u = (u i ) i∈N n is an n-dimensional sequence over a field K. For i ∈ N n , we let [x i ]u = u i = [x i ],
if no ambiguity on u shall arise. We extend linearly this [ ] operator to K[x].

C-relations are straightforwardly extended to n-dimensional sequences [START_REF] Saints | Algebraic-geometric codes and multi-dimensional cyclic codes: Theory and algorithms for decoding using Gröbner bases[END_REF]Definition 21]: let S be a finite subset of N n , let α = (αs)s∈S be nonzero, then u = (u i ) i∈N n satisfies the linear relation defined by α if for all i ∈ N n , s∈S αs x i+s = s∈S αs u i+s = 0.

(3)

Operator [ ] allows us to make a one-to-one correspondence between C-relations satisfied by u and polynomials in K[x]: relation s∈S αs us = 0 is associated to polynomial s∈S αs x s . According to Sakata [START_REF] Sakata | Finding a minimal set of linear recurring relations capable of generating a given finite two-dimensional array[END_REF] and to Fitzpatrick and Norton [START_REF] Fitzpatrick | Finding a basis for the characteristic ideal of an n-dimensional linear recurring sequence[END_REF], a C-recursive sequence is defined as follows.

Definition 1. An n-dimensional sequence u = (u i ) i∈N n is C-recursive if from
a nonzero finite number of initial terms us, s ∈ S, and a finite number of C-relations, without any contradiction, one can compute any term of the sequence.

From the polynomial viewpoint, Definition 1 can be restated as sequence u is C-recursive if and only if set

I = {P ∈ K[x], ∀ i ∈ N n , P x i = 0} is a zero-dimensional ideal.
A Gröbner basis of this ideal can be computed using the so-called BMS [START_REF] Sakata | Extension of the Berlekamp-Massey algorithm to N Dimensions[END_REF] and Scalar-FGLM [5, Algorithm 3] algorithms, see also [START_REF] Berthomieu | Linear Algebra for Computing Gröbner Bases of Linear Recursive Multidimensional Sequences[END_REF] for an extended version of [START_REF] Berthomieu | Linear Algebra for Computing Gröbner Bases of Linear Recursive Multidimensional Sequences[END_REF].

C-recursive sequences are also exactly the only sequences such that their generating series i∈N n u i x i is a rational fraction in x whose denominator factors into univariate polynomials, see for instance [START_REF] Berthomieu | Linear Algebra for Computing Gröbner Bases of Linear Recursive Multidimensional Sequences[END_REF]Propostion 19] For a sequence u not C-recursive, both algorithms, called on the table of terms u i , |i| ≤ d, return a (d + 1)-truncated Gröbner basis. Sequences such as the binomial one b =

i 1 i 2 i∈N 2
may satisfy a C-relation but not be C-recursive.

For instance, b satisfies equation ( 1), Pascal's rule, while its generating series is the rational fraction

1 /(1 -x 1 -x 1 x 2 )
whose denominator is irreducible.

The greater class of P-recursive sequences satisfy P-relations (linear recurrence relations with polynomials coefficients). These relations are analogous to relations as equation ( 3) but with αs a polynomial in i. For instance, both equations (2a) and (2b) satisfied by b are P-relations. Definition 2. An n-dimensional sequence u = (u i ) i∈N n is P-recursive if from a nonzero finite number of initial terms u i , i ∈ S, and a finite number of P-relations with polynomials coefficients, without any contradiction, one can compute any term of the sequence.

Tuple of sequences

In this section, we aim to deal with multiple sequences u (1) , . . . , u (m) all n-dimensional over a common field K. In particular, we shall ask them to satisfy C-relations, of type

[P1] u (1) + • • • + [Pm] u (m) = 0.
We extend the operator [ ] to vector of polynomials: for all P = (P1, . . . , Pm)

∈ K[x] m , [P] = m =1 [P ] u ( )
. This yields a one-to-one correspondence between vectors of polynomials in K[x] m and linear relations. We now aim to find the K

[x]-submodule I = {(P1, . . . , Pm) ∈ K[x] m , ∀ i ∈ N n , [P1 x i ] u (1) + • • • + [Pm x i ] u (m) = 0} of K[x] m . We shall denote (e1, . . . , em) the canonical basis of K[x] m . Let u (1) = (3 i 1 ) i∈N 2 and u (2) = (3 i 1 + 2 i 2 ). The ideal of relations of u (1) is x2 -1, x1 -3 , while this of u (2) is x 2 2 -3 x2 + 2, x1 + 2 x2 -5 . Therefore, we need to know u (1) 0,0 , u (2) 0,0 , u (2)
0,1 to compute any term of both sequences. Yet, the module of relations satisfied by (u (1) , u (2) ) is spanned by {(x2 -1, 0), (x1 -3, 0), (1, x2 -2), (-2, x1 -1)} and we only need to know u

(1) 0,0 , u (2)
0,0 to compute any term.

P-recursive sequences

For P-recursive sequences, we can naturally create several auxiliary sequences such that they satisfy C-relations. For instance, with sequence b, we introduce the three auxiliary sequences: b

(0) = b, b (1) = i1 i 1 i 2 i∈N and b (2) = i2 i 1 i 2 i∈N
. Equation ( 1) is still satisfied by b (0) and both equations (2a) and (2b) can be restated as C-relations satisfied by these 3 auxiliary sequences:

∀ i ∈ N 2 , b (1)
i 1 +1,i 2 + b (2) i 1 +1,i 2 = b (1)
i 1 ,i 2 + b (0) i 1 ,i 2 , (4a) b (2)
i 1 ,i 2 +1 = b (1) i 1 ,i 2 -b (2)
i 1 ,i 2 . (4b)
Hence, their module of relations contains Span{(x1 3 . More generally, for any P-recursive sequence, one would have to introduce as many auxiliary sequences as needed depending on the number of indices and the allowed degree for the polynomial coefficients.

x2 -x2 - 1, 0, 0), (1, 1 -x1, x1), (0, -1, x2 + 1)} ⊆ K[x]
Usually a new set of variables t = (t1, . . . , tn) is introduced, such that [t k x i ] = i k u i , and thus [t j x i ] = i j u i . This yields commutation rules t k x k = x k (t k +1) and t k x = x t k for k = . We denote by K t, x the skew-polynomial ring with these commutation rules. We also let K(t) x be the ring of skew-polynomials in x whose coefficients are rational fractions in t.

In this setting, P-relations are in a one-to-one correspondence with polynomials in K t, x and a P-recursive sequence is a sequence such that the set of relations

I = {P ∈ K(t) x , ∀ i, j ∈ N n , P t j x i = 0} is a 0-dimensional right ideal, i.e. for all k, 1 ≤ k ≤ n, I ∩ K(t) x k = {0}. For b, this yields the right ideal spanned by (t1 -t2) x1 -(t1 + 1), t2 x2 + t2 -t1 ⊆ K(t) x .
In general, if one considers the right ideal I = {P ∈ K t, x , ∀ i, j ∈ N n , P t j x i = 0}, then for all k, there exists d > 0 and j such that t j x d k is the leading monomial of a polynomial in the ideal.

Lattice-based relations

We now deal with sequences satisfying C-relations when evaluated in points i ∈ Λ, a sublattice of N n .

The one-dimensional sequence u = (ui) i∈N defined by for all i ∈ N, u2i = u

(0) i = i! and u2i+1 = u (1) i = i! + 1 satisfies for all i ∈ N u2i+3-u2i+2-u2i+1+u2i = u (1) i+1 -u (0) i+1 -u (1) i +u (0) i = 0. (5)
but not ui+3 -ui+2 -ui+1 + ui = 0 for odd i. Nonetheless, computing equation 5 is of interest since with (u2i) i∈N and u1, it allows us to compute all the odd-indexed terms of u.

As in Section 2.1, we can write equation 5:

u (1) i+1 -u (0) i+1 -u (1) i + u (0) i = [(x + 1) (e2 -e1) x i ] = 0.
More generally, let (λ1, . . . , λn) be a basis of Λ, S be a finite subset of N n such that sequence u = (u i ) i∈N n satisfies

s∈S αs u i+s = 0, ∀ i ∈ Λ. ( 6
)
Then, one can create auxiliary sequences u (f ) = (u i+f ) i∈Λ for all f ∈ F , the fundamental domain of Λ to find the relation of equation 6.

Mimicking the P-recursive case, we introduce a new variable t such that for all i ∈ N n , t

x i = x i 1 λ 1 • • • x in λn t and [t x i ] = u i 1 λ 1 +•••+in λn , we can extend these relations to K t, x by linearity. For i ∈ Λ, we let i = c1 λ1 + • • • + cn λn, therefore x i = t x c 1 1 • • • x cn n .
Then, one obtains relation

s∈S αs x i+s = s∈S αs x s t x c 1 1 • • • x cn n = 0.

GUESSING RELATIONS FOR TUPLES OF SEQUENCES

This section is devoted to the design of an FGLM-like algorithm in the module K[x] m to compute the Gröbner basis of the submodule of relations of n-dimensional tables u (1) , . . . , u (m) defined over a common field K.

For the paper to be self-contained, we recall the notation and definitions of [START_REF] Berthomieu | Linear Algebra for Computing Gröbner Bases of Linear Recursive Multidimensional Sequences[END_REF]. We fix an admissible order ≺ on x, T is the ordered set of terms made from x. We extend ≺ to monomial vectors of type x i e , x i e with x i e

x i e if x i ≺ x i or if x i = x i and ≤ .

For P ∈ K[x] m , we denote by T (P) the set of terms appearing in P and by LT(P) the maximum of T (P) for ≺.

As we try to design an algorithm, we cannot expect to check that a relation is satisfied for all indices i ∈ N n . During the execution of the algorithm, we shall deal with a finite subset of indices U ⊂ N n , and we aim to find relations valid for those indices:

∀ i ∈ U , m =1 s∈S α ,s u ( ) i+s = 0. From the polynomial viewpoint, this translates to ∀ i ∈ U , m =1 s∈S α ,s x i+s e = 0.
Definition 3. Let U be a finite subset of N n . We say that a vector of polynomials P = (P1, . . . , Pm)

∈ K[x] m is valid up to U if [P x i ] = m =1 P x i u ( ) = 0 for all i ∈ U .

Staircases

Since we want to generate minimal relations, it is mandatory to determine a set of terms in T for which there is no C-relations. Useful staircases were defined for C-recursive tables in [5, Definition 6]. The following definitions extends it for multiple tables. Definition 4. Let T1, . . . , Tm, U be finite subsets of T such that T1, . . . , Tm ⊆ U and # U = m =1 # T . We say that a set S = S1 e1 ∪ • • • ∪ Sm em, with S1 ⊆ T1, . . . , Sm ⊆ Tm, is a useful staircase with V ⊆ U wrt. u (1) , . . . , u (m) , T = T1 e1 ∪ • • • ∪ Tm em, U , and

≺ if # V = # S and 1≤ ≤m s∈S β ,s [s v e ] = 0, ∀ v ∈ V
implies that β ,s = 0 for all , s, sets S1, . . . , Sm, V are maximal for the inclusion and minimal for ≺. We compare two ordered sets for ≺ by seeing them as tuples of their elements and then comparing them lexicographically.

As we are in a more general case, useful staircases still fail in general to be true staircases in K[x] m , see [START_REF] Berthomieu | Linear Algebra for Computing Gröbner Bases of Linear Recursive Multidimensional Sequences[END_REF]Example 3]. However, this problem is easily solved by stabilizing the staircase, i.e. adding to each set of the tuple all the divisors of its elements.

Algorithm for finding the relations

We design an algorithm for checking that a finite subset S ⊆ T is a useful staircase with V ⊆ U wrt. u (1) , . . . , u (m) , T, U and ≺. In the following example, we compute the useful staircase for two univariate tables.

Example 1. Let us consider u (1) = (3 i 1 ) i∈N 2 and u (2) = (3 i 1 +2 i 2 ) i∈N 2 . Let T1 = T2 = {1, x2} and U = {1, x2, x1, x 2
2 }, then we have the following potential relations

    e 1 e 2
x 2 e 1 x 2 e 2 1 1 2 1 3

x 2 1 3 1 5

x 1 3 4 3 5

x 2 2 1 5 1 9         α1 β1 αx 2 βx 2     = 0.
From the matrix viewpoint, a useful staircase is a set of independent columns with as many low-labeled columns as possible wrt. independent rows with as many low-labeled rows as possible. This matrix has rank 2: the first two columns (resp. the first two rows) are independent. Therefore, S = {e1, e2} is the useful staircase with V = {1, x2}. Thus, we can find relations (-α1 + x2) e1 -β1 e2 and -α 1 e1 + (-β 1 + x2) e2 with a linear system:

1 2 1 3 α1 α 1 β1 β 1 = 1 3 1 5 .
This yields (x2 -1) e1 = e1 + (x2 -2) e2 = 0.

We are now in a position to define the structured matrix associated to lists of terms T1, . . . , Tm ⊆ U ⊆ T .

Definition 5. Let T1, . . . , Tm, U be finite lists of terms such that T1, . . . , Tm ⊆ U ⊆ T and # U = m =1 # T . The multi-Hankel matrix HU,T has the same columns and rows as the concatenation of matrices HU,T 1 , . . . , HU,T m , where each HU,T is a multi-Hankel matrix for table u ( ) as defined in [START_REF] Berthomieu | Linear Algebra for Computing Gröbner Bases of Linear Recursive Multidimensional Sequences[END_REF]Definition 7]:

    ••• s∈ T 1 ••• ••• ••• s∈Tm ••• . . . . . . . . . . . . • • • . . . . . . . . . v∈U • • • [s v] u (1) • • • • • • • • • [s v] u (m) • • • . . . . . . . . . . . . • • • . . . . . . . . .     . If each S ⊆ T , V ⊆ U , # V = m =1 # S

and HV,S has full rank, then m

=1 S e is a useful staircase with V .

Furthermore, finding a monic polynomial vector P of a given support such that P x i = 0, for all x i ∈ U is equivalent to solving a linear system HU,T α+H U,LT(P) = 0, where T is the support of P without its leading term.

Proposition 1. Let T1, . . . , Tm, U be finite subsets of T . If the finite sets S1, . . . , Sm, V are such that S = m =1 S e is a useful staircase with V wrt. u (1) , . . . , u (m) , T, U and ≺, then det HV,S = 0, rank HV,S = rank HU,T .

If S 1 ⊆ S1, . . . , S m ⊆ Sm, V ⊆ V are finite sets of terms, then rank H V ,S ≤ rank HV,S.

Extension of Scalar-FGLM

As presented in [START_REF] Berthomieu | Linear Algebra for Computing Gröbner Bases of Linear Recursive Multidimensional Sequences[END_REF][START_REF] Berthomieu | Linear Algebra for Computing Gröbner Bases of Linear Recursive Multidimensional Sequences[END_REF], Scalar-FGLM is an algorithm whose input is a table u, an integer d and a monomial order ≺ refined by the total degree. If u is C-recursive of order less than d, then the algorithm returns a Gröbner basis of the ideal of relations of u. Otherwise it returns a (d + 1)truncated Gröbner basis of an ideal such that the relations induced by the Gröbner basis are only valid up to degree d. This extension is natural: the input shall be a vector of n-dimensional tables (u (1) , . . . , u (m) ), a vector of integers (d1, . . . , dm) and a monomial order ≺ on K[x] which we extend to K[x] m as specified earlier.

Algorithm 1, Multi-Scalar-FGLM, computes the module of relations of a list of tables u (1) , . . . , u (m) up to degrees d1, . . . , dm.

Algorithm 1: Multi-Scalar-FGLM. Return G.

Example 2. Let us unroll the algorithm on a modified version of Example 1, with tables u (1)

= (3 i 1 ) i∈N 2 , u (2) = (3 i 1 + 2 i 2 ) i∈N 2 and (d1, d2) = (1, 1). First, T = {1, x2, x1} e1 ∪ {1, x2, x1} e2, so that U = {1, x2, x1, x 2 2 , x1 x2, x 2 1 } and matrix HU,T is        e 1 e 2 x 2 e 1 x 2 e 2 x 1 e 1 x 1 e 2 1
1 2 1 3 3 4

x 2 1 3 1 5 3 5

x 1 The first two independent columns are e1, e2 while the first two independent rows are 1, x2. This means that the useful staircase is S = {e1, e2} with V = {1, x2}. The true staircase is then S = S and we initialize L with {x2, x1, x 2 2 , x1 x2, x 2 1 } e1 ∪ {x2, x1, x 2 2 , x1 x2, x 2 1 } e2. For = 1, we set τ = x2, solve the system HV,S α + ( 11 ) yielding P = (x2 -1) e1 and update L to {x1,

x 2 1 } e1 ∪ {x2, x1, x 2 2 , x1 x2, x 2 1 } e2.
Then, we set τ = x1, solve the system HV,S α + ( 33 ) yielding P = (x1-3) e1 and update L to {x2, x1, x 2 2 , x1 x2, x 2 1 } e2. For = 2, we set τ = x2, solve the system HV,S α + ( 35 ) yielding P = e1 + (x2 -2) e2 and update L to {x1, x 2 1 } e2. Then, we set τ = x1, solve the system HV,S α + ( 45 ) yielding P = -2 e1 + (x1 -1) e2 and update L to ∅.

The algorithm returns {P, P , P , P }.

GUESSING RELATIONS FOR P-RECURSIVE SEQUENCES

As stated in Section 2.2, one could use Algorithm 1 to deal with P-recursive tables. For instance, if one suspects a sequence u to satisfy relations of type s∈S αs(i + s) u i+s = 0 with αs ∈ K[i] of degree at most δ, then one would just need to call the algorithm with all the u (j) = (i j u i ) i∈N , s.t. |j| ≤ δ, as input tables.

A criterion for P-recursive sequences

The main drawback of this strategy is the output of the algorithm: there are polynomials representing redundant relations if the sequence is P-recursive but not C-recursive. For instance, Algorithm 1 called on b (0) , b (1) , b (2) , of Section 2.2, with degree bounds (1, 1, 1) in x returns both relations (4a) and (4b) together with six extra relations.

Yet, the first two are enough to compute every term of the sequence starting from 0 0 = 1. These superfluous relations come from the fact that we forget the extra structure of the module of relations, i.e. it is also a right ideal of K t, x for t = (t1, . . . , tn), t k x = x t k if k = and t k x k = x k (t k + 1).

Therefore, it is possible to extend Multi-Scalar-FGLM to take into account this structure. Whenever a relation t j+δ x i+d = r∈R,s∈S αr,s t j+r x i+s (i + d) j+δ u i+d = r∈R,s∈S αr,s (i + s) j+r u i+s is true up to a certain threshold for i, then its multiples by i k + d k are also true up to the same threshold.

Proposition 2. Let u be a P-recursive table over K. Let T (resp. U ) be a finite subset of monomials in K t, x (resp. in K[x]). If t j x i is not in the useful staircase of matrix t j x i+i t j x i ∈T,x i ∈U then neither is any multiple t j+j x i .

Proof. If t j x i is not in the useful staircase, a relation starting with u i+k exists. Multiplying it by i k + d k yields a new relation whose associated polynomial has a leading term t k t j x i , hence column t k t j x i can also be removed.

Let us notice that this is analogous to the behavior of FGLM for polynomial ideals in K[x] [START_REF] Faugère | Efficient Computation of Zero-dimensional Gröbner Bases by Change of Ordering[END_REF]: whenever a relation x d = s≺d αs x s is found, then automatically for all i ∈ N n , x i+d = s≺d αs x i+s . This is also analogous to Buchberger -Möller algorithms wherein the finitely many points can be interpreted as the finitely many table terms we have. We illustrate this criterion.

Example 3. We consider u = (i!) i∈N and skew-polynomial ring K t, x . Columns will be labeled by T = {1, t, x, t x} and rows by U = {1, x, x 2 , x 3 }. This yield the rank 3 matrix

HU,T =     1 t x t x 1 1 0 1 1 x 1 1 2 4
x 2 2 4 6 18

x 3 6 18 24 96     .
Column labeled x is the first dependent one meaning we find relations [x i+1 -(t + 1)

x i ] = ui+1 -(i + 1) ui = 0 for all i ∈ {0, . . . , 3}. This leads us to determine that (i + 1) ui+1

+ (i 2 + 2 i + 1) ui = 0 = [t x i+1 + (t 2 + 2 t + 1) x i ] for all i ∈ {0, . . . , 3} without computing [t 2 ], [t 2 x], [t 2 x 2 ], [t 2 x 3 ],
i.e. column t 2 . Therefore, column t x is independent from the previous ones but is not in the useful staircase.

In this context, a useful staircase is the maximal subset S ⊆ T such that columns of HU,S are all linearly independent and for all t j x i ∈ T \ S and j ∈ N n , t j+j x i ∈ S.

To compare monomials in K t, x in order to build our matrix, we proceed as follows. We let ≺ be a monomial ordering refined by the total degree on t and x with the convention that for any k, , t k ≺ x .

The following algorithm, P-recursive-FGLM, computes the linear recurrence relations with polynomial coefficients of an input P-recursive table. Thanks to the criterion given in Proposition 2, the linear systems that we solve at the end are smaller than in Algorithm 1.

Algorithm 2: P-recursive-FGLM.

Input: A table u = (u i ) i∈N n with coefficients in K, (δ, d)
given bounds and ≺ a monomial ordering. Output: A reduced (δ, d)-truncated Gröbner basis wrt. ≺ of the ideal of relations of u in K t, x . Let T be the ordered list of monomials in t, x of bidegree at most (δ, d). Let U be the ordered list of the first # T monomials in x. Build the matrix H U,T . Find S the useful staircase. Find V the set of first linearly independent rows of H U,S . S := Stabilize(S). L := T \ S .

// list of next terms to study G := { }.

// the future Gröbner basis

While L = ∅ do τ := min≺(L). Find α s.t. H V,S α + H V,{τ } = 0. G := G ∪ τ +
V = {1, x2, x1, x 2 2 , x1 x2}.              1 t 2 t 1 x 2 x 1 t 2 x 2 t 2 x 1 t 1 x 2 t 1 x 1 1 1 0 0 1 1 1 0 0 1 x 2 1 1 0 1 /2 1 1 1 0 1 x 1 1 0 1 1 2 1 0 1 4 x 2 2 1 /2 1 0 1 /6 1 /2 1 /2 1 0 1 /2 x 1 x 2 1 1 1 1 /
             (7)
We set L = {x1, t2 x2, t2 x1, t1 x1} and let τ = x1. We then find G = {x1 -t1 -1} and L is updated to {t2 x2}.

We set τ = t2 x2, find G = {x1 -t1 -1, t2 x2 -1} and update L to ∅. Remark 3. It is sometimes more convenient to obtain the ideal of relations, not in K t, x but in K(t) x , the ring of skew-polynomials in x with rational fractions in t coefficients. In this setting, we can mimic Scalar-FGLM algorithm by setting L to the set of all monomials in t, x of bidegree at most (δ, d + 1) which are not in S . In a way, this would ensure that the staircase of the Gröbner basis obtained at the end is closed. Then, we need to change also how to update L: since now each t k is invertible, whenever τ = t j x i , we remove every monomial t j x i+i .

Adaptive algorithm

If one is concerned about the number of queries to table u, it is necessary not to visit monomials yielding relations which are multiple of already computed relations. Therefore, we need to have an adaptive approach of the algorithm: we shall discover the staircase term by term, see also Adaptive Scalar-FGLM [5, Section 5]. As in Remark 3, we assume that we look for a Gröbner basis of the ideal in K(t) x .

We start from S = ∅, V = ∅ and try to increase them both: if H V ∪{x i },S∪{t j x i } has a greater rank than HV,S, then S and V are updated accordingly. Otherwise, we have found a new relation whose leading term is t j x i . Therefore, we do not update S nor V and shall never suggest any multiple of x i (resp. x i ) as a new term in S (resp. V ). This yields Algorithm 3.

Input: Tables u ( 1 )

 1 , . . . , u (m) over K of dimension n, (d 1 , . . . , dm) given bounds and ≺ a monomial order. Output: A reduced (d 1 + 1, . . . , dm + 1)-truncated Gröbner basis wrt. ≺ of their module of relations. Let T = T d 1 e 1 , ∪ • • • ∪ T dm em where T d is the ordered list of monomials of degree at most d. Let U be the ordered list of the first # T monomials in x. Build the matrix H U,T . Extract a submatrix of maximal rank by finding S = S 1 e 1 ∪ • • • ∪ Sm em and V such that S is the useful staircase for V , i.e. s.t. rank H U,T = rank H V,S . S := Stabilize(S). L := {T d 1 +1 \ S 1 , . . . , T dm+1 \ S m }. // list of next terms G := { }. // the future Gröbner basis For from 1 to m do While L[ ] = ∅ do τ := min≺(L[ ]) and remove τ from L[ ]. Find α s.t. H V,S α + H V,{τ e } = 0. G := G ∪ τ e + s∈S αs s . Sort L[ ] by increasing order (wrt. ≺) and remove multiples of LT(G).

Example 5 .

 5 Algorithm 2 on the binomial table b returns the three relations (1), (2a) and (2b) as {t2 x2-(t1-t2), (t1-t2) x1 -(t1 + 1), x1 x2 -x2 -1} for degree bounds at least (1, 2) in t and x and ordering DRL(t2 ≺ t1 ≺ x2 ≺ x1). Therefore during the execution of the algorithm, the least monomials of list L were t2 x2, t1 x1 and x1 x2.Over K(t) x , we would have to consider that x1 x2 is a multiple of t2 x2, when updating L. Therefore, Pascal's rule is never computed. This is consistent with the fact that g1 = t2 x2 -t1 + t2 and g2 = (t1 -t2) x1 -(t1 + 1) form a Gröbner basis in K t, x , that (x1 x2 -x2 -1) ∈ I = g1, g2 but that g1 + g2 x2 = (t1 -t2)(x1 x2 -x2 -1) ∈ I. Hence, x1 x2 -x2 -1 is a member of the ideal I K(t) x .

  Example 4. Let us unroll Algorithm 2 on table ( i 1 ! /i 2 !) i∈N 2 with δ = d = 1 and ordering DRL(t2 ≺ t1 ≺ x2 ≺ x1). As x1 is not in this set, we need to remove both t2 x1 and t1 x1 by Proposition 2. Hence the useful staircase is S = {1, t2, t1, x2, t1 x2}. The set of the first 5 independent rows is

	s∈S αs s . Sort L by increasing order (wrt. ≺) and remove
	multiples of LT(G).
	Return G.
	We have T = {1, t2, t1, x2, x1, t2 x2, t2 x1, t1 x2, t1 x1} and
	we let U = {1, x2, x1, x 2 2 , x1 x2, x 2 1 , x 3 2 , x1 x 2 2 , x 2 1 x2}. This
	yields matrix HU,T described in equation 7.
	The first linearly independent columns are 1, t2, t1, x2, t2 x1,
	t1 x2, t1 x1.

  /6 1 /2 0 1 /24 1 /6 1 /6

					2	2	1	2	1 /2	4
	x 2 1	2	0	4	2	6	2	0	4	18
	x 3 2	1 1 /2	0	1 /6
	x 1 x 2 2	1 /2 1 1 /2 1 /6	1	1 /2	2	1 /6	2
	x 2 1 x 2	2	2	4	1	6	2	6	2	18

Acknowledgments

We thank the anonymous referees for their careful reading and their helpful comments, C. Banderier, A. Bostan, F. Chyzak and B. Salvy for valuable discussions. This work has been partly supported by the French National Research Agency ANR-11-BS02-0013 HPAC project.

MIXED APPROACH WITH GRÖBNER BASES COMPUTATIONS

When dealing with a sequence u = (u i ) i∈N n , not all its ideal of relations I in K t, x is necessarily interesting. For instance, whenever one wants to know if u is P-recursive, it is enough to compute for all k, 1 ≤ k ≤ n, a polynomial Algorithm 3: Adaptive P-recursive-FGLM.

Input: A table u = (u i ) i∈N n with coefficients in K, (δ, d) given bounds and ≺ a monomial ordering. Output: A reduced Gröbner basis of a zero-dimensional ideal of degree ≥ d. L := {1}, M := {1}.

// list of next terms to study S := { }, V := { }.

// S is a useful staircase for V G := { }.

// leading terms of the Gröbner basis H := { }.

// deleted rows While L = ∅ and M = ∅ do τ := min≺ L, χ := min≺ M . If H V ∪{χ},S∪{τ } is full rank then S := S ∪ {τ }, V := V ∪ {χ}, L := L ∪ {t 1 τ, . . . , tn τ, x 1 τ, . . . , xn τ | deg t t i τ ≤ δ} M := M ∪ {x 1 χ, . . . , xn χ} Sort L (resp. M ) by increasing order (wrt. ≺) and remove duplicates and multiples of G (resp. H ).

with T the set of monomials in t, x of bidegree at most (deg t τ,

Return S, V and G.

Else G := G ∪ {τ }, H := H ∪ {χ} and remove multiples of τ in L and of χ in M .

Error "Run Algorithm 2".

in I whose leading monomial is t j x d k with j any and d > 0. We recall that this ensures that I K(t) x is 0-dimensional. Therefore, only a subideal of I may be needed to certify that u is P-recursive.

In Adaptive P-recursive-FGLM algorithm, as soon as a polynomial P1 with leading monomial LT(P1) is found, all the t j x i LT(P1) are discarded so that the multiples of P1 are not computed. Therefore, whenever polynomials P2, . . . , Pr are found, we can discard all the right-linear combinations of P1, . . . , Pr, hence the ideal they span. In this algorithm, if P1, . . . , Pr are exactly the first output polynomials, then we know that they are in the final Gröbner basis and that their leading monomials are the least first. By computing the Gröbner basis of J = P1, . . . , Pr ⊆ K t, x , we can expect to be able to find polynomials with leading monomials not in LT(P1), . . . , LT(Pr) . On the one hand, these polynomials are redundant since they are deduced from P1, . . . , Pr, i.e. the relations they represent are deduced by those induced by P1, . . . , Pr. On the other hand, we may find polynomials ensuring that J is zero-dimensional and thus that I is.

Let us assume that when computing the Gröbner basis of J, we find another polynomial Pr+1. Let us assume furthermore that, seen as monomials in x, LT(P1), . . . , LT(Pr) ≺ LT(Pr+1) and that Pr+1 is also in the Gröbner basis of I. Then, we computed Pr+1 without new queries to the table.

Remark 4. The Gröbner basis computation can be done either in K t, x or in K(t) x . In the latter case, as soon as a relation with leading monomial x i is found, no relation with leading monomial x i+i shall be computed. However, in the former case, it is possible that we find new relations with leading monomials x i+i , hence with higher degrees in x, but with also lower degrees in t.

In the black-box model, this approach can be very beneficial, since the computation of even one element of u can be costly. Therefore, we can modify Adaptive P-recursive-FGLM so that it takes these Gröbner bases computations into account: first we compute the polynomials as soon as they are found, instead of just storing their leading monomials in G . Then, whenever # G is big enough, we replace G by the true Gröbner basis of the ideal spanned by G and we remove all the multiples of LT(g), for g ∈ G, in L.

Computation of Gröbner bases in an Ore algebra can be done with the adaptations of Buchberger and F4 [START_REF] Faugère | A new efficient algorithm for computing Gröbner bases (F4)[END_REF] algorithms in computer algebra systems [START_REF] Levandovskyy | Non-commutative Computer Algebra for polynomial algebras: Gröbner bases, applications and implementation[END_REF], such as Maple or Plural the noncommutative extension of Singular.

Let us illustrate this behavior with an example.

Example 6. We consider Gessel's walks and sequence g = (gn,i,j) (n,i,j)∈N 3 denoting the number of walks of length n ending at (i, j), with steps in {(-1, 0), (-1, -1), (1, 0), (1, 1)}. The kernel method, see [START_REF] Bostan | On 3-dimensional lattice walks confined to the positive octant[END_REF] for instance, ensures that g is Precursive if and only if so are both g = (gn,0,j) (n,j)∈N 2 and g = (gn,i,0

We apply this mixed approach to find a 0-dimensional subideal of relations of table g up to degrees [START_REF] Berthomieu | Linear Algebra for Computing Gröbner Bases of Linear Recursive Multidimensional Sequences[END_REF][START_REF] Berthomieu | Linear Algebra for Computing Gröbner Bases of Linear Recursive Multidimensional Sequences[END_REF]. Variables t0, x0 shall be associated with n, while t2, x2 shall be associated with j. The first computed polynomial has leading monomial t0 t Hence we can suspect that g is P-recursive.

Given a 3D-space walk sequence w = (w n,i,j,k ) (n,i,j,k)∈N 3 counting the number of walks of length n ending at (i, j, k), the kernel method still applies. Thanks to the C-relation satisfied by said walk, one can prove, see [START_REF] Bostan | On 3-dimensional lattice walks confined to the positive octant[END_REF]Section 4.3] that the generating series W (x) = n,i,j,k≥0 w n,i,j,k x n 0 x i 1 x j 2 x k 3 is a K(x)-linear combination of W1 = W (x0, 0, x2, x3), W2, W3, W12 = W (x0, 0, 0, x3), W13, W23, W123 = W (x0, 0, 0, 0) and 1. That is, there is a C-relation between w, (w n,0,j,k ) (n,j,k)∈N 3 , . . . , (w n,0,0,k ) (n,k)∈N 2 , . . . , (wn,0,0,0) n∈N .

In the next example, we deal with the 3D-space walk w whose steps are in {(-1, -1, -1), (-1, -1, 1), (-1, -1, 0), (1, 0, 0)}, see [START_REF] Bostan | On 3-dimensional lattice walks confined to the positive octant[END_REF]Proposition 6].

Example 7. We use the same notation as in Example 6 with the extra convention that t1, x1 (resp. t3, x3) are associated with i (resp. k).

Running Multi-Scalar-FGLM on w and all seven tables obtained from it when possibly setting parameters i, j, k to 0 yields the kernel equation

Increasing gradually to column degree (2, 4) and row degree 30, we obtain a 1448 × 350 matrix and 8 leading terms. A Gröbner basis computation yields 25 leading terms including t0 t1 t2 x 2 1 . This closes the staircase on the x1 side. We increase now the column degree to [START_REF] Benoit | The Dynamic Dictionary of Mathematical Functions (DDMF)[END_REF][START_REF] Berlekamp | Nonbinary BCH decoding[END_REF] and the row degree to 40, obtain a 2959 × 600 matrix and find an extra 25 relations. Unfortunately this does not close the staircase on either side.

A Gröbner basis computation, ideally truncated in total degree 9, returns 127 relations. Among which, there are leading monomials t 3 0 t2 t 2 3 x 4 0 and t0 t 2 1 t2 x 4 2 . In total, 19 monomials are purely in t and one of the xi. Hence, the staircase is closed in K(t) x .

In Table 1, we detail the linear algebra part for W3. These computations have been integated in and done with FGb [START_REF] Faugère | FGb: A Library for Computing Gröbner Bases[END_REF].

Degree

COMPLEXITY ANALYSIS

This section is devoted to the study of the complexity of both Multi-Scalar-FGLM and P-recursive-FGLM algorithms. This comes down to the study of the structure of the multi-Hankel matrix.

Algorithm for tuples of tables

For the Multi-Scalar-FGLM algorithm, the multi-Hankel matrix HU,T is the concatenation of m multi-Hankel matrices HU,T . From [6, Section 6], we know that for an ndimensional table and LEX(x1 ≺ • • • ≺ xn) ordering, if the leading monomials are x d 1 1 , . . . , x d 2 2 , . . . , x dn n , then the multi-Hankel matrix created by Scalar-FGLM is a submatrix of a multilevel block Hankel matrix of depth n, or n-multiblock Hankel matrix which is heavily structured. Definition 6. A multilevel block Hankel matrix of depth 0 is a scalar. For any n ∈ N, a multilevel block Hankel matrix of depth n + 1 is a block Hankel matrix whose blocks are multilevel block Hankel matrices of depth n.

By construction in Algorithm 1, each T of T = T1 e1 ∪ • • • ∪ Tn en is a subset of U , therefore HT ,T e is a submatrix of the generalized multi-Hankel matrix handled by the algorithm. Let S = S1 e1 ∪ • • • ∪ Sm em and V be such that S is a useful staircase with V . A necessary condition for a monomial s to be in S is that column s e is linearly independent from the previous ones in HT ,T e . Thus, row s e is also independent from the previous ones in HU,T . Hence, s ∈ S implies that s ∈ V and each S is a subset of V .

Let Σ = m =1 S , then Σ ⊆ V . Let M be the least integer such that V ⊆ W = Σ ∪ M k=2 w k Σ. Now, HV,S is a submatrix of HW,Σ e 1 ,∪•••∪Σ em and HW,Σ e 1 ,∪•••∪Σ em is a block matrix whose blocks Hw k Σ,Σ are all n-multiblock:

Proposition 5. Given a tuple of m tables. Let µ be the number of polynomials in the Gröbner basis output by Multi-Scalar-FGLM. For 1 ≤ k ≤ n, let d k be the greatest powers in x k appearing in the Gröbner basis. Let M be the number of blocks of size d1 • • • dn appearing on each column of the computed matrix by the algorithm. Then the number of operations in the base field to compute the Gröbner basis is no more than

Proof. Applying Hankel displacement operator makes each block Hw k Σ,Σ have displacement rank d2 • • • dn, see [START_REF] Berthomieu | Linear Algebra for Computing Gröbner Bases of Linear Recursive Multidimensional Sequences[END_REF]Proposition 17]. Therefore, HW,Σ e 1 ∪•••∪Σ em have displacement rank M m d2 • • • dn. Thanks to [START_REF] Bostan | Solving Toeplitzand Vandermonde-like Linear Systems with Large Displacement Rank[END_REF], solving a linear system with this matrix can be done in

. Finally, we solve µ systems.

Algorithms for P-recursive tables

For j ∈ N n , let Σ j be the set of all monomials t j x i in the useful staircase. By the column removal criterion of Proposition 2, we know that Σ j+j ⊆ Σ j for any j ∈ N n . Hence for all j, Σ j ⊆ Σ0 = Σ.

The remaining part of the analysis is the same as for Multi-Scalar-FGLM called on m = n+δ δ = O(n δ ) tables. This yields the following corollary. Corollary 6. Let δ be the maximal degree in t allowed to compute the relations of a table. Let µ be the number of polynomials in the Gröbner basis output by P-recursive-FGLM. For 1 ≤ k ≤ n, let d k be the greatest powers in x k appearing in the Gröbner basis. Let M be the number of blocks of size d1 • • • dn appearing on each column of the computed matrix by the algorithm. Then the number of operations in the base field to compute the Gröbner basis is no more than

Concerning the number of table queries, Algorithm 2, Precursive-FGLM, can reach elements of u far away from the useful staircase whenever the degree bound in x is too big. In Algorithm 3, Adaptive P-recursive-FGLM, the matrix elements are [t j x i+i ] = (i + i ) j u i+i with t j x i ∈ S, x i ∈ Σ and x i ∈ V . Hence, the number of table queries is the cardinal of Σ V = {τ χ | τ ∈ Σ, χ ∈ V }.

Embedding V into Σ ∪ M k=2 w k Σ ensures that Σ V ⊆ (2 Σ) ∪ M k=2 w k (2 Σ), where 2 Σ = {τ τ | τ, τ ∈ Σ}. Hence, # (Σ V ) ≤ M # (2 Σ) ≤ M (# Σ) 2 /2 in the worst case.

Whenever Σ is included in a parallelotope of size C # Σ, then Ruzsa's Theorem [START_REF] Ruzsa | Generalized arithmetical progressions and sumsets[END_REF]Theorem 1.1] ensures that there exists c such that # (2 Σ) ≤ c # Σ, hence # (Σ V ) ≤ c M # Σ.


est destinée au dépôt

Seismic zonation of Hispaniola Island for hazard assessment purposes

Monique Terrier, Didier Bertil, Myriam Belvaux

To cite this version:

___________________________________________________________________________

The island of Hispaniola, which includes Haiti and the Dominican Republic, is a strongly active seismic zone at the level of the island arc of the Greater Antilles, on the boundary of the North American and Caribbean plates. The geodetic data and the geological reconstructions indicate to this place a convergence of plates of the order of 2 cm per year on average. The earthquake on January 10, 2010, which has destroyed the capital of Haiti, Port-au-Prince, came to add to the long list of destructive earthquakes which hit Hispaniola during the last centuries. The seismic hazard is however marred strong uncertainties due to a complex seismotectonic environment and insufficient knowledge of the observed seismicity. This seismic zoning is the fruit of a sismotectonique synthesis realized within the framework of assessment of hazard and seismic risk in the capital of the Dominican Republic, Santo Domingo (260 km2; 3 million people) (project PNUD). Initially, a unified catalog of earthquakes was established, collecting historical and instrumental data. The oldest known earthquake dated 1562. More than 25 epicenters of earthquakes of magnitude Mw at least equal to 6,5 are located at the level of Hispaniola or nearby immediate, among whom the earthquakes of magnitude superior to 7,5, of 1691, on 1842, on 1887 and 1946. To best clarify the focal characteristics of historical earthquakes, the study sought to find the most contemporary documentary sources describing the seismic events.

Although there are uncertainties on the magnitude and location, the catalogue of seismicity is considered complete from since 1945 for above 5.0 magnitude.

Figure 1: Seismic zonation of Hispaniola Island (Caraïbes)

Mapping active faults of Hispaniola was carried out on the basis of scientific publications, digital models of terrain, and geological maps (including maps to 1: 50 000 available on the Dominican Republic). The information on the geometry of the faults, the kinematics of their current movement, and the characteristics of their geologic evolution is integrated into the GIS. Based on this sismotectonique interpretation, we distinguish 3 types of seismic sources: intercrustaux domains, major active faults, and the North American subduction zone. Each seismic source is described according to the mechanism and current rate of deformation, the associated seismicity, the estimated maximal magnitude and the corresponding return period. The parameters of seismic activity of the sources used for probabilistic seismic hazard calculations are inferred mainly from these characteristics. This zoning bases itself on a state of the knowledge of the active faults or potentially active of this part of the Caribbean margin in the form of synthetic seismotectonic covering the whole of the island of Hispaniola and its margins. It is applied to the evaluation of the seismic hazard of the urban agglomeration of Santo Domingo, whose results will be the subject of a later publication. This study demonstrates that a model of seismic sources based on the most complete possible description of the potentially active faults is essential in the regions where the lack of data of seismicity impacted strongly seismic hazard calculations.

zonation of Hispaniola Island for hazard assessment purposes

Monique Terrier, Didier Bertil, Myriam Belvaux. Seismic zonation of Hispaniola Island for hazard assessment purposes. 35th International Geological Congress : IGC 2016, Aug 2016, Cape Town, South Africa. hal-01314288 th International Geological Congress -27 august -4 September 2016 -Cape Town, South Africa Seismic Terrier, M., Bertil, D. and Belvaux, M. BRGM, Orléans, France


Study of the

___________________________________________________________________________

This study is part of the Tandem project (Tsunamis in the Atlantic and the English ChaNnel: Definition of the Effects through numerical Modeling, 2014-2018), dedicated to the appraisal of coastal effects due to tsunami waves on the French coastlines. Within this framework, the present abstract concerns 1) the inventory of the possible seismic sources located between the Celtic Sea and the Cadix gulf and 2) the proposal of the worst-case source scenarios.

The study gathered and analysed a large number of information sources (databases, publications, maps) necessary for the realization of a seismotectonic synthesis of this part of the western area of the European plate. The map of active or potentially active faults is accompanied by a table of information indicating: the geometry of the fault system (segmentation, length, dip), its geological features, earthquakes attached to the structure, the current kinematics (movement type, strain rate), the main bibliographical references.

Figure 1 -Selected seismic sources (line in purple) for the assessment of the tsunami susceptibility of the French Atlantic coast and English channel

From this seismotectonic synthesis, several scenarios for the simulation of tsunami have been chosen according to the size of the seismic sources (sufficient to generate at least 6magnitude earthquakes), their location (at sea), and their current activity. The scenarios of seismic rupture proposed are, a priori, the most penalizing for the French Coast along the Atlantic facade and English channel.

Each zone is characterized by: 1) the nature of the movement, 2) the major earthquakes associated with the structure, 3) the rate of deformation, 4) the maximun possible magnitude estimated on the basis of the paleoseismic data and/or on the geometry of the fault, 5) the geometry of the coseismic rupture for the estimated maximum magnitude (and corresponding uncertainties).

The first modelling results are presented. Upcoming results for this part of the Tandem project will concern the interpretation of the results in terms of susceptibility of the French coast to tsunamis generated by seismic ruptures.

  possible tsunamigenic seismic sources for the French coast (Atlantic facade and English Channel) Monique Terrier, Anne Lemoine To cite this version: Monique Terrier, Anne Lemoine. Study of the possible tsunamigenic seismic sources for the French coast (Atlantic facade and English Channel). 35th International Geological Congress : IGC 2016, Aug 2016, Cape Town, South Africa. hal-01314299 Study of the possible tsunamigenic seismic sources for the French coast (Atlantic facade and English Channel) Terrier, M., Lemoine, A., Pedreros, R. BRGM, Orléans, France


Vedran Vukotić
email: vedran.vukotic@irisa.fr

Christian Raymond
email: christian.raymond@irisa.fr

Guillaume Gravier
email: guillaume.gravier@irisa.fr

Bidirectional Joint Representation Learning with Symmetrical Deep Neural Networks for Multimodal and Crossmodal Applications

Keywords: neural networks, deep learning, representation, embedding, multimodal, crossmodal, retrieval, video retrieval, video hyperlinking, image and text, autoencoder, bidirectional learning, tied weights, shared weights

Common approaches to problems involving multiple modalities (classification, retrieval, hyperlinking, etc.) are early fusion of the initial modalities and crossmodal translation from one modality to the other. Recently, deep neural networks, especially deep autoencoders, have proven promising both for crossmodal translation and for early fusion via multimodal embedding. In this work, we propose a flexible crossmodal deep neural network architecture for multimodal and crossmodal representation. By tying the weights of two deep neural networks, symmetry is enforced in central hidden layers thus yielding a multimodal representation space common to the two original representation spaces. The proposed architecture is evaluated in multimodal query expansion and multimodal retrieval tasks within the context of video hyperlinking. Our method demonstrates improved crossmodal translation capabilities and produces a multimodal embedding that significantly outperforms multimodal embeddings obtained by deep autoencoders, resulting in an absolute increase of 14.14 in precision at 10 on a video hyperlinking task (α = 10 -4 ).

INTRODUCTION

Deep neural networks have been long know to produce meaningful data representations [START_REF] Hinton | Reducing the dimensionality of data with neural networks[END_REF], either as deep belief networks, autoencoders or a combination of both. More recently, deep neural networks have been successfully deployed in tasks requiring consideration of multiple modalities. These tasks vary from retrieval [START_REF] Feng | Cross-modal retrieval with correspondence autoencoder[END_REF][START_REF] Srivastava | Learning representations for multimodal data with deep belief nets[END_REF][START_REF] Lu | Semantic retrieval of personal photos using a deep autoencoder fusing visual features with speech annotations represented as word/paragraph vectors[END_REF], ranking [START_REF] Weston | Large scale image annotation: learning to rank with joint word-image embeddings[END_REF] and classification tasks [START_REF] Cha | Multimodal sparse representation learning and applications[END_REF][START_REF] Ngiam | Multimodal deep learning[END_REF] to generative tasks [START_REF] Srivastava | Learning representations for multimodal data with deep belief nets[END_REF]. Data ICMR 2016, June 6, 2016, New York, USA. c 2016 Copyright held by the owner/author(s) often consist of bimodal pairs such as images and tags [START_REF] Cha | Multimodal sparse representation learning and applications[END_REF], images and speech [START_REF] Feng | Cross-modal retrieval with correspondence autoencoder[END_REF][START_REF] Lu | Semantic retrieval of personal photos using a deep autoencoder fusing visual features with speech annotations represented as word/paragraph vectors[END_REF], audio and video [START_REF] Ngiam | Multimodal deep learning[END_REF], but the systems exploiting them are not necessary bounded to those pairs. In all generality, methods for fusing modalities are often required when working with multimodal data. The most common approach consists in creating a joint multimodal representation by embedding every single-modal representations into a common representation space. There are two main groups of such approaches:

1. Multimodal approaches create a joint representation of the initially disjoint modalities or otherwise merge the initial modalities without necessarily providing a bidirectional mapping of the initial representation spaces to the new representation space and back. These approaches are typically used in retrieval and classification tasks where translating back from the multimodal representation to the single-modal ones is not required.

2. Crossmodal approaches focus on bidirectional mapping of the initial representations [START_REF] Feng | Cross-modal retrieval with correspondence autoencoder[END_REF], often by also creating a joint representation space in the process of doing so. They are able to map from one modality to another and back, as well as representing them in a joint representation space. These approaches can be used where crossmodal translation is required (e.g., multimodal query expansion, crossmodal retrieval) in addition to classification tasks.

In this work, we present a novel deep neural network architecture for crossmodal mapping and multimodal embedding. The seminal idea of the approach is to keep separate deep neural networks for each modality while tying the weights of the middle layers between the neural networks so as to yield a common multimodal representation. In this setting, the common middle layer acts as a common multimodal representation space that is attainable from either one of the modalities and from which we can attain either one of the modalities. As a proof of concept, experimental evidence of the benefits of the architecture is given by multimodal query expansion and multimodal retrieval tasks within the context of video hyperlinking. Extensions to other multimodal and crossmodal tasks remains straightforward.

METHODOLOGY

In this section, we analyze two methods for creating joint multimodal representations: multimodal autoencoders and

Multimodal Autoencoders

Two typical autoencoders are shown in Figure 1. The first (left) one illustrates a common approach that consists in concatenating the representations [START_REF] Ngiam | Multimodal deep learning[END_REF][START_REF] Lu | Semantic retrieval of personal photos using a deep autoencoder fusing visual features with speech annotations represented as word/paragraph vectors[END_REF] of the two modalities and training the autoencoder to reconstruct the data presented as input. The hidden layer in the middle is then used to obtain a joint multimodal representation (multimodal embedding).

The second (right) architecture is quite similar but has separate inputs and outputs (one for each modality) and separate hidden layers. One hidden layer in common is used for creating a joint multimodal representation. Sometimes, one modality is sporadically removed from the input to make the autoencoder learn to represent both modalities from one. The activations of the hidden layer are used as a multimodal joint representation. This enables autoencoders to also provide crossmodal mapping [START_REF] Ngiam | Multimodal deep learning[END_REF] in addition to a joint representation.

Bidirectional Representation Learning -Deep Neural Networks with Tied Weights

Previously described multimodal autoencoders include fully connected hidden layers that interact with both modalities. In other terms, the two modalities are available on the input and/or the two are reconstructed. The middle hidden layer represents a multimodal embedding and crossmodal translation is possible through this layer. However, mixing modalities is less optimal than directly mapping from one modality to another and back.

We propose a variation of the second autoencoder that implements a deep neural network translating from one modality to the other and back. Learning is performed in both directions: one modality is presented as an input and the other as the expected output while at the same time the second one is presented as input and the first one as expected output. This is equivalent to using two separate deep neural networks and tying them (sharing specific weight variables) to make them symmetrical, as illustrated in Figure 2. Implementation-wise the variables representing the weights are shared across the two networks and are in fact the same variables. Learning of the two crossmodal mappings is then performed simultaneously and they are forced to be as close as possible to each other's inverses by the symmetric architecture in the middle. A joint representation in the middle of the two crossmodal mappings is also formed while learning.

Formally, let h (j) i denote (the activation of) a hidden layer at depth j in network i (i = 1, 2; one for each modality), xm i the feature vector for modality i and o i the output of the network for modality i. In turn, for each network, W (j) i denotes the weight matrix of layer j and b (j) i the bias vector. Finally, we assume that each layer admits f as an activation function. The architecture is then defined by:

h (1) 1 = f (W (1) 1 × xm 1 + b (1) 1 ) h (1) 2 = f (W (1) 2 × xm 2 + b (1) 2 ) h (2) 1 = f (W (2) × h (1) 1 + b (2) 1 ) h (2) 2 = f (W (3) T × h (1) 2 + b (2) 2 ) h (3) 1 = f (W (3) × h (2) 1 + b (3) 1 ) h (3) 2 = f (W (2) T × h (2) 2 + b (3) 2 ) o1 = f (W (4) 1 × h (3) 1 + b (4) 1 ) o2 = f (W (4) 2 × h (3) 2 + b (4)
2 ) It is important to note that in the above equations, the weight matrices W (3) and W (2) are used twice due to weight tying, for computing h

(2) 1 , h (3) 2 and h (2) 2 , h (3) 1 respectively.
Training is performed by applying gradient descent to minimize the mean squared error of (o1, xm 2 ) and (o2, xm 1 ) thus effectively minimizing the reconstruction error in both directions and creating a joint representation in the middle, where both representations are projected.

Given such an architecture, crossmodal translation is done straightforwardly by presenting the first modality as xm i and obtaining the output in the representation space of the second modality as o i . A multimodal embedding is obtained by presenting one or both modalities (xm 1 and/or xm 2 ) at their respective inputs and reading the central hidden layers and/or h

(2)

1 . If both modalities are available, a concatenation of their embedding is used. If only one modality is available its embedding is duplicated. This allows for easy comparison, independently of modality availability.

EXPERIMENTS

Dataset and Initial Representations

We tested the symmetric bidirectional architecture within the framework of the video hyperlinking task using the Me-diaEval 2014 dataset and the respective groundtruth that was collected as part of the challenge [START_REF] Eskevich | The search and hyperlinking task at MediaEval[END_REF]. In this task, there are two main concepts: anchors and targets. Anchors represent segments of interest within videos that a user would like to know more about. Targets represent potential segments of interests that might or might not be related with a specific anchor. The goal is to hyperlink relevant targets for each anchor by using multimodal approaches. For each video, multiple data and modalities are available. In this work, we used two modalities: the automatic transcripts of the audio track and KU Leuven visual concepts. In practice, targets are not given and have to be defined automatically before assessing their relevance to each of the 30 anchors provided. Evaluation of the relevance is thus done post hoc on Amazon Mechanical Turk (AMT). In this paper, we consider a set of targets made of the top-10 targets that each participating team proposed for each anchor, along with the relevance judgments from AMT. In total, the dataset consists of 30 anchors, 10,809 targets and a ground truth with 12,340 anchor-target pairs (either related or unrelated). Interestingly, among the anchor and target segments, not all have both transcripts and visual concepts available. Table 1 illustrates the different cases of modality availability within the dataset. The task consists of using multimodal information to rank the targets by relevance for each anchor and comparing their relevance with the previously established groundtruth.

We chose to represent the transcripts and visual concepts of each anchor and target with a Word2Vec skip-gram model with hierarchical sampling [START_REF] Mikolov | Distributed representations of words and phrases and their compositionality[END_REF], a representation size of 100 and a window size of 5. The visual concepts were sorted previous to learning and the representations of the words and concepts found within a segment were averaged [START_REF] Campr | Comparing semantic models for evaluating automatic document summarization[END_REF]. This option worked best for our task (see Table 2) and additionally, our work is focused on cross-modality and joint multimodal representations, thus the choice of the input representations is not crucial. A standard cosine distance is used in all the experiments.

We implemented the two autoencoders described in Section 2.1 in Keras 1 and our proposed bidirectional symmetrical deep neural network in Lasagne 2 . All embeddings have a dimension of 100 as larger dimensions did not bring any significant improvement. The autoencoders we tested had architectures with 200-100-200 hidden layers, as had the symmetrical bidirectional networks. Other sizes were also tested but performed worse or not better and were not included in the comparison. Since this is unsupervised learning, the learning was performed on the part of the dataset that contains both transcripts and visual concepts and tested on the whole dataset. Table 3 reports the performance of the different methods. Using only transcripts yields a precision at 10 (P@10) of 58.67 % while using only visual concepts yields 50.00 %. Combining the two modalities in different linear combinations yields an improved P@10 to 61.32 %.

Crossmodal Query Expansion

Bidirectional symmetrical architectures enable crossmodal expansion where a missing modality is filled in by translating from the other one. If the transcript is not available for a segment, it is generated from the visual concepts and conversely. Using query expansion so that all segments have all modalities, we obtain 55.75 % for the visual concepts and 58.16 % for the transcripts. The first difference is significant (α = 0.001) improvement, while the second is not. This is due to the relatively small number of samples with one missing modality, so filling the missing modalities does not have a big influence.

Multimodal Embedding

A multimodal representation created with the two autoencoders (AE) of Sec. 2.1 yields 57.40 % (left part of Fig. 1) and 59.60 % (right part of Fig. 1). With symmetrical bidirectional networks, the two modalities are handled in the following manner to generate a multimodal embedding: 1) the modalities are both present for both the target and the anchor: all four representations are projected into the new multimodal representation space 2) one segment (anchor or target) has just one modality while the other has both: everything is projected into the new multimodal representation space, the one available modality of one segment being compared against the two modalities of the other 3) segments have only one modality: both modalities are projected into the new multimodal representation space.

When embedding only the transcripts to the new joint multimodal representation space with our proposed architecture, we obtain a P@10 of 70.43 %. Embedding only the visual concepts yields 54.92 %. The biggest improvement is achieved when embedding the two modalities, which yields 73.74 %, a significant (α = 10 -4 ) improvement over autoencoders of 14.14 in P@10.

Figure 3 illustrates the first two principal components of the two modalities, in their original representation spaces (left) and the first two principal components of their representation in the joint embedding space. The 30 anchors that are shown display a more ordered structure in the embedding space then in their original space, indicating that projections to a common representation space are successful.

Figure 4 shows a comparison of the two possible usage of bidirectional symmetrical architectures: for multimodal expansion (as a crossmodal translator) and as a multimodal embedder that translates the two modalities to a joint representation space that has been trained in a bidirectional way and that represents a common space for the crossmodal translation. The graph shows P@10 values for 10k training iterations. For reference, also the three baseline values are included (original transcripts, visual concepts and their linear combination). Multimodal query expansion improves P@10, not by a large margin due to the small number of samples missing one modality but it proves successful crossmodal translation. Embedding into the new joint representation space improves P@10 with a bigger margin.

CONCLUSIONS

Our proposed deep neural network architecture with tied weights is trained bidirectionally from the first modality to the second and from the second modality to the first and creates a crossmodal mapping between the two representation spaces that can be successfully used in multimodal query expansion. Due to its enforced symmetry, a joint multimodal embedding is also created that further improves multimodal data representation. We have shown that projecting into this space is significantly better than translating between the initial spaces, even when using only one modality. Projecting both modalities into the joint multimodal representation space yields the best results on our task. We expect that this architecture may provide similar improvements for other tasks that would be the subject of future works.

ACKNOWLEDGMENTS

Partially funded via the CominLabs excellence laboratory financed by the National Research Agency under reference ANR-10-LABX-07-01.

Figure 2 :

 2 Figure 2: Proposed architecture: training is done in both directions; a shared representation is created by tying the weights (sharing the variables) and enforcing symmetry in the central part

Figure 3 :

 3 Figure 3: Mapping the two modalities of the 30 anchors (first two components): a) in the original representation spaces b) after embedding them

Figure 4 :

 4 Figure 4: Comparison of P@10 for different training epochs (averages across multiple runs) and the baseline (original representations)

Table 1 :

 1 Available modalities for anchor/target pairs

					Targets
				Both	Transcripts Visual c. None
	Anchors	Both Transcripts Visual c. None	10,525 0 957 0	236 0 15 0	528 0 68 0	11 0 0 0
	h	(2) 1

Table 2 :

 2 Comparison of initial (single-modal) representations. Precision at 10 (%) obtained with TF-IDF, Word2Vec and Paragraph Vectors

	Representation	Transcripts Visual concepts
	LSI & TF-IDF	23.67	18.33
	Avg. Word2Vec	58.67	50.00
	PV-DM	45.00	45.33
	PV-DBOW	41.67	48.33
	Table 3: Comparison of the tested methods: preci-
	sion at 10 (%) and standard deviation
	Method	P@10	σ
		Baseline
	Only transcripts		58.67	-
	Only visual concepts		50.00	-
	Linear combination		61.32	3.1
	Autoencoders
	Embedded tran. and v. c. with AE #1	57.40	1.24
	Embedded tran. and v. c. with AE #2	59.60	0.65
	Bidirectional DNNs with tied weights
	Embedded transcripts		70.43	0.46
	Embedded visual concepts	54.92	0.99
	Embedded transcripts and visual c. 73.74 0.82
	Expanded transcripts		58.16	0.24
	Expanded visual concepts	55.75	0.13
	Expanded transcripts and visual c.	62.35	0.25

1 http://keras.io 2 https://github.com/Lasagne/Lasagne