Elsevier

Journal of Symbolic Computation

Volume 109, March–April 2022, Pages 1-30
Journal of Symbolic Computation

Polynomial-division-based algorithms for computing linear recurrence relations

https://doi.org/10.1016/j.jsc.2021.07.002Get rights and content

Abstract

Sparse polynomial interpolation, sparse linear system solving or modular rational reconstruction are fundamental problems in Computer Algebra. They come down to computing linear recurrence relations of a sequence with the Berlekamp–Massey algorithm. Likewise, sparse multivariate polynomial interpolation and multidimensional cyclic code decoding require guessing linear recurrence relations of a multivariate sequence.

Several algorithms solve this problem. The so-called Berlekamp–Massey–Sakata algorithm (1988) uses polynomial additions and shifts by a monomial. The Scalar-FGLM algorithm (2015) relies on linear algebra operations on a multi-Hankel matrix, a multivariate generalization of a Hankel matrix. The Artinian Gorenstein border basis algorithm (2017) uses a Gram-Schmidt process.

We propose a new algorithm for computing the Gröbner basis of the ideal of relations of a sequence based solely on multivariate polynomial arithmetic. This algorithm allows us to both revisit the Berlekamp–Massey–Sakata algorithm through the use of polynomial divisions and to completely revise the Scalar-FGLM algorithm without linear algebra operations.

A key observation in the design of this algorithm is to work on the mirror of the truncated generating series allowing us to use polynomial arithmetic modulo a monomial ideal. It appears to have some similarities with Padé approximants of this mirror polynomial.

As an addition from the paper published at the ISSAC conference, we give an adaptive variant of this algorithm taking into account the shape of the final Gröbner basis gradually as it is discovered. The main advantage of this algorithm is that its complexity in terms of operations and sequence queries only depends on the output Gröbner basis.

All these algorithms have been implemented in Maple and we report on our comparisons.

Introduction

The Berlekamp–Massey algorithm (BM ), introduced by Berlekamp (1968) and Massey (1969) is a fundamental algorithm in Coding Theory (Bose and Ray-Chaudhuri, 1960; Hocquenghem, 1959), and Computer Algebra. It allows one to perform efficiently sparse polynomial interpolation, sparse linear system solving or modular rational reconstruction.

In 1988, Sakata extended the BM algorithm to dimension n. This algorithm, known as the Berlekamp–Massey–Sakata algorithm (BMS ), can be used to compute a Gröbner basis of the zero-dimensional ideal of the relations satisfied by a sequence (Sakata, 1988, Sakata, 1990, Sakata, 2009). Analogously to dimension 1, the BMS algorithm allows one to decode cyclic codes in dimension n>1, an extension of Reed–Solomon's codes. Furthermore, the latest versions of the Sparse-FGLM algorithm rely heavily on the efficiency of the BMS algorithm to compute the change of ordering of a Gröbner basis (Faugère and Mou, 2011, Faugère and Mou, 2017).

In dimension 1, it is well known that the BM algorithm can be seen in a matrix form requiring to solve a linear Hankel system of size D, the order of the recurrence, see Jonckheere and Ma (1989), or the Levinson–Durbin method (Levinson, 1947; Wiener, 1964). If we let M(D) be a cost function for multiplying two polynomials of degree D, for instance M(D)O(DlogDloglogD) (Cantor and Kaltofen, 1991; Cooley and Tukey, 1965), then solving a linear Hankel system of size D comes down to performing a truncated extended Euclidean algorithm called on two polynomials of degree D (Blackburn, 2006; Brent et al., 1980; Dornstetter, 1987). More precisely, it can be done in O(M(D)logD) operations.

In Berthomieu et al., 2015, Berthomieu et al., 2017, the authors present the Scalar-FGLM algorithm, extending the matrix version of the BM algorithm for multidimensional sequences. It consists in computing the relations of the sequence through the computation of a maximal submatrix of full rank of a multi-Hankel matrix, a multivariate generalization of a Hankel matrix. Then, it returns the minimal Gröbner basis G of the ideal of relations satisfied by the sequence. These notions are recalled in Section 2. If we denote by S the staircase defined by G and T the input set of monomials containing

, then the complexity of the Scalar-FGLM algorithm is O((#T)ω), where 2ω3 is the linear algebra exponent. However, we do not know how to exploit the multi-Hankel structure to improve this complexity.

The Artinian Gorenstein border basis algorithm (AGbb ) was presented in Mourrain (2017) for computing a border basis B of the ideal of relations. It extends the algorithm of Berthomieu et al. (2015) using polynomial arithmetic allowing it to reach the better complexity O((#S+#B)#S#T) with the above notation.

Another viewpoint is that computing linear recurrence relations can be seen as computing Padé approximants of a truncation of the generating series i1,,in0wi1,,inx1i1xnin. In Fitzpatrick and Flynn (1992), the authors extend the extended Euclidean algorithm for computing multivariate Padé approximants. Given a polynomial P and an ideal B, find polynomials F and C such that P=FCmodB, where the leading monomials of F and C satisfy some constraints.

It is also worth noticing that we now know that both the BMS and the Scalar-FGLM algorithms are not equivalent, see Berthomieu and Faugère (2020), i.e. it is not possible to tweak one algorithm to mimic the behavior of the other. However, if the input sequence is linear recurrent and sufficiently many sequence terms are visited, then both algorithms compute a Gröbner basis of the zero-dimensional ideal of relations.

In the whole paper, we assume that the input sets of the Scalar-FGLM algorithm are the sets of all the monomials less than a given monomial and that these sets are finite. In order to improve the complexity of the algorithm, we will use polynomial arithmetic in all the operations. Even though they are not equivalent, this reduces the gap between the BMS and the Scalar-FGLM algorithms and provides a unified presentation.

In Section 3, we present the BM, the BMS and the Scalar-FGLM algorithms in a unified polynomial viewpoint. Using the mirror of the truncated generating series is a key ingredient letting us perform the computations modulo a specific monomial ideal B: a vector in the kernel of a multi-Hankel matrix is a polynomial C such that where P is the mirror of the truncated generating series,

denotes the leading monomial and tC is a monomial associated to C.

One interpretation of this is the computation of multivariate Padé approximants FC of P modulo B with different constraints than in Fitzpatrick and Flynn (1992) since we require that

is in a given set of terms and
satisfies equation (1).

This polynomial point of view allows us to design the Polynomial Scalar-FGLM algorithm (Algorithm 1) in Section 4 based on multivariate polynomial divisions. It computes, in a sense, a generating set of polynomials whose product with P modulo B must satisfy equation (1). If they do not, by polynomial divisions, we make new ones until finding minimal polynomials satisfying this constraint. It is worth noticing that in dimension 1, we recover the truncated extended Euclidean algorithm applied to the mirror polynomial of the generating series of the input sequence, truncated in degree D, and xD+1. All the examples are available on Berthomieu and Faugère (2021).

Our main result is Theorem 21, a simplified version of which is

Theorem 1

Let w be a sequence,be a total degree monomial ordering and a be a monomial. Let us assume that the reduced Gröbner basis G of the ideal of relations of w forand its staircase S satisfy

and for all ga, s=maxσa{σ,σga}, we have max(S)s. Then, the Polynomial Scalar-FGLM algorithm applied to w,and a terminates and computes a minimal Gröbner basis of the ideal of relations of w forin O(#S(#S+#G)#{σ,σa}) operations in the base field.

Let us also remark that the complexity bound is based on naive multivariate polynomial arithmetic and that this algorithm can benefit from improvements made in this domain.

In applications such as Gröbner bases change of orderings through the Sparse-FGLM algorithm, sequence queries are costly (Faugère and Mou, 2017). In Berthomieu et al. (2015), an adaptive variant of the Scalar-FGLM algorithm was designed aiming to minimize the number of sequence queries to recover the relations.

In Section 5, we show how we can transform the Adaptive Scalar-FGLM algorithm of Berthomieu et al. (2015) into an algorithm using polynomial arithmetic. This algorithm is output sensitive and probabilistic, like the Adaptive Scalar-FGLM algorithm is. That is, its main advantage is that its complexity only depends on the sizes of the computed staircase and Gröbner basis.

Theorem 2 see Theorem 26

Let w be a sequence,be a total degree monomial ordering.

Let us assume that calling the Adaptive Polynomial Scalar-FGLM algorithm on w andyields the Gröbner basis G and its staircase S.

Then, the Adaptive Polynomial Scalar-FGLM algorithm performs at most O((#S+#G)2#2S) operations in the base field and

table queries to recover G, where for a set T, 2T={tt,t,tT}.

Finally, in Section 6, we compare the Polynomial Scalar-FGLM algorithm with our implementations of the BMS, the Scalar-FGLM and the AGbb algorithms. Our algorithm performs always fewer arithmetic operations than the others starting from a certain size. Even for an example family favorable towards the BMS algorithm, our algorithm performs better.

Although we have compared the numbers of arithmetic operations, it would be beneficial to have an efficient implementation. This would be the first step into designing a similar efficient algorithm for computing linear recurrence relations with polynomial coefficients, extending the Beckermann–Labahn algorithm (Beckermann and Labahn, 1994) for computing multivariate Hermite–Padé approximants.

Amongst the changes from the ISSAC version of the paper (Berthomieu and Faugère, 2018), the main additions are a complete redesign of the Scalar-FGLM algorithm through polynomial arithmetic in Section 3.2.2 and a full description of the Adaptive Polynomial Scalar-FGLM algorithm, an adaptive variant of the Polynomial Scalar-FGLM algorithm using polynomial divisions as well, in Section 5.2. Generically, one could expect to make one relation with leading monomial mxi2 through a division of polynomials related to relations with leading monomials m and mxi. Yet, the naive approach given in Berthomieu and Faugère (2018, Section 5) could not do so as it does not perform any division. The Adaptive Polynomial Scalar-FGLM algorithm visits the monomials in the same order as the Adaptive Scalar-FGLM algorithm to recover the relations and replaces any linear algebra computations by polynomial ones, see Berthomieu et al. (2015). We also give the complexity of this algorithm in terms both of the number of operations and the number of sequence queries.

Furthermore, one of the main obstructions to the design of this adaptive variant is that at each step, some polynomials are updated. This update process adds terms supposed to be small with respect to the ordering. Yet, their leading terms were not stable during this update process.

Lastly, in Section 3, we now more clearly define what tC should be in Equation (1). This is a key point in the Polynomial Scalar-FGLM algorithm and a more complete description is available in Proposition 13.

Section snippets

Notation

We give a brief description of classical notation used in the paper.

From matrices to polynomials

Before detailing the unified polynomial viewpoint, we recall the linear algebra viewpoint of the BM, the BMS and the Scalar-FGLM algorithms.

A division-based algorithm

The goal is now to design an algorithm based on polynomial division to determine all the Cg for |-minimal g such that

is small enough, where Fg=PT+UCgmodB and F˜g is obtained from Fg as in Proposition 13.

We start with two sets of terms

and
so that
. We initialize B=(B1,,Bn)=(x1D1+1,,xnDn+1), RB1=[B1,0],,RBn=[Bn,0] and R1=[PT+U,1].

For any monomial g not in the ideal spanned by B and Rg=[Fg,Cg]=[PT+UCgmodB,Cg], by Proposition 13, Cg is a valid relation if g

An adaptive variant

In some applications, the actual size of the staircase, or at least an upper bound thereof, is known. While it provides an early termination criterion for the BMS, Scalar-FGLM and Polynomial Scalar-FGLM algorithms, this might fail to drastically reduce the number of table queries. Indeed, for the

ordering, whether the set of leading monomials of the Gröbner basis is {xn,,x2,x1d} or all the monomials of degree d: {xnd,xn1xnd1,,x1xnd1,,x1d}, the BMS algorithm requires to visit all

Experiments

In this section, we report on the number of counted arithmetic operations performed by the different algorithms for computing the Gröbner basis of the ideal of relations of some table families. They are counted using naive multiplications. To do so, we have a counter that is incremented by 1 each time a product in the base field is performed. For advanced functions that were not recoded, we use a formula. For instance, after a call to NormalForm(Fm,[F1,,Fr]), assuming the associated quotients

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

We thank the anonymous referees for their careful reading and their helpful comments to improve this paper. The authors are supported by the joint ANR-FWF ANR-19-CE48-0015 ECARP project, the ANR grants ANR-18-CE33-0011 Sesame and ANR-19-CE40-0018 De Rerum Natura projects, the PGMO grant P-2019-0018 CAMiSAdo and the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement N. 813211 (POEMA).

References (30)

  • E. Berlekamp

    Nonbinary BCH decoding

    IEEE Trans. Inf. Theory

    (1968)
  • J. Berthomieu et al.

    Linear algebra for computing Gröbner bases of linear recursive multidimensional sequences

  • J. Berthomieu et al.

    A polynomial-division-based algorithm for computing linear recurrence relations

  • J. Berthomieu et al.

    Experiments

  • S.R. Blackburn

    Fast rational interpolation, Reed-Solomon decoding, and the linear complexity profiles of sequences

    IEEE Trans. Inf. Theory

    (2006)
  • Cited by (1)

    • Guessing Gröbner bases of structured ideals of relations of sequences

      2022, Journal of Symbolic Computation
      Citation Excerpt :

      Given sufficiently many terms, the first two return a Gröbner basis of the ideal of relations while the third one returns a border basis of this ideal. Furthermore, in Berthomieu and Faugère (2018, 2022) the authors designed an algorithm extending both the Berlekamp–Massey–Sakata and the Scalar-FGLM algorithms using polynomial arithmetic and in Berthomieu and Faugère (2016), they extended the Scalar-FGLM algorithm for guessing relations with polynomial coefficients. However, none of these algorithms were designed to take the structure of the table terms into account.

    1

    Laboratoire d'Informatique de Paris 6, Sorbonne Université, boîte courrier 169, 4 place Jussieu, F-75252 Paris Cedex 05, France.

    View full text