Polynomial-division-based algorithms for computing linear recurrence relations
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 , 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 be a cost function for multiplying two polynomials of degree D, for instance (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 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 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 and T the input set of monomials containing , then the complexity of the Scalar-FGLM algorithm is , where 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 of the ideal of relations. It extends the algorithm of Berthomieu et al. (2015) using polynomial arithmetic allowing it to reach the better complexity 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 . 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 , 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 is a monomial associated to C.
One interpretation of this is the computation of multivariate Padé approximants 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 . 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 of the ideal of relations of w for ≺ and its staircase S satisfy and for all , , we have . 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 for ≺ in operations in the base field.
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 and ≺ yields the Gröbner basis and its staircase S. Then, the Adaptive Polynomial Scalar-FGLM algorithm performs at most operations in the base field and table queries to recover , where for a set T, .
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 through a division of polynomials related to relations with leading monomials m and . 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 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 for |-minimal g such that is small enough, where and is obtained from as in Proposition 13.
We start with two sets of terms and so that . We initialize , and .
For any monomial g not in the ideal spanned by B and , by Proposition 13, is a valid relation if
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 or all the monomials of degree d: , 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 , 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)
- et al.
Linear algebra for computing Gröbner bases of linear recursive multidimensional sequences
J. Symb. Comput.
(2017) - et al.
In-depth comparison of the Berlekamp–Massey–Sakata and the scalar-FGLM algorithms: the adaptive variants
J. Symb. Comput.
(2020) - et al.
On a class of error correcting binary group codes
Inf. Control
(1960) - et al.
Fast solution of Toeplitz systems of equations and computation of Padé approximants
J. Algorithms
(1980) - et al.
Sparse FGLM algorithms
J. Symb. Comput.
(2017) - et al.
A Gröbner basis technique for Padé approximation
J. Symb. Comput.
(1992) - et al.
A simple Hankel interpretation of the Berlekamp-Massey algorithm
Linear Algebra Appl.
(1989) Finding a minimal set of linear recurring relations capable of generating a given finite two-dimensional array
J. Symb. Comput.
(1988)Extension of the Berlekamp-Massey algorithm to N dimensions
Inf. Comput.
(1990)- et al.
A uniform approach for the fast computation of matrix-type Pade approximants
SIAM J. Matrix Anal. Appl.
(1994)
Nonbinary BCH decoding
IEEE Trans. Inf. Theory
Linear algebra for computing Gröbner bases of linear recursive multidimensional sequences
A polynomial-division-based algorithm for computing linear recurrence relations
Experiments
Fast rational interpolation, Reed-Solomon decoding, and the linear complexity profiles of sequences
IEEE Trans. Inf. Theory
Cited by (1)
Guessing Gröbner bases of structured ideals of relations of sequences
2022, Journal of Symbolic ComputationCitation 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.