Fast algorithm for border bases of Artinian Gorenstein algebras

TL;DR

This paper introduces an efficient algorithm to compute border bases of Artinian Gorenstein algebras from multi-index sequences, extending the Berlekamp-Massey-Sakata algorithm with improved complexity for applications like coding theory and tensor decomposition.

Contribution

The paper presents a novel algorithm that computes border bases of Artinian Gorenstein algebras from sequences, enhancing existing methods with better efficiency and broader applications.

Findings

Algorithm computes generators of recurrence relations efficiently.

Provides explicit border basis and multiplication tables for the algebra.

Demonstrates practical effectiveness through benchmarks.

Abstract

Given a multi-index sequence $$\sigma$$, we present a new efficient algorithm to compute generators of the linear recurrence relations between the terms of $$\sigma$$. We transform this problem into an algebraic one, by identifying multi-index sequences, multivariate formal power series and linear functionals on the ring of multivariate polynomials. In this setting, the recurrence relations are the elements of the kerne l$I$\sigma$$ of the Hankel operator $H$\sigma$$ associated to $$\sigma$$. We describe the correspondence between multi-index sequences with a Hankel operator of finite rank and Artinian Gorenstein Algebras. We show how the algebraic structure of the Artinian Gorenstein algebra $A$\sigma$$ associated to the sequence $$\sigma$$ yields the structure of the terms $\sigma$$\alpha$ for all $$\alpha$ $\in$ N n$. This structure is explicitly given by a border basis of $A$\sigma$$, which is presented as a quotient of the polynomial ring $K[x 1 ,. .. , xn$] by the kernel $I$\sigma$$ of the Hankel operator $H$\sigma$$. The algorithm provides generators of $I$\sigma$$ constituting a border basis, pairwise orthogonal bases of $A$\sigma$$ and the tables of multiplication by the variables in these bases. It is an extension of Berlekamp-Massey-Sakata (BMS) algorithm, with improved complexity bounds. We present applications of the method to different problems such as the decomposition of functions into weighted sums of exponential functions, sparse interpolation, fast decoding of algebraic codes, computing the vanishing ideal of points, and tensor decomposition. Some benchmarks illustrate the practical behavior of the algorithm.

Equations

$$p=\sum _{\alpha \in A\subset \text{mathbbm}{N}^n}{p}_{\alpha }{\mathbf{x}}^{\alpha }\in \text{mathbbm}K[\mathbf{x}]\mapsto ⟨\sigma \mid p⟩=\sum _{\alpha \in A\subset \text{mathbbm}{N}^n}{p}_{\alpha }{\sigma }_{\alpha }.$$

$$⟨p\star \sigma \mid q⟩$$

$$\text{mathbbm}K[\mathbf{x}]\times \text{mathbbm}K[\mathbf{x}]$$

$$(p,q)$$

$$\sigma =({\sigma }_{\alpha }{)}_{\alpha \in \text{mathbbm}{N}^n}\in \text{mathbbm}{K}^{\text{mathbbm}{N}^n}\mapsto \sigma (\mathbf{z})=\sum _{\alpha \in \text{mathbbm}{N}^n}{\sigma }_{\alpha }{\mathbf{z}}^{\alpha }\in \text{mathbbm}K[[\mathbf{z}]].$$

$$\sigma =\sum _{i=1}^r\sum _{\alpha \in A_i}\frac{{\omega }_{i,\alpha }{\mathbf{z}}^{\alpha }}{{\prod }_{j=1}^n(1-{\xi }_{i,j}{z}_j{)}^{1+{\alpha }_j}}$$

$$H_{\sigma }:\text{mathbbm}K[\mathbf{x}]$$

$$p=\sum _{\beta \in B}{p}_{\beta }{\mathbf{x}}^{\beta }$$

$$[H_{\sigma }]$$

$$H_{\sigma }^{V,V^{\prime }}:V$$

$$p(\mathbf{x})$$

$$[H_{\sigma }^{B,B^{\prime }}]=(⟨\sigma |b_j{b}_i^{\prime }⟩{)}_{1⩽i,j⩽r}.$$

$$[H_{\sigma }^{B,B^{\prime }}]$$

$$I^{\perp }=\{\sigma \in \text{mathbbm}K[\mathbf{x}{]}^{\ast }\mid \forall p\in I,⟨\sigma |p⟩=0\}.$$

$$0\to I\to \text{mathbbm}K[\mathbf{x}]\stackrel{H_{\sigma }}{\to }{\mathcal{A}}^{\ast }\to 0$$

$$\langle p_{\beta},q_{\beta^{\prime}}\rangle_{\sigma}=\left\{\begin{array}[]{ll}1&\operatorname{if}\beta=\beta^{\prime}\\ 0&\operatorname{otherwise}.\end{array}\right.$$

$$p_{\alpha }=\mathrm{proj}({\mathbf{x}}^{\alpha },[p_{{}_{\beta }}{]}_{\beta \in {\mathbf{b}}_{\prec \alpha }},[m_{\beta }{]}_{\beta \in {\mathbf{b}}_{\prec \alpha }})$$

$$q_{\alpha }=m_{\alpha }-\sum _{\beta \in {\mathbf{b}}_{\prec \alpha }}⟨p_{\beta },m_{\alpha }{⟩}_{\sigma }{q}_{\beta }.$$

$$⟨m_{\beta }{⟩}_{\beta \in {\mathbf{b}}_{\preccurlyeq \alpha }}$$

$$⟨p_{\alpha },q_{\alpha }{⟩}_{\sigma }$$

$$p_{\alpha }={\mathbf{x}}^{\alpha }+\sum _{\beta \prec \alpha }{\lambda }_{\beta }{p}_{\beta }={\mathbf{x}}^{\alpha }+b_{\alpha }$$

$$H_{\sigma}^{V,V^{\prime}}=\left(\begin{array}[]{cc}L_{r}&0\\ \ast&0\end{array}\right)$$

$$M_k:=(⟨x_k{p}_{{\beta }_j},q_{{\beta }_i}{⟩}_{\sigma }{)}_{1⩽i,j⩽r}=(⟨\sigma |x_k{p}_{{\beta }_j}{q}_{{\beta }_i}⟩{)}_{1⩽i,j⩽r}.$$

$$⟨p_{\beta },{\mathbf{x}}^{\gamma }{⟩}_{\sigma }=⟨p_{\beta },{\mathbf{x}}^{\delta }{\mathbf{x}}^{{\gamma }^{\prime }}{⟩}_{\sigma }=⟨p_{\beta },p_{\delta }{\mathbf{x}}^{{\gamma }^{\prime }}{⟩}_{\sigma }-⟨p_{\beta },b_{\delta }{\mathbf{x}}^{{\gamma }^{\prime }}{⟩}_{\sigma }.$$

$$[{\xi }_{i_j}^{\alpha }{]}_{|\alpha |⩽d,1⩽j⩽r}({\omega }_{i_j})=({\sigma }_{\alpha }{)}_{|\alpha |⩽d}.$$

$${\sigma }_{\alpha }=h({\zeta }_1^{{\alpha }_1},\dots ,{\zeta }_n^{{\alpha }_n})=\sum _{i=1}^r{\omega }_i({\zeta }_1^{{\gamma }_{i,1}}{)}^{{\alpha }_1}\cdots ({\zeta }_n^{{\gamma }_{i,n}}{)}^{{\alpha }_n}$$

$$t=\sum _{{\alpha }_0+{\alpha }_1+\cdots +{\alpha }_n=d}{t}_{\alpha }\left(\genfrac{}{}{0ex}{}{d}{\alpha }\right)x_0^{{\alpha }_0}{x}_1^{{\alpha }_1}\dots x_n^{{\alpha }_n}$$

$$t=\sum _{i=1}^r{\omega }_i({\xi }_{i,0}{x}_0+{\xi }_{i,1}{x}_1+\cdots +{\xi }_{i,n}{x}_n{)}^d$$

$${\sigma }_{\alpha }:=t_{(d-{\alpha }_1\cdots -{\alpha }_n,{\alpha }_1,\dots ,{\alpha }_n)}=\sum _{i=1}^r{\omega }_i{\xi }_{i,1}^{{\alpha }_1}\cdots {\xi }_{i,n}^{{\alpha }_n}=\sum _{i=1}^r{\omega }_i{\xi }_i^{\alpha }$$

$$x_1{p}_1$$

$$x_1{p}_{x_1}$$

$$x_1{p}_{x_2}$$

$$M_{1}=\left[\begin{array}[]{ccc}\frac{5}{4}&-\frac{5}{16}&0\\ 1&\frac{91}{20}&\frac{96}{25}\\ 0&-1&\frac{1}{5}\end{array}\right].$$

Full text

1 Introduction

Discovering hidden structures from probing or sampling is a problem which appears in many contexts and in many applications. An interesting instance of this general problem is recovering the structure of a sequence of values, from the knowledge of some of its terms. It consists in guessing any term of the sequence from the first known terms. A classical way to tackle this problem, which goes back to Bernoulli, is to find linear recurrence relations between the first terms of a sequence, to compute the roots of the associated characteristic polynomial and to deduce the expression of any term of the sequence from these roots.

In this paper, we consider the structure discovering problem for multi-index sequences $\sigma=(\sigma_{\alpha})_{\alpha\in\mathbbm{N}^{n}}\in\mathbbm{K}^{\mathbbm{N}^{n}}$ of values in a field $\mathbbm{K}$. Given a finite set of values $\sigma_{\alpha}$ for $\alpha\in\boldsymbol{a}\subset\mathbbm{N}^{n}$, we want to guess a formula for the general terms of the sequence $\sigma$. An important step of this approach is to compute characteristic polynomials of the sequence $\sigma=(\sigma_{\alpha})_{\alpha\in\mathbbm{N}^{n}}$. They correspond to multi-index recurrence relations with constant coefficients between the terms of $\sigma$. The ideal of these recurrence relation polynomials define an Artinian Gorenstein algebra. We present a fast algorithm to compute a border basis of this ideal from the first terms of the sequence $\sigma$. This method also yields a basis of the Artinian Gorenstein algebra as well as its multiplicative structure.

Related works

The approach that we present is related to Prony’s method in the univariate case and to its variants [33], [29], [17], [8], … and to the more recent extensions in the multivariate case [1], [27], [18], [32]. Linear algebra tools are used to determine a basis of the quotient algebra $\mathcal{A}_{\sigma}$ or to compute an $H$-basis for the presentation of $\mathcal{A}_{\sigma}$. An analysis of the complexity of these approaches yields bounds in $\mathcal{O}(\tilde{s}^{3})$ (or $\mathcal{O}(\tilde{s}^{\omega})$), where $\tilde{s}$ is the size of the Hankel matrices involved in these methods, typically the number $\binom{d^{\prime}+n}{n}$ of monomials in degree at most $d^{\prime}<\frac{d}{2}$ if the terms of the sequence are known up to the degree $d$. The problem is also related to Padé approximants, well investigated in the univariate case [11], [3], [35], but much less developed in the multivariate case [28], [13].

Finding recurrence relations is a problem which is also well developed in the univariate case. Berlekamp [5] and Massey [20] proposed an efficient algorithm to compute such recurrence relations, with a complexity in $\mathcal{O}(r^{2})$ where $r$ is the size of the minimal recurrence relation. Exploiting further the properties of Hankel matrices, the complexity of computing recurrence relations can be reduced in the univariate case to $\tilde{\mathcal{O}}(r)$.

Sakata extended Berlekamp-Massey algorithm to multi-index sequences, computing a Gröbner basis of the polynomials in the kernel of a multi-index Hankel matrix [31]. See also [30] for an analysis and overview of the algorithm. The computation of multivariate linear recurrence relations have been further investigated, e.g. in [16] and more recently in [7], where the coefficients of the Gröbner basis are computed by solving multi-index Hankel systems.

Contributions

We translate the structure discovering problem into an algebraic setting, by identifying multi-index sequences of values, generating formal power series and linear functionals on the ring of polynomials. Through this identification, we associate to a multi-index sequence ${\sigma}$, a Hankel operator $H_{\sigma}$ which kernel $I_{\sigma}$ defines an Artinian Gorenstein Algebra $\mathcal{A}_{\sigma}$ when $H_{\sigma}$ is of finite rank. We present a new efficient algorithm to compute the algebraic structure of $\mathcal{A}_{\sigma}$, using the first terms $\sigma_{\alpha}$ for $\alpha\in\boldsymbol{a}\subset\mathbbm{N}^{n}$. The structure $\mathcal{A}_{\sigma}$ is described by a border basis of the ideal $I_{\sigma}$.

This algorithm is an extension of the Berlekamp-Massey-Sakata (BMS) algorithm. It computes border bases of the recurrence relations, which are more general than Gröbner bases. They also offer a better numerical stability [25] in the solving steps required to address the decomposition problem. The algorithm, based on a Gram-Schmidt orthogonalisation process, is simplified. The complexity bound also improves the previously known bounds for computing such recurrence relations. We show that the arithmetic complexity of computing a border basis is in $\mathcal{O}((r+\delta)rs)$ where $r$ is the number of roots of $I_{\sigma}$ (counted with multiplicities), $\delta$ is the size of the border of the monomial basis and $s$ is the number of known terms of the sequence $\sigma$.

The algorithm outputs generators of the recurrence relations, a monomial basis, an orthogonal basis and the tables of multiplication by the variables in this basis of $\mathcal{A}_{\sigma}$. The structure of the terms of the sequence $\sigma$ can be deduced from this output, by applying classical techniques for solving polynomial systems from tables of multiplication. We show how the algorithm can be applied to different problems such as the decomposition of functions into weighted sums of exponential functions, sparse interpolation, fast decoding of algebraic codes, vanishing ideal of points, and tensor decomposition.

Notation

Let $\mathbbm{K}$ be a field, $\bar{\mathbbm{K}}$ its algebraic closure, $\mathbbm{K}[x_{1},\ldots,x_{n}]=\mathbbm{K}[\boldsymbol{x}]$ be the ring of polynomials in the variables $x_{1},\ldots,x_{n}$ with coefficients in the field $\mathbbm{K}$, $\mathbbm{K}[[y_{1},$ $\ldots,y_{n}]]=\mathbbm{K}[\boldsymbol{y}]$ be the ring of formal power series in the variables $y_{1},\ldots,y_{n}$ with coefficients in $\mathbbm{K}$. We denote by $\mathbbm{K}^{\mathbbm{N}^{n}}$ the set of sequences $\sigma=(\sigma_{\alpha})_{\alpha\in\mathbbm{N}^{n}}$  of numbers $\sigma_{\alpha}\in\mathbbm{K}$, indexed by $\mathbbm{N}^{n}$. $\forall\alpha=(\alpha_{1},\ldots,\alpha_{n})\in\mathbbm{N}^{n},$ $\alpha!=\prod_{i=1}^{n}\alpha_{i}!$, $\boldsymbol{x}^{\alpha}=\prod_{i=1}^{n}x_{i}^{\alpha_{i}}$. The monomials in $\mathbbm{K}[\boldsymbol{x}]$ are the elements of the form $\boldsymbol{x}^{\alpha}$ for $\alpha\in\mathbbm{N}^{n}$. For a set $B\subset\mathbbm{K}[\boldsymbol{x}]$, $B^{+}=\cup_{i=1}^{n}x_{i}B\cup B$, $\partial B=B^{+}\setminus B$. A set $B$ of monomials of $\mathbbm{K}[\boldsymbol{x}]$ is connected to $1$, if $1\in B$ and for $\boldsymbol{x}^{\beta}\in B$ different from $1$, there exists $\boldsymbol{x}^{\beta^{\prime}}\in B$ and $i\in[1,n]$ such that $\boldsymbol{x}^{\beta}=x_{i}\boldsymbol{x}^{\beta^{\prime}}$. For $F\subset\mathbbm{K}[\boldsymbol{x}]$, $\langle F\rangle$ is the vector space of $\mathbbm{K}[\boldsymbol{x}]$ spanned by $F$ and $(F)$ is the ideal generated by $F$. For $V,V^{\prime}\subset\mathbbm{K}[\boldsymbol{x}]$, $V\cdot V^{\prime}$ is the set of products of an element of $V$ by an element of $V^{\prime}$.

2 Polynomial-Exponential series

In this section, we recall the correspondence between sequences $\sigma=(\sigma_{\alpha})_{\alpha\in\mathbbm{N}^{n}}\in\mathbbm{K}^{\mathbbm{N}^{n}}$ associated to polynomial-exponential series and Artinian Gorenstein Algebras.

2.1 Duality

A sequence $\sigma=(\sigma_{\alpha})_{\alpha\in\mathbbm{N}^{n}}\in\mathbbm{K}^{\mathbbm{N}^{n}}$ is naturally associated to a linear form operating on polynomials, that is, an element of $\operatorname{Hom}_{\mathbbm{K}}(\mathbbm{K}[\boldsymbol{x}],\mathbbm{K})=\mathbbm{K}[\boldsymbol{x}]^{\ast}$, as follows:

[TABLE]

This correspondence is bijective since a linear form $\sigma\in\mathbbm{K}[\boldsymbol{x}]^{\ast}$ is uniquely defined by the sequence $\sigma_{\alpha}=\langle\sigma\mid\boldsymbol{x}^{\alpha}\rangle$ for $\alpha\in\mathbbm{N}^{n}$. The coefficients $\sigma_{\alpha}=\langle\sigma\mid\boldsymbol{x}^{\alpha}\rangle$ for $\alpha\in\mathbbm{N}^{n}$ are also called the moments of $\sigma$. Hereafter, we will identify $\mathbbm{K}^{\mathbbm{N}^{n}}$ with $\mathbbm{K}[\boldsymbol{x}]^{\ast}=\operatorname{Hom}_{\mathbbm{K}}(\mathbbm{K}[\boldsymbol{x}],\mathbbm{K})$.

The dual space $\mathbbm{K}[\boldsymbol{x}]^{\ast}$ has a natural structure of $\mathbbm{K}[\boldsymbol{x}]$-module, defined as follows: $\forall\sigma\in\mathbbm{K}[\boldsymbol{x}]^{\ast},\forall p,q\in\mathbbm{K}[\boldsymbol{x}]$,

[TABLE]

We check that $\forall\sigma\in\mathbbm{K}[\boldsymbol{x}]^{\ast},\forall p,q\in\mathbbm{K}[\boldsymbol{x}]$, $(pq)\star\sigma=p\star(q\star\sigma)$. See e.g. [15], [21] for more details.

For any $\sigma\in\mathbbm{K}[\boldsymbol{x}]^{\ast}$, the inner product associated to $\sigma$ on $\mathbbm{K}[\boldsymbol{x}]$ is defined as follows:

[TABLE]

Sequences in $\mathbbm{K}^{\mathbbm{N}^{n}}$ are also in correspondence with series in $\mathbbm{K}[[\boldsymbol{z}]]$, via the so-called $\boldsymbol{z}$-transform:

[TABLE]

If $\mathbbm{K}$ is a field of characteristic [math], we can identify the sequence $\sigma=(\sigma_{\alpha})_{\alpha\in\mathbbm{N}^{n}}\in\mathbbm{K}^{\mathbbm{N}^{n}}$with the series $\sigma(\boldsymbol{y})=\sum_{\alpha\in\mathbbm{N}^{n}}\sigma_{\alpha}\frac{\boldsymbol{y}^{\alpha}}{\alpha!}\in\mathbbm{K}[[\boldsymbol{y}]].$ Using this identification, we have  $\forall p\in\mathbbm{K}[\boldsymbol{x}],\forall\sigma\in\mathbbm{K}[[\boldsymbol{y}]]$, $p\star\sigma(\boldsymbol{y})=p(\partial_{y_{1}},\ldots,\partial_{y_{n}})(\sigma(\boldsymbol{y}))$.

Through these identifications, the dual basis of the monomial basis $(\boldsymbol{x}^{\alpha})_{\alpha\in\mathbbm{N}^{n}}$ is $(\boldsymbol{z}^{\alpha})_{\alpha\in\mathbbm{N}^{n}}$ in $\mathbbm{K}[[\boldsymbol{z}]]$ and $\left(\frac{\boldsymbol{y}^{\alpha}}{\alpha!}\right)_{\alpha\in\mathbbm{N}^{n}}$ in $\mathbbm{K}[[\boldsymbol{y}]]$.

Among the elements of $\operatorname{Hom}(\mathbbm{K}[\boldsymbol{x}],\mathbbm{K})$, we have the evaluation $\boldsymbol{e}_{{}_{\xi}}:p(\boldsymbol{x})\in\mathbbm{K}[\boldsymbol{x}]\mapsto p(\xi)\in\mathbbm{K}$ at a point $\xi\in\mathbbm{K}^{n}$, which corresponds to the sequence $(\xi^{\alpha})_{\alpha\in\mathbbm{N}^{n}}$ or to the series $\boldsymbol{e}_{{}_{\xi}}(\boldsymbol{z})=\sum_{\alpha\in\mathbbm{N}^{n}}\xi^{\alpha}\boldsymbol{z}^{\alpha}=\hskip 2.5pt\frac{1\hskip 1.75pt}{\prod_{i=1}^{n}(1-\xi_{i}z_{i})}\in\mathbbm{K}[[\boldsymbol{z}]]$, or to the series $\boldsymbol{e}_{{}_{\xi}}(\boldsymbol{y})=\sum_{\alpha\in\mathbbm{N}^{n}}\xi^{\alpha}\frac{\boldsymbol{y}^{\alpha}}{\alpha!}=\hskip 2.5pte^{\xi_{1}y_{1}+\cdots+\xi_{n}y_{n}}=e^{\langle\xi,\mathbf{y}\rangle}$ in $\mathbbm{K}[[\boldsymbol{y}]]$. These series belong to the more general family $\mathcal{{POLYEXP}}$ of polynomial-exponential series $\sigma=\sum_{i=1}^{r}\omega_{i}\boldsymbol{e}_{\xi_{i}}\in\mathbbm{K}[[\boldsymbol{y}]]$ with $\xi_{i}\in\mathbbm{K}^{n},\omega_{i}\in\mathbbm{K}[\boldsymbol{y}]$. This set corresponds in $\mathbbm{K}[[\boldsymbol{z}]]$ to the set of series of the form

[TABLE]

with $\xi_{i}\in\mathbbm{K}^{n},\omega_{i,\alpha}\in\mathbbm{K},\alpha\in A_{i}\subset\mathbbm{N}^{n}$ and $A_{i}$ finite.

Definition 2.1

For a subset $D\subset\mathbbm{K}[[\boldsymbol{y}]]$, the inverse system generated by $D$ is the vector space spanned by the elements $p\star\delta$ for $\delta\in D$, $p\in\mathbbm{K}[\boldsymbol{x}]$, that is, by the elements in $D$ and all their derivatives.

*For $\omega\in\mathbbm{K}[\boldsymbol{y}]$, we denote by ${\mu}(\omega)$ the dimension of the inverse system of $\omega$, generated by $\omega$ and all its derivatives. For $\sigma=\sum_{i=1}^{r}\omega_{i}\boldsymbol{e}_{\xi_{i}}\in\mathcal{{POLYEXP}}(\boldsymbol{y})$, ${\mu}(\sigma)=\sum_{i=1}^{r}{\mu}(\omega_{i})$. *

2.2 Hankel operators

The external product $\star$ allows us to define a Hankel operator as a multiplication operator by a dual element $\in\mathbbm{K}[\boldsymbol{x}]^{\ast}$:

Definition 2.2

The Hankel operator associated to an element $\sigma\in\mathbbm{K}[\boldsymbol{x}]^{\ast}=\mathbbm{K}^{\mathbbm{N}^{n}}$ is

[TABLE]

*Its kernel is denoted $I_{\sigma}$. We say that the series $\sigma$ has finite rank $r\in\mathbbm{N}$ if $\operatorname{rank}H_{\sigma}=r<\infty$. *

As $\forall p,q\in\mathbbm{K}[\boldsymbol{x}]$,  $pq\star\sigma=p\star(q\star\sigma)$, we easily check that $I_{\sigma}=\ker H_{\sigma}$ is an ideal of $\mathbbm{K}[\boldsymbol{x}\mathbf{}]$ and that $\mathcal{A}_{\sigma}=\mathbbm{K}[\boldsymbol{x}]/I_{\sigma}$ is an algebra.

Given a sequence $\sigma=(\sigma_{\alpha})_{\alpha\in\mathbbm{N}^{n}}\in\mathbbm{K}^{\mathbbm{N}^{n}}$, the kernel of $H_{\sigma}$ is the set of polynomials $p=\sum_{\beta\in B}p_{\beta}\boldsymbol{x}^{\beta}$ such that $\sum_{\beta\in B}p_{\beta}\sigma_{\alpha+\beta}=0$ for all $\alpha\in\mathbbm{N}^{n}$. This kernel is the set of linear recurrence relations of the sequence $\sigma=(\sigma_{\alpha})_{\alpha\in\mathbbm{N}^{n}}$.

Remark 2.3

The matrix of the operator $H_{\sigma}$ in the bases $(\boldsymbol{x}^{\alpha})_{\alpha\in\mathbbm{N}^{n}}$ and its dual basis $(\boldsymbol{z}^{\alpha})_{\alpha\in\mathbbm{N}^{n}}$ is

[TABLE]

The coefficients of $[H_{\sigma}]$ depend only the sum of the multi-indices indexing the rows and columns, which explains why it is called a Hankel operator.

In the reconstruction problem, we are dealing with truncated series with known coefficients $\sigma_{\alpha}$ for $\alpha$ in a subset $\boldsymbol{a}$ of $\mathbbm{N}^{n}$. This leads to the definition of truncated Hankel operators.

Definition 2.4

For  two vector spaces $V,V^{\prime}\subset\mathbbm{K}[\boldsymbol{x}]$ and  $\sigma\in\langle V\cdot V^{\prime}\rangle^{\ast}\subset\mathbbm{K}[\boldsymbol{x}]^{\ast}$, the truncated Hankel operator on $(V,V^{\prime})$, denoted by $H_{\sigma}^{V,V^{\prime}}$, is the following map:

[TABLE]

If $B=\{b_{1},\ldots,b_{r}\}$ (resp. $B^{\prime}=\{b_{1}^{\prime},\ldots,b_{r}^{\prime}\}$) is a basis of $V$ (resp. $V^{\prime}$), then the matrix of the operator $H_{\sigma}^{V,V^{\prime}}$ in $B$ and the dual basis of $B^{\prime}$ is

[TABLE]

If $B$ and $B^{\prime}$ are monomial sets, we obtain the so-called truncated moment matrix of $\sigma$:

[TABLE]

(identifying a monomial $\boldsymbol{x}^{\beta}$ with its exponent $\beta$). These structured matrices share with the classical univariate Hankel matrices many interesting properties (see e.g. in [23]).

2.3 Artinian Gorenstein algebra

A $\mathbbm{K}$-algebra $\mathcal{A}$ is Artinian if $\dim_{\mathbbm{K}}(\mathcal{A})<\infty$. It can be represented as the quotient $\mathcal{A}=\mathbbm{K}[\boldsymbol{x}]/I$ of a polynomial ring $\mathbbm{K}[\boldsymbol{x}]$ by a (zero-dimension) ideal $I\subset\mathbbm{K}[\boldsymbol{x}]$.

A classical result states that the quotient algebra $\mathcal{A}=\mathbbm{K}[\boldsymbol{x}]/I$ is finite dimensional, i.e. Artinian, iff $\mathcal{V}_{\bar{\mathbbm{K}}}(I)$ is finite, that is, $I$ defines a finite number of (isolated) points in $\bar{\mathbbm{K}}^{n}$ (see e.g. [12][Theorem 6] or [14][Theorem 4.3]).

The dual $\mathcal{A}^{\ast}=\operatorname{Hom}_{\mathbbm{K}}(\mathcal{A},\mathbbm{K})$ of $\mathcal{A}=\mathbbm{K}[\boldsymbol{x}]/I$ is naturally identified with the sub-space

[TABLE]

A Gorenstein algebra is defined as follows:

Definition 2.5

*A $\mathbbm{K}$-algebra $\mathcal{A}$ is Gorenstein if $\exists\sigma\in\mathcal{A}^{\ast}=\operatorname{Hom}_{\mathbbm{K}}(\mathcal{A},\mathbbm{K})$ such that $\forall\rho\in\mathcal{A}^{\ast},\exists a\in\mathcal{A}$ with $\rho=a\star\sigma$ and $a\star\sigma=0$ implies $a=0$. *

In other words, $\mathcal{A}=\mathbbm{K}[\boldsymbol{x}]/I$ is Gorenstein iff $\mathcal{A}^{\ast}=\{p\star\sigma\mid p\in\mathbbm{K}[\boldsymbol{x}]\}=\operatorname{im}H_{\sigma}$ and $p\star\sigma=0$ implies $p\in I$. Equivalently, $\mathcal{A}=\mathbbm{K}[\boldsymbol{x}]/I$ is Gorenstein iff there exists $\sigma\in\mathbbm{K}[\boldsymbol{x}]^{\ast}$ such that we have the exact sequence:

[TABLE]

so that $H_{\sigma}$ induces an isomorphism between $\mathcal{A}=\mathbbm{K}[\boldsymbol{x}]/I$ and $\mathcal{A}^{\ast}$. In other words, a Gorenstein algebra $\mathcal{A}$ is the quotient of a polynomial ring by the kernel of a Hankel operator, or equivalently by an ideal of recurrence relations of a multi-index sequence.

An Artinian Gorenstein can thus be described by an element $\sigma\in\mathbbm{K}[\boldsymbol{x}]^{\ast}$, such that $\operatorname{rank}H_{\sigma}=\dim\mathcal{A}^{\ast}=\dim\mathcal{A}$ is finite. In the following, we will assume that the Artinian Gorenstein algebra is given by such an element $\sigma\in\mathbbm{K}[\boldsymbol{x}]^{\ast}\equiv\mathbbm{K}^{\mathbbm{N}^{n}}$. The corresponding algebra will be $\mathcal{A}_{\sigma}=\mathbbm{K}[\boldsymbol{x}]/I_{\sigma}$ where $I_{\sigma}=\ker H_{\sigma}$.

By a multivariate generalization of Kronecker’s theorem [22][Theorem 3.1], the sequences $\sigma$ such that $\operatorname{rank}H_{\sigma}=r<\infty$ are the polynomial-exponential series $\sigma\in\mathcal{{POLYEXP}}$ with ${\mu}(\sigma)=r$.

The aim of the method we are presenting, is to compute the structure of the Artinian Gorenstein algebra $\mathcal{A}_{\sigma}$ from the first terms of the sequence $\sigma=(\sigma_{\alpha})_{\alpha\in\mathbbm{N}^{n}}$. We are going to determine bases of $\mathcal{A}_{\sigma}$ and generators of the ideal $I_{\sigma}$, from which we can deduce directly the multiplicative structure of $\mathcal{A}_{\sigma}$.

The following lemma gives a simple way to test the linear independency in $\mathcal{A}_{\sigma}$ using truncated Hankel matrices (see [22][Lemma 3.3]):

Lemma 2.6

*Let $\sigma\in\mathbbm{K}[\boldsymbol{x}]^{\ast}$, $B=\{b_{1},\ldots,b_{r}\}$, $B^{\prime}=\{b_{1}^{\prime},\ldots,$ $b_{r}^{\prime}\}\subset\mathbbm{K}[\boldsymbol{x}]$. If the matrix $H_{\sigma}^{B,B^{\prime}}=(\langle\sigma|b_{i}b_{j}^{\prime}\rangle)_{\in B,\beta^{\prime}\in B^{\prime}}$ is invertible, then $B$ (resp. $B^{\prime}$) is linearly independent in $\mathcal{A}_{\sigma}$. *

This lemma implies that if $\dim\mathcal{A}_{\sigma}=r<+\infty$, $|B|=|B^{\prime}|=r=\dim\mathcal{A}_{\sigma}$ and $H_{\sigma}^{B,B^{\prime}}$ is invertible, then $(\boldsymbol{x}^{\beta})_{\beta\in B}$ and $(\boldsymbol{x}^{\beta^{\prime}})_{\beta^{\prime}\in B^{\prime}}$ are bases of $\mathcal{A}_{\sigma}$.

Given a Hankel operator $H_{\sigma}$ of finite rank $r$, it is clear that the truncated operators will have at most rank $r$. We are going to use the so-called flat extension property, which gives conditions under which a truncated Hankel operator of rank $r$ can be extended to a Hankel operator of the same rank (see [19] and extensions [10], [6], [22]).

Theorem 2.7

*Let $V,V^{\prime}\subset\mathbbm{K}[\boldsymbol{x}]$ be vector spaces connected to $1$, such that $x_{1},\ldots,x_{n}\in V$ and let $\sigma\in\langle V\cdot V^{\prime}\rangle^{\ast}$. Let $B\subset V$, $B^{\prime}\subset V^{\prime}$ such that $B^{+}\subset V,B^{\prime+}\subset V^{\prime}$. If $\operatorname{rank}H_{\sigma}^{V,V^{\prime}}=\operatorname{rank}H_{\sigma}^{B,B^{\prime}}=r$, then there is a unique extension $\tilde{\sigma}\in\mathbbm{K}[[\boldsymbol{y}]]$ such that $\tilde{\sigma}$ coincides with $\sigma$ on $\langle V\cdot V^{\prime}\rangle$ and $\operatorname{rank}H_{\tilde{\sigma}}=r$. In this case, $\tilde{\sigma}\in\mathcal{{POLYEXP}}$ with $r={\mu}(\sigma)$ and $I_{\tilde{\sigma}}=(\ker H_{\sigma}^{B^{+},B^{\prime}})$. *

2.4 Border bases

We recall briefly the definition of border basis and the main properties, that we will need. Let $B$ be a monomial set of $\mathbbm{K}[\boldsymbol{x}]$.

Definition 2.8

*A rewriting family $F$ for a (monomial) set $B$ is a set of polynomials $F=\{f_{i}\}_{i\in\boldsymbol{i}}\subset\mathbbm{K}[\boldsymbol{x}]$ such that $f_{i}=\boldsymbol{x}^{\alpha_{i}}+b_{i}$ with $b_{i}\in\langle B\rangle$, $\alpha_{i}\in\partial B$, $\alpha_{i}\neq\alpha_{j}$ if $i\neq j$. The rewriting family $f$ is complete if $(\boldsymbol{x}^{\alpha_{i}})_{i\in\boldsymbol{i}}=\partial B$. *

The monomial $\boldsymbol{x}^{\alpha_{i}}$ is called the leading monomial of $f_{i}$ and denoted $\gamma(f_{i}).$

Definition 2.9

*A family $F\subset\mathbbm{K}[\boldsymbol{x}]$ is a border basis with respect to $B$ if it is a complete rewriting family for $B$ such that $\mathbbm{K}[\boldsymbol{x}]=\langle B\rangle\oplus(F)$. *

This means that any element of $\mathbbm{K}[\boldsymbol{x}]$ can be projected along the ideal $I=(F)$ onto a unique element of $\langle B\rangle$. In other words, $B$ is a basis of the quotient algebra $\mathcal{A}=\mathbbm{K}[\boldsymbol{x}]/I$.

Let $B^{[0]}=B$ and for $k\in\mathbbm{N}$, $B^{[k+1]}=(B^{[k]})^{+}$. If $1\in B$, then for any $p\in\mathbbm{K}[\boldsymbol{x}]$, there exist $k\in\mathbbm{N}$, such that $p\in\langle B^{[k]}\rangle$.

For a complete rewriting family $F$ with respect to a monomial set $B$ containing $1$, a projection $\pi_{F}$ of $\mathbbm{K}[\boldsymbol{x}]$ on $\langle B\rangle$ can be defined recursively on the set of monomials $m$ of $\mathbbm{K}[\boldsymbol{x}]$ by

• •

if $m\in B$, $\pi_{F}(m)=m$;

• •

if $m\in\partial B$, $\pi_{F}(m)=m-f$ where $f$ is the (unique) polynomial in $F$ for which $\gamma(f)=m$,

• •

if $m\in B^{[k+1]}-B^{[k]}$ for $k>1$, there exists $m^{\prime}\in B^{[k]}$ and $i_{0}\in[1,n]$ such that $m=x_{i_{0}}m^{\prime}$. Let $\pi_{F}(m)=\pi_{F}(x_{i_{0}}\pi_{F}(m^{\prime}))$.

This map defines a projector from $\mathbbm{K}[\boldsymbol{x}]$ onto $\langle B\rangle$. The kernel of $\pi_{F}$ is contained in the ideal $(F)$. The family $F$ is a border basis iff $\ker(\pi_{F})=(F)$.

Checking that a complete rewriting family is a border basis reduces to checking commutation properties. This leads to efficient algorithms to compute a border basis. For more details, see [24], [25], [26].

A special case of border basis is when the leading term $\gamma(f)$ of $f\in F$ is the maximal monomial of $f$ for a monomial ordering $\succ$. Then $F$ is a Gröbner basis of $I$ for this monomial ordering $\succ$.

A border basis $F$ with respect to a monomial set $B$ gives directly the tables of multiplication $M_{i}$ by the variables $x_{i}$ in the basis $B$. For a monomial $b\in B$, $M_{i}(b)=\pi_{F}(x_{i}b)=x_{i}b-f$ with$f\in F$ such that $\gamma(f)=x_{i}b$ if $x_{i}b\in\partial B$ and $f=0$ otherwise.

3 Border bases of series

Given the first terms $\sigma_{\alpha}$ for $\alpha\in\boldsymbol{a}$ of the sequence $\sigma=(\sigma_{\alpha})_{\alpha\in\mathbbm{N}^{n}}\in\mathbbm{K}^{\mathbbm{N}^{n}}$, where $\boldsymbol{a}\subset\mathbbm{N}^{n}$ is a finite set of exponents, we are going to compute a basis of $\mathcal{A}_{\sigma}$ and generators of $I_{\sigma}$. We assume that the monomial set $\boldsymbol{x}^{\boldsymbol{a}}=\{\boldsymbol{x}^{\alpha},\alpha\in\boldsymbol{a}\}$ is connected to 1.

3.1 Orthogonal bases of $\mathcal{A}_{\sigma}$

An important step in the decomposition method consists in computing a basis $B$ of $\mathcal{A}_{\sigma}$. In this section, we describe how to compute a monomial basis $B=\{\boldsymbol{x}^{\beta}\}$ and two other bases $\boldsymbol{p}=(p_{\beta})$ and $\boldsymbol{q}=(q_{\beta})$, which are pairwise orthogonal for the inner product $\langle\cdot,\cdot\rangle_{\sigma}$:

[TABLE]

To compute these pairwise orthogonal bases, we will use a projection process, similar to Gram-Schmidt orthogonalization process. The main difference is that we compute pairs $p_{\beta},q_{\beta}$ of orthogonal polynomials. As the inner product $\langle\cdot,\cdot\rangle_{\sigma}$ may be isotropic, the two polynomials $p_{\beta},q_{\beta}$ may not be equal, up to a scalar. For a polynomial $f$ and two families of polynomials $\boldsymbol{p}=[p_{1},\ldots,p_{l}]$, $\boldsymbol{m}=[m_{1},\ldots,m_{l}]$, we will use the following procedure $\operatorname{proj}(f,\boldsymbol{p},\boldsymbol{m})$.

Algorithm 1 corresponds to the Modified Gram-Schmidt algorithm, when the scalar product is definite positive. It is known to have a better numerical behavior than the direct Gram-Schmidt orthogonalization process [34][Lecture 8]. It computes the polynomial $\operatorname{proj}(f,\boldsymbol{p},\boldsymbol{m})$ characterized by the following lemma.

Lemma 3.1

*If $\langle p_{i},m_{j}\rangle_{\sigma}=0$ if $j<i$  and $\langle p_{i},m_{i}\rangle_{\sigma}=1$, there is a unique polynomial $g$ such that $g=f-\sum_{i=1}^{l}\lambda_{i}p_{i}$ with $\lambda_{i}\in\mathbbm{K}$ and $\langle g,m_{i}\rangle_{\sigma}=0$ for $i=1,\ldots,l$. *

Proof 3.2

We prove by induction on the index $i$ of the loop that $g$ is orthogonal to $[m_{1},\ldots,m_{i}]$. For $i=1$, $g=f-\langle f,m_{1}\rangle_{\sigma}p_{1}$ is such that $\langle g,m_{1}\rangle_{\sigma}=\langle f,m_{1}\rangle_{\sigma}-\langle f,m_{1}\rangle_{\sigma}\langle p_{1},m_{1}\rangle_{\sigma}$ $=0$.

*If the property is true at step $k\leqslant l$, i.e. $\langle g,m_{i}\rangle_{\sigma}=0$ for $i<k$, then $g^{\prime}=g-\langle g,m_{k}\rangle_{\sigma}p_{k}$ is such that $\langle g,m_{i}\rangle_{\sigma}-\langle g,m_{k}\rangle_{\sigma}\langle p_{k},m_{i}\rangle_{\sigma}=\langle g,m_{i}\rangle_{\sigma}=0$ by induction hypothesis. By construction, $\langle g^{\prime},m_{k}\rangle=\langle g,m_{k}\rangle_{\sigma}-\langle g,m_{k}\rangle_{\sigma}\langle p_{k},m_{k}\rangle_{\sigma}=0$ and the induction hypothesis is true for $k$. As the matrix $(\langle p_{\operatorname{ij}},m_{i}\rangle_{\sigma})_{1\leqslant i,j\leqslant l}$ is invertible, there exists a unique polynomial of the form $g=f-\sum_{j=1}^{l}\lambda_{j}p_{j}$, such that $\langle g,m_{i}\rangle_{\sigma}=0$ for $i=1,\ldots,l$, which concludes the proof of the lemma. *

Algorithm 2 for computing a border basis of $\mathcal{A}_{\sigma}$ proceeds inductively, starting from $\boldsymbol{p}=[],\boldsymbol{m}=[],\boldsymbol{b}=[]$, extending the basis $\boldsymbol{p}$ with a new polynomial $p_{\alpha}$, orthogonal to the vector space spanned by $\boldsymbol{m}$ for the inner product $\langle\cdot,\cdot\rangle_{\sigma}$, extending $\boldsymbol{m}$ with a new monomial $m_{\alpha}$, such that $\langle p_{\alpha},m_{\alpha}\rangle_{\sigma}=1$ and $\langle p_{\beta},m_{\alpha}\rangle=0$ for $\beta\in\boldsymbol{b}$ and extending $\boldsymbol{b}$ with $\alpha$.

The main difference with Algorithm 4.1 in [22] is the projection procedure and the list of monomials $\boldsymbol{s}$, $\boldsymbol{t}$ used to generate new monomials and to perform the projections. The lists $\boldsymbol{b},\boldsymbol{d},\boldsymbol{c},\boldsymbol{s},\boldsymbol{t}$ are lists of exponents, identified with monomials.

We verify that at each loop of the algorithm, the lists $\boldsymbol{b}$, $\boldsymbol{d}$ and $\boldsymbol{s}$ are disjoint and $\boldsymbol{b}\cup\boldsymbol{d}\cup\boldsymbol{s}=\boldsymbol{a}$. We also verify that $m_{\alpha}$ are monomials up to a scalar, that the set of their exponents is $\boldsymbol{c}$, that $\boldsymbol{c}$ and $\boldsymbol{t}$ are disjoint and that $\boldsymbol{c}\cup\boldsymbol{t}=\boldsymbol{a}$.

The algorithm uses the function $\operatorname{next}(\boldsymbol{b},\boldsymbol{d},\mathbf{c},\boldsymbol{s})$, which computes the set of monomials $\mathfrak{n}$ in $\partial\boldsymbol{b}\cap\boldsymbol{s}$, which are not in $\boldsymbol{d}$ and such $\mathfrak{n}\cdot\mathbf{c}\subset\mathbf{a}=\boldsymbol{b}\cup\boldsymbol{d}\cup\boldsymbol{s}$.

We denote by $\prec$ the order induced by the treatment of the monomials of $\boldsymbol{a}$ in the loops of the algorithm, so that the monomials treated at the $l^{\operatorname{th}}$ loop are smaller than the monomials in $\boldsymbol{n}$ at the $(l+1)^{\operatorname{th}}$ loop. For $\alpha\in\boldsymbol{a}$, we denote by $\boldsymbol{b}_{\prec\alpha}$ the list of monomial exponents $\beta\in\boldsymbol{b}$ with $\beta\prec\alpha$ and by $B_{\prec\alpha}$ the vector space spanned by these monomials. For $\alpha\in\boldsymbol{b}$, let $\boldsymbol{b}_{\preccurlyeq\alpha}=\boldsymbol{b}_{\prec\alpha}\cup[\alpha]$.

The following properties are also satisfied during this algorithm:

Lemma 3.3

*For $\alpha\in\boldsymbol{b}$, we have $\forall\beta\in\boldsymbol{b}_{\prec\alpha}$, $\langle p_{\alpha},m_{\beta}\rangle_{\sigma}=0$ and $\langle p_{\alpha},m_{\alpha}\rangle_{\sigma}=1$. For $\alpha\in\boldsymbol{d}$, $\langle p_{\alpha},\boldsymbol{x}^{\gamma}\rangle_{\sigma}=0$ for all $\gamma\in\boldsymbol{a}$ such that $\boldsymbol{x}^{\gamma}p_{\alpha}\in\langle\boldsymbol{x}^{\boldsymbol{a}}\rangle$. *

Proof 3.4

By construction,

[TABLE]

is orthogonal to $m_{\beta}$ for $\beta\in\boldsymbol{b}_{\prec\alpha}$. We consider two exclusive cases: $\alpha\in\boldsymbol{b}$ and $\alpha\in\boldsymbol{d}$.

• •

If $\alpha\in\boldsymbol{b}$, then there exists $\boldsymbol{x}^{\gamma}\in\boldsymbol{s}$ such that $\langle p_{\alpha},\boldsymbol{x}^{\gamma}\rangle_{\sigma}\neq 0$. Thus $m_{\alpha}=\frac{1}{\langle p_{\alpha},\boldsymbol{x}^{\gamma}\rangle_{\sigma}}\boldsymbol{x}^{\gamma}$ is such that $\langle p_{\alpha},m_{\alpha}\rangle_{\sigma}=1$. By construction, $\langle p_{\alpha},m_{\beta}\rangle_{\sigma}=0$ for $\beta\in\boldsymbol{b}_{\prec\alpha}$.

• •

If $\alpha\in\boldsymbol{d}$, then there is no $\boldsymbol{x}^{\gamma}\in\boldsymbol{s}$ such that $\langle p_{\alpha},\boldsymbol{x}^{\gamma}\rangle_{\sigma}\neq 0$ and $\boldsymbol{x}^{\gamma}p_{\alpha}\in\langle\boldsymbol{x}^{\boldsymbol{a}}\rangle$. Thus $p_{\alpha}$ is orthogonal to $\boldsymbol{x}^{\gamma}$ for all $\gamma\in\boldsymbol{s}$ with $\boldsymbol{x}^{\gamma}p_{\alpha}\in\langle\boldsymbol{x}^{\boldsymbol{a}}\rangle$. By construction, $p_{\alpha}$ is orthogonal to $m_{\beta}$ for $\beta\in\boldsymbol{b}$. As $\boldsymbol{b}\cup\boldsymbol{s}=\boldsymbol{a}$,  $\langle p_{\alpha},\boldsymbol{x}^{\gamma}\rangle_{\sigma}=0$ for all $\gamma\in\boldsymbol{a}$ such that $\boldsymbol{x}^{\gamma}p_{\alpha}\in\langle\boldsymbol{x}^{\boldsymbol{a}}\rangle$.

*This concludes the proof of this lemma. *

Lemma 3.5

*For $\alpha\in\boldsymbol{b}$, $\langle m_{\beta}\rangle_{\beta\in\boldsymbol{b}_{\preccurlyeq\alpha}}=\langle q_{\beta}\rangle_{\beta\in\boldsymbol{b}_{\preccurlyeq\alpha}}$ and the bases $\boldsymbol{p}=[p_{\beta}]_{\beta\in\boldsymbol{b}_{\preccurlyeq\alpha}},\boldsymbol{q}=[q_{\beta}]_{\beta\in\boldsymbol{b}_{\preccurlyeq\alpha}}$ are pairwise orthogonal. *

Proof 3.6

We prove it by induction on $\alpha$. If $\alpha=\boldsymbol{0}$ is not in b, then $\sigma_{\alpha}=0$ for all $\alpha\in\boldsymbol{a}$, $\boldsymbol{p}$ and $\boldsymbol{q}$ are empty and the property is satisfied. If $\alpha=\boldsymbol{0}$ is in $\boldsymbol{b}$, then $p_{\alpha}=1$ and $q_{\alpha}=m_{\alpha}$ is such that $\langle p_{\alpha},m_{\alpha}\rangle_{\sigma}=1$. The property is true for $\alpha=\boldsymbol{0}$.

Suppose that it is true for all $\beta\in\boldsymbol{b}_{\prec\alpha}$. By construction, the polynomial $q_{\alpha}=\operatorname{proj}(m_{\alpha},[q_{{}_{\beta}}]_{\beta\in\boldsymbol{b}_{\prec\alpha}},$ $[p_{\beta}]_{\beta\in\boldsymbol{b}_{\prec\alpha}})$ is orthogonal to $p_{\beta}$ for $\beta\prec\alpha$. By induction hypothesis, $[p_{\beta}]_{\beta\in\boldsymbol{b}_{\prec\alpha}},\boldsymbol{q}=[q_{\beta}]_{\beta\in\boldsymbol{b}_{\prec\alpha}}$ are pairwise orthogonal, thus

[TABLE]

By the induction hypothesis, we deduce that

[TABLE]

By Lemma 3.3, $p_{\alpha}$ is orthogonal to $m_{\beta}$ for $\beta\in\boldsymbol{b}_{\prec\alpha}$ and thus to $q_{\beta}$ for $\beta\in\boldsymbol{b}_{\prec\alpha}$. We deduce that

[TABLE]

*This shows that $[p_{\beta}]_{\beta\in\boldsymbol{b}_{\preccurlyeq\alpha}}$ and $\boldsymbol{q}=[q_{\beta}]_{\beta\in\boldsymbol{b}_{\preccurlyeq\alpha}}$ are pairwise orthogonal and concludes the proof by induction. *

Lemma 3.7

*At the $l^{\operatorname{th}}$ loop of the algorithm, the polynomials $p_{\alpha}$ for $\alpha\in\boldsymbol{n}$ are of the form $p_{\alpha}=\boldsymbol{x}^{\alpha}+b_{\alpha}$ with $b_{\alpha}\in B_{\prec\alpha}$. *

Proof 3.8

We prove by induction on the loop index $l$ that we have  $p_{\alpha}=\boldsymbol{x}^{\alpha}+b_{\alpha}$ with $b_{\alpha}\in B_{\prec\alpha}$.

The property is clearly true for $l=0$, $\alpha=\boldsymbol{0}$ and $p_{\alpha}=1=\boldsymbol{x}^{\boldsymbol{0}}$. Suppose that it is true for any $l^{\prime}<l$ and consider the $l^{\operatorname{th}}$ loop of the algorithm. The polynomial $p_{\alpha}$ is constructed by projection of $\boldsymbol{x}^{\alpha}$ on $\langle p_{\alpha}\rangle_{\beta\in\boldsymbol{b}}$ orthogonally to $\langle m_{\beta}\rangle_{\beta\in\boldsymbol{b}}$ where $\boldsymbol{b}=\boldsymbol{b}_{\prec\alpha}$. By induction hypothesis, $p_{\beta}=\boldsymbol{x}^{\beta}+b_{\beta}$ with $b_{\beta}\in B_{\prec\beta}\subset B_{\prec\alpha}$. Then by Lemma 3.1, we have

[TABLE]

*with $\lambda_{\beta}\in\mathbbm{K}$, $b_{\alpha}\in B_{\prec\alpha}$. Thus, the induction hypothesis is true for $l$, which concludes the proof. *

3.2 Quotient algebra structure

We show now that the algorithm outputs a border basis of an Artinian Gorenstein algebra $\mathcal{A}_{\tilde{\sigma}}$ for an extension $\tilde{\sigma}$ of $\sigma$, when all the border relations are computed, that is, when $\boldsymbol{d}=\partial\boldsymbol{b}$.

Theorem 3.9

Let $\boldsymbol{b}=[\beta_{1},\ldots,\beta_{r}]$, $\boldsymbol{c}=[\gamma_{1},\ldots,\gamma_{r}]$, $\boldsymbol{p}=[p_{\beta_{1}},\ldots,p_{\beta_{r}}]$, $\boldsymbol{q}=[q_{\beta_{1}},\ldots,q_{\beta_{r}}]$ and $\boldsymbol{k}=[p_{\alpha_{1}},\ldots,p_{\alpha_{s}}]$ be the output of Algorithm 2. Let $V=\langle\boldsymbol{x}^{\boldsymbol{b}^{+}}\rangle$. If $\boldsymbol{d}=\partial\boldsymbol{b}$ and $\boldsymbol{c}^{+}\subset\boldsymbol{b}^{\prime}$ connected to $1$ such that $\boldsymbol{x}^{\boldsymbol{b}^{+}}\cdot\boldsymbol{x}^{\boldsymbol{b}^{\prime}}=\boldsymbol{x}^{\boldsymbol{a}}$ then $\sigma$ coincides on $\langle\boldsymbol{x}^{\boldsymbol{a}}\rangle$ with a series $\tilde{\sigma}\in\mathbbm{K}[[\boldsymbol{y}]]$ such that

• •

$\operatorname{rank}H_{\tilde{\sigma}}=r$,

• •

$(\boldsymbol{p},\boldsymbol{q})$* are pairwise orthogonal bases of $\mathcal{A}_{\tilde{\sigma}}$ for the inner product $\langle\cdot,\cdot\rangle_{\tilde{\sigma}}$,*

• •

The family $\boldsymbol{k}=\{p_{\alpha},\alpha\in\partial\boldsymbol{b}\}$ is a border basis of the ideal $I_{\tilde{\sigma}}$, with respect to $\boldsymbol{x}^{\boldsymbol{b}}$.

• •

The matrix of multiplication by $x_{k}$ in the basis $\boldsymbol{p}$ (resp. q) of $\mathcal{A}_{\tilde{\sigma}}$ is $M_{k}:=(\langle\sigma|x_{k}p_{\beta_{j}}q_{\beta_{i}}\rangle)_{1\leqslant i,j\leqslant r}$ (resp. $M_{k}^{t}$).

Proof 3.10

By construction, $\boldsymbol{x}^{\boldsymbol{b}^{+}}$ is connected to $1$. Let $V=\langle\boldsymbol{x}^{\boldsymbol{b}^{+}}\rangle$ and $V^{\prime}=\langle\boldsymbol{x}^{\boldsymbol{b}^{\prime}}\rangle$. As $\boldsymbol{b}^{+}=\boldsymbol{b}\cup\boldsymbol{d}$, a basis of $V$ is formed by the monomials $\boldsymbol{x}^{\boldsymbol{b}}$ and the polynomials $p_{\alpha}=\boldsymbol{x}^{\alpha}+b_{\alpha}$ with $b_{\alpha}\in\langle\boldsymbol{x}^{\boldsymbol{b}}\rangle$ for $\alpha\in\boldsymbol{d}$. The matrix of $H_{\sigma}^{V,V^{\prime}}$ in this basis of $V$ and a basis of $V^{\prime},$ which first elements are $m_{\beta_{1}},\ldots,m_{\beta_{r}}$, is of the form

[TABLE]

where $L_{r}$ is a lower triangular invertible matrix of size $r$. The kernel of $H_{\sigma}^{V,V^{\prime}}$ is generated by the polynomials $p_{\alpha}$ for $\alpha\in\boldsymbol{d}$.

By Theorem 2.7, $\sigma$ coincides on $V\cdot V^{\prime}=\langle\boldsymbol{x}^{\boldsymbol{a}}\rangle$ with a series $\tilde{\sigma}$ such that $\boldsymbol{x}^{\boldsymbol{b}}$ is a basis of $\mathcal{A}_{\bar{\sigma}}=\mathbbm{K}[\boldsymbol{x}]/I_{\tilde{\sigma}}$ and $I_{\tilde{\sigma}}=(\ker H_{\tilde{\sigma}}^{V,V^{\prime}})=(p_{\alpha})_{\alpha\in\boldsymbol{d}}$.

By Lemma 3.7, $p_{\alpha}=\boldsymbol{x}^{\alpha}+b_{\alpha}$ with $\alpha\in\partial\boldsymbol{b}$ and $b_{\alpha}\in\langle\boldsymbol{x}^{\boldsymbol{b}}\rangle$. Thus $(p_{\alpha})_{\alpha\in\partial\boldsymbol{b}}$ is a border basis with respect to $\boldsymbol{x}^{\boldsymbol{b}}$ for the ideal $I_{\tilde{\sigma}}$, since $\boldsymbol{x}^{\boldsymbol{b}}$ is a basis of of $\mathcal{A}_{\bar{\sigma}}$. This shows that $\operatorname{rank}H_{\tilde{\sigma}}=\dim\mathcal{A}_{\tilde{\sigma}}=|\boldsymbol{b}|=r$.

By Lemma 3.5, $(\boldsymbol{p},\boldsymbol{q})$ are pairwise orthogonal for the inner product $\langle\cdot,\cdot\rangle_{\sigma}$, which coincides with $\langle\cdot,\cdot\rangle_{\tilde{\sigma}}$ on $\langle\boldsymbol{x}^{\boldsymbol{a}}\rangle$. Thus they are pairwise orthogonal bases of $\mathcal{A}_{\tilde{\sigma}}$ for the inner product $\langle\cdot,\cdot\rangle_{\tilde{\sigma}}$.

As we have $x_{k}p_{\beta_{j}}\equiv\sum_{i=1}^{r}\langle x_{k}p_{\beta_{j}},q_{\beta_{i}}\rangle_{\sigma}p_{\beta_{i}}$, the matrix of multiplication by $x_{k}$ in the basis $\boldsymbol{p}$ of $\mathcal{A}_{\tilde{\sigma}}$ is

[TABLE]

*Exchanging the role of p and q, we obtain $M_{k}^{t}$ for the matrix of multiplication by $x_{k}$ in the basis q. *

Lemma 3.11

*If $\prec$ is a monomial ordering and if at the end of the algorithm $\boldsymbol{d}=\partial\boldsymbol{b}$ and $\boldsymbol{c}^{+}\subset\boldsymbol{b}^{\prime}$ connected to $1$ with $\boldsymbol{x}^{\boldsymbol{b}^{+}}\cdot\boldsymbol{x}^{\boldsymbol{b}^{\prime}}=\boldsymbol{x}^{\boldsymbol{a}}$, then $\boldsymbol{b}=\boldsymbol{c}$ and $\boldsymbol{k}$ is a Gröbner basis of the ideal $I_{\sigma}$ for the monomial ordering. *

Proof 3.12

If $\prec$ is a monomial ordering, then the polynomials $p_{\alpha}=\boldsymbol{x}^{\alpha}+b_{\alpha}$, $\alpha\in\partial\boldsymbol{b}$ are constructed in such a way that their leading term is $\boldsymbol{x}^{\alpha}$. Therefore the border basis $\boldsymbol{k}=(p_{\alpha})_{\alpha\in\partial\boldsymbol{b}}$ is also a Gröbner basis.

By construction, $\boldsymbol{c}$ is the set of monomials $\gamma\in\boldsymbol{a}$ such that $\langle p_{\beta},\boldsymbol{x}^{\gamma}\rangle_{\sigma}\neq 0$ for some $\beta\in\boldsymbol{b}$. Suppose that $\gamma\in\boldsymbol{c}$ is not in $\boldsymbol{b}$. Then $\boldsymbol{x}^{\gamma}\in(\boldsymbol{x}^{\boldsymbol{d}})$ and there is $\delta\in\boldsymbol{d}$ and $\gamma^{\prime}\in\boldsymbol{a}$ such that $\gamma=\delta+\gamma^{\prime}$. As $p_{\delta}\in\boldsymbol{k}$, we have $\langle p_{\delta},\boldsymbol{x}^{\alpha}\rangle_{\sigma}=0$ for $\alpha\in\boldsymbol{a}$ such that $p_{\delta}\boldsymbol{x}^{\alpha}\in\langle\boldsymbol{x}^{\boldsymbol{a}}\rangle$. By Lemma 3.7, $p_{\delta}=\boldsymbol{x}^{\delta}+b_{\delta}$ with $b_{\delta}\in\boldsymbol{b}_{\prec\delta}$ with $\boldsymbol{x}^{\delta}\succ b_{\delta}$.

[TABLE]

*We have $\langle p_{\beta},p_{\delta}\boldsymbol{x}^{\gamma^{\prime}}\rangle_{\sigma}=\langle p_{\delta},p_{\beta}\boldsymbol{x}^{\gamma^{\prime}}\rangle_{\sigma}=0$ since $p_{\delta}\in\boldsymbol{k}$ and $p_{\delta}p_{\beta}\boldsymbol{x}^{\gamma^{\prime}}\in\langle\boldsymbol{x}^{\boldsymbol{a}}\rangle$. As $\gamma$ is the first monomial of $\boldsymbol{a}$ such that $\langle p_{\beta},\boldsymbol{x}^{\gamma}\rangle_{\sigma}\neq 0$ and $b_{\delta}\boldsymbol{x}^{\gamma^{\prime}}\prec\boldsymbol{x}^{\delta+\gamma^{\prime}}=\boldsymbol{x}^{\gamma}$, we have $\langle p_{\beta},b_{\delta}\boldsymbol{x}^{\gamma^{\prime}}\rangle_{\sigma}$, which implies that $\langle p_{\beta},\boldsymbol{x}^{\gamma}\rangle_{\sigma}=0$. This is in contradiction with the hypothesis $\langle p_{\beta},\boldsymbol{x}^{\gamma}\rangle_{\sigma}\neq 0$, therefore $\gamma\in\boldsymbol{b}$. We deduce that $\boldsymbol{c}\subset\boldsymbol{b}$ and the equality holds since the two sets have the same cardinality. *

Notice that to construct a minimal reduced Gröbner basis of $I_{\tilde{\sigma}}$ for the monomial ordering $\prec$, it suffices to keep the elements $p_{\alpha}\in\boldsymbol{k}$ with $\alpha$ minimal for the component-wise partial ordering.

3.3 Complexity

Let $s=|\boldsymbol{a}|$ and $r=|\boldsymbol{b}|$, $\delta=|\partial\boldsymbol{b}|$. As $\boldsymbol{b}\subset\boldsymbol{a}$ and the monomials in $\partial\boldsymbol{b}$ are the product of a monomial in $\boldsymbol{b}$ by one of the variables $x_{1},\ldots,x_{n}$, we have $r\leqslant s$ and $\delta\leqslant nr$.

Proposition 3.13

*The complexity of the algorithm to compute the bases $\boldsymbol{p}$ and $\boldsymbol{q}$ is $\mathcal{O}((r+\delta)rs)$. *

Proof 3.14

*At each step, the computation of $p_{\alpha}$ (resp. $q_{\alpha}$) requires $\mathcal{O}(r^{2})$ arithmetic operations, since the support of the polynomials $p_{\beta}$, $q_{\beta}$ $(\beta\in\boldsymbol{b})$ is in $\boldsymbol{b}$ and $|\boldsymbol{b}|\leqslant r$. Computing $\langle\boldsymbol{x}^{\gamma},p_{\alpha}\rangle_{\sigma}$ for all $\gamma\in\boldsymbol{t}$ requires $\mathcal{O}(rs)$ arithmetic operations. As the number of polynomials $p_{\alpha}$ is at most $|\boldsymbol{b}^{+}|=r+\delta$, the total cost for computing $\boldsymbol{p}$ and $\boldsymbol{q}$ is thus in $\mathcal{O}((r+\delta)(r^{2}+rs))=\mathcal{O}((r+\delta)rs)$. *

As $\delta\leqslant nr$, the complexity of this algorithm is in $\mathcal{O}(nr^{2}s)$.

The algorithm is connected to the Berlekamp-Massey-Sakata algorithm, which computes a Gröbner basis for a monomial ordering $\prec$. In the BMS algorithm, a minimal set $\mathcal{F}$ of recurrence polynomials valid for the monomials smaller that a given monomial $m$ is computed. A monomial basis $\boldsymbol{b}^{\ast}$ generated by all the divisors of some corner elements is constructed. The successor $m^{+}$ of the monomial $m$ for the monomial ordering $\prec$ is considered and the family $\mathcal{F}$ of valid recurrence polynomials is updated by computing their discrepancy at the monomial $m^{+}$ and by cancelling this discrepancy, if necessary, by combination with one lower polynomial [30].

Let $\delta$ be the size of the border $\partial\boldsymbol{b}^{\ast}$ of the monomial basis $\boldsymbol{b}^{\ast}$computed by BMS algorithm. At each update, there are at most $\delta$ polynomials in $\mathcal{F}$. Let $s^{\prime}$ be the maximum number of their non-zero terms. Then the update of $\mathcal{F}$ requires $\mathcal{O}(\delta s^{\prime})$ arithmetic operations. The number of updates is bounded by the number $r+\delta$ of monomials in $\boldsymbol{b}^{+}$. Checking the discrepency of a polynomial in $\mathcal{F}$ for all the monomials in $\boldsymbol{x}^{\boldsymbol{a}}$ requires $\mathcal{O}(s^{\prime}s)$ arithmetic operations. Thus, the total cost of the BMS algorithm is in $\mathcal{O}((r+\delta)\delta s^{\prime}+\delta s^{\prime}s)$. As the output polynomials in the Gröbner basis are not necessarily reduced, the maximal number of terms $s^{\prime}\leqslant s$ can be of the same order than $s$. Thus the complete complexity of BMS algorithm is in $\mathcal{O}(\delta s^{2})=\mathcal{O}(nrs^{2})$, which is an order of magnitude larger than the bound of Proposition 3.13, assuming that $r\ll s$.

The method presented in [7] for computing a Gröbner basis of the recurrence polynomials computes the rank of a Hankel matrix of size the number $\tilde{s}$ of monomials of degree $\leqslant d$ for a bound $d$ on the degree of the recurrence relations. It deduces a monomial basis $\boldsymbol{b}$ stable by division and obtains the valid recurrence relations for the border monomials by solving a linear Hankel system of size $r$. Thus the complexity is in $\mathcal{O}(\delta r^{\omega}+\tilde{s}^{\omega})$ where $2.3\leqslant\omega\leqslant 3$. It is also larger than the bound of Proposition 3.13. This bound could be improved by exploiting the rank displacement of the structured matrices involved in this method [9], but the known bounds on the displacement rank of the matrices involved in the computation do not improve the bound of Proposition 3.13.

4 Examples

4.1 Multivariate Prony method

Given a function $h(u_{1},\ldots,u_{n})=\sum_{i=1}^{r}\omega_{i}e^{\zeta_{i,1}u_{1}+\cdots+\zeta_{i,n}u_{n}}$, the problem is to compute its decomposition as a weighted sum of exponentials, from values of $h$. The method proposed by G. Riche de Prony for sums of univariate exponential functions consists in sampling the function at regularly spaced values [2]. In the multivariate extension of this method, the function is sampled on a grid in $\mathbbm{R}^{n}$, for instance $\mathbbm{N}^{n}$. The decomposition is computed from a subset of the multi-index sequence of evaluation $\sigma_{\alpha}=h(\alpha_{1},\ldots,\alpha_{n})$ for $\alpha=(\alpha_{1},\ldots,\alpha_{n})\in\mathbbm{N}^{n}$. The ideal $I_{\sigma}$ associated to this sequence is the ideal defining the points $\xi_{i}=(e^{\zeta_{i,1}},\ldots,e^{\zeta_{i,1n}})$. To compute this decomposition, we apply the border basis algorithm to the sequence $\sigma_{\alpha}$ for $|\alpha|\leqslant d$ with $d$ high enough, and obtain a border basis of the ideal $I_{\sigma}$ defining the points $\xi_{1},\ldots,\xi_{r}\in\mathbbm{K}^{n}$, a basis of $\mathcal{A}_{\sigma}$ and the tables of multiplication in this basis. By applying the decomposition algorithm in [22], we deduce the points $\xi_{i}=(e^{\zeta_{i,1}},\ldots,e^{\zeta_{i,1n}})$. Taking the log of their coordinates $\log(\xi_{i,j})=\zeta_{i,j}$ yields the coordinates of the frequencies $\zeta_{i}$.

4.2 Fast decoding of algebraic-geometric codes

Let $\mathbbm{K}$ be a finite field. We consider an algebraic-geometric code $C$ obtained by evaluation of polynomials in $\mathbbm{K}[x_{1},\ldots,x_{n}]\overset{}{}$ of degree $\leqslant d$ at points $\xi_{1},\ldots,\xi_{l}\in\mathbbm{K}^{n}$. It is a finite vector space in $\mathbbm{K}^{l}$. We use the words of the orthogonal code $C^{\perp}=\{(m_{1},\ldots,m_{l})\mid m\cdot c=m_{1}c_{1}+\cdots+m_{l}c_{l}=0\}$ for the transmission of information. Suppose that an error $\omega=(\omega_{1},\ldots,\omega_{l})$ occurs in the transmission of a message $m=(m_{1},\ldots,m_{l})$ so that the message $m^{\ast}=m+\omega$ is received. Let $\omega_{i_{1}},\ldots,\omega_{i_{r}}$ be the non-zero coefficients of the error vector $\omega$. To correct the message $r$, we use the moments or syndromes $\sigma_{\alpha}=(\xi_{1}^{\alpha},\ldots,\xi_{l}^{\alpha})\cdot m^{\ast}=(\xi_{1}^{\alpha},\ldots,\xi_{l}^{\alpha})\cdot\omega=\sum_{j=1}^{r}w_{i_{j}}\xi_{i_{j}}^{\alpha}$ for $\alpha=(\alpha_{1},\ldots,\alpha_{n})\in\mathbbm{N}^{n}$ with $|\alpha|\leqslant d$. We compute generators of the set of error-locator polynomials, that is, the polynomials vanishing at the points $\xi_{i_{1}},\ldots,\xi_{i_{r}}$ and deduce the weights or errors $\omega_{i_{j}}$ by solving the Vandermonde system

[TABLE]

The points $\xi_{i_{j}}$ correspond to the position of the errors and $\omega_{i_{j}}$ to their amplitude. By applying the border basis algorithm, we obtain a border basis of  the ideal of error-locator polynomials, from which we deduce the position and amplitude of the errors.

4.3 Sparse interpolation

Given a sparse polynomial $h(u_{1},\ldots,u_{n})=\sum_{i=1}^{r}\omega_{i}u_{1}^{\gamma_{i,1}}$ $\cdots$ $u_{n}^{\gamma_{i,n}}$, which is a weighted sum of $r$ monomials with non-zero weights $\omega_{i}\in\mathbbm{K}$, the problem is to compute the exponents $(\gamma_{i,1},\ldots,\gamma_{i,n})\in\mathbbm{N}^{n}$ of the monomials and the weights $\omega_{i}$, from evaluations of the blackbox functions $h$.  The approach, proposed initially in [4], [36], consists in evaluating the function at points of the form $(\zeta_{1}^{k},\ldots,\zeta_{n}^{k})$ for some values of $\zeta_{1},\ldots\zeta_{n}\in\mathbbm{K}$ and to apply univariate Prony-type methods or Berlekamp-Massey algorithms to the sequence $\sigma_{k}=h(\zeta_{1}^{k},\ldots,\zeta_{n}^{k})$, for $k\in\mathbbm{N}$. The approach can be extended to multi-index sequences $(\sigma_{\alpha})_{\alpha\in\mathbbm{N}^{n}}$ by computing the terms

[TABLE]

for $\alpha=(\alpha_{1},\ldots,\alpha_{n})\in\mathbbm{N}^{n}$.  It can also be extended to sequences constructed from polylog functions [22]. By applying the border basis algorithm to the multi-index sequence $\sigma_{\alpha}=h(\zeta_{1}^{\alpha_{1}},\ldots,\zeta_{n}^{\alpha_{n}})$ for $|\alpha|\leqslant d$ with $d\in\mathbbm{N}$ high enough, we obtain generators of the ideal $I_{\sigma}$ defining the points $\xi_{i}=(\zeta_{1}^{\gamma_{i,1}},\ldots,\zeta_{n}^{\gamma_{i,n}})$ and deduce the weights $\omega_{i}$, $i=1,\ldots,r$. By computing the log of the coordinates of the points $\xi_{i}$, we deduce the exponent vectors $\gamma_{i}=(\gamma_{i,1},\ldots,\gamma_{i,n})\in\mathbbm{N}^{n}$ for $i=1,\ldots,r$.

4.4 Tensor decomposition

Given a homogeneous polynomial

[TABLE]

of degree $d\in\mathbbm{N}$ with $t_{\alpha}\in\mathbbm{K}$, $\binom{d}{\alpha}=\frac{d!}{\alpha_{0}!\cdots\alpha_{n}!}$, we want to a decomposition of $t$ as sum of powers of linear forms:

[TABLE]

with a minimal $r$, $\omega_{i}\neq 0$ and $(\xi_{i,0},\ldots,\xi_{i,n})\neq 0$. By a change of variables, we can assume that $\xi_{i,0}\neq 0$ in such a decomposition, and by dividing each linear form by $\xi_{i,0}$ and multiplying $\omega_{i}$ by $\xi_{i,0}^{d}$, we can even assume that $\xi_{i,0}=1$. Then by expansion of the powers of the linear forms and by identification of the coefficients, we obtain

[TABLE]

for $\alpha=(\alpha_{1},\ldots,\alpha_{n})\in\mathbbm{N}^{n}$ with $|\alpha|\leqslant d$. We apply the border basis algorithm to this sequence, in order to obtain generators of the ideal $I_{\sigma}$ defining the points $\xi_{1},\ldots,\xi_{r}\in\mathbbm{K}^{n}$ and providing the weights $\omega_{i}$. If the number of terms $r$ is small enough compared to the number of terms $\sigma_{\alpha}$, then the set of border relations are complete and it is possible to compute the decomposition (2).

4.5 Vanishing ideal of points

Given a set of points $\Xi=\{\xi_{1},\ldots,\xi_{r}\}\subset\mathbbm{K}^{n}$, we want to compute polynomials defining these points, that is, a set of generators of the ideal of polynomials vanishing on $\Xi$. For that purpose, we choose non-zero weights $w_{i}\in\mathbbm{K}$, a degree $d\in\mathbbm{N}$ and we compute the sequence of moments $\sigma_{\alpha}=\sum_{i=1}^{r}\omega_{i}\xi^{\alpha}$ for $|\alpha|\leqslant d$. The generating series $\sigma$ associated to these moments define an Artinian Gorenstein algebra $\mathcal{A}_{\sigma}=\mathbbm{K}[\boldsymbol{x}]/I_{\sigma}$, where $I_{\sigma}$ is the ideal of polynomials vanishing $\Xi$ [22]. This ideal $I_{\sigma}$ defines the points $\xi_{i}$  with multiplicity $1$. The idempotents $\{\boldsymbol{u}_{i}\}_{i=1\ldots r}$ associated to the points $\Xi$ form a family of interpolation polynomials at these points: $\boldsymbol{u}_{i}(\xi_{j})=0$ if $i\neq j$ and $\boldsymbol{u}_{i}(\xi_{i})=1$. They are the common eigenvectors of the multiplication operators in $\mathcal{A}_{\sigma}$. By applying the border basis algorithm to the sequence $\sigma_{\alpha}$ for $|\alpha|\leqslant d$ with $d$ high enough, we obtain generators of the ideal $I_{\sigma}$ defining the points $\xi_{1},\ldots,\xi_{r}\in\mathbbm{K}^{n}$, a basis of $\mathcal{A}_{\sigma}$ and the tables of multiplication in this basis. By computing the eigenvectors of a generic combination of the multiplication tables by a variable, we obtain a family of interpolation polynomials at the roots $\Xi$.

4.6 Benchmarks

We present some experimentations of an implementation of Algorithm 2111available at https://gitlab.inria.fr/mourrain/PolyExp in the programming language Julia222https://julialang.org/. The arithmetic operations are done in the finite field $\mathbbm{Z}/32003\mathbbm{Z}$. We choose $r$ random points $\xi_{i}$ with $n$ coordinates in $\mathbbm{Z}/32003\mathbbm{Z}$, take the sequence of moments $\sigma_{\alpha}=\sum_{i=1}^{r}\xi_{i}^{\alpha}$ up for $|\alpha|\leqslant d$  with weights equal to $1$. Figure 4.6 shows the timing (in sec.) to compute the border basis, checking the validity of the recurrence relations up to degree $d$. The computation is done on a MacOS El Capitan, 2.8 GHz Intel Core i7, 16 Go.

Fig. 4.6: Vanishing ideal of random points.

The timing is approximately linear in the number $r$ of points, with a slope increasing quadratically in $n$.

Examples

Example .1

We consider the sequence $\sigma\in\mathbbm{K}^{\mathbbm{N}}$ such that $\sigma_{d_{1}}=1$ and $\sigma_{i}=0$ for $0\leqslant i\neq d_{1}\leqslant d$ and $d_{1}<d$.

In the first step of the algorithm, we take $p_{0}=1$ and compute the first $\gamma\in[0,\ldots,d]$ such that $\langle x^{\gamma},p_{1}\rangle_{\sigma}$ is not zero. This yields $m_{0}=x^{d_{1}}$ and $\boldsymbol{b}=[0]$, $\boldsymbol{c}=[d_{1}]$.

In a second step, we have $p_{1}=x-\langle x,m_{1}\rangle_{\sigma}p_{0}=x$. The first $\gamma\in[0,\ldots,d]\setminus\{d_{1}\}$ such that $\langle x^{i},p_{1}\rangle_{\sigma}$ is not zero yields $\boldsymbol{b}=[0,1]$, $\boldsymbol{c}=[d_{1},d_{1}-1]$, $m_{1}=x^{d_{1}-1}$.

We repeat this computation until $\boldsymbol{b}=[0,\ldots,d_{1}]$, $\boldsymbol{c}=[d_{1},d_{1}-1,\ldots,1]$ with $m_{i}=x^{d_{1}-i}$, $p_{i}=x^{i}$ for $i=0,\ldots,d_{1}$.

In the following step, we have $p_{d_{1}+1}=\operatorname{proj}(x^{d_{1}+1},\boldsymbol{p},\boldsymbol{m})=x^{d_{1}+1}-\langle x^{d_{1}+1},m_{1}\rangle_{\sigma}p_{1}-\cdots-\langle x^{d_{1}+1},m_{d_{1}}\rangle_{\sigma}p_{d_{1}}=x^{d_{1+1}}$ such that $\langle x^{d_{1}+1},x^{j}\rangle_{\sigma}=0$ for $0\leqslant j\leqslant d$. The algorithm stops and outputs $\boldsymbol{b}=[1,\ldots,x^{d_{1}}]$, $\boldsymbol{c}=[x^{d_{1}},x^{d_{1}-1},\ldots,1]$, $\boldsymbol{k}=[x^{d_{1}+1}]$.

Example .2

We consider the function h(u_{1},u_{2})={\color[rgb]{0,1,0}{\color[rgb]{0,1,0}2+3}\hskip 2.5pt}\cdot{\color[rgb]{0,0,1}2^{u_{1}}2^{u_{2}}{\color[rgb]{0,1,0}{\color[rgb]{0,1,0}-}}3^{u_{1}}}. Its associated generating series is $\sigma=\sum_{\alpha\in\mathbbm{N}^{2}}h(\alpha)\boldsymbol{z}^{\alpha}=4+5z_{1}+7z_{2}+5z_{1}^{2}+11z_{1}z_{2}+13z_{2}^{2}+\cdots$.

At the first step, we have $\mathbf{x}^{\boldsymbol{b}}=[1]$, $\boldsymbol{p}=[1]$, $\boldsymbol{q}=\left[\frac{1}{4}\right]$. At the second step, we compute $\mathbf{x}^{\boldsymbol{b}}=[1,x_{1},x_{2}]$, $\boldsymbol{p}=[1,x_{1}-\frac{5}{4},x_{2}-\frac{9}{5}x_{1}-4]=[p_{1},p_{x_{1}},p_{x_{2}}]$ and $\boldsymbol{q}=\left[\frac{1}{4}p_{1},-\frac{4}{5}p_{x_{1}},\frac{5}{24}p_{x_{2}}\right]$. At the next step, we obtain $\boldsymbol{k}=[],\boldsymbol{d}=[x_{1}^{2},x_{1}x_{2},x_{2}^{2}]$.

[TABLE]

The matrix of multiplication by $x_{1}$ in the basis $\boldsymbol{p}$ is

[TABLE]

Its eigenvalues are {\color[rgb]{0,0,1}[1,2,3]} and the corresponding matrix of eigenvectors is

[TABLE]

that is, the polynomials $U(x)=[2-\frac{1}{2}\hskip 2.5ptx_{1}-\frac{1}{2}\hskip 2.5ptx_{2},-1+x_{2},\frac{1}{2}\hskip 2.5ptx_{1}-\frac{1}{2}\hskip 2.5ptx_{2}]$. By computing the Hankel matrix

[TABLE]

we deduce the weights {\color[rgb]{0,1,0}{\color[rgb]{0,1,0}2,3,-1}} and the frequencies {\color[rgb]{0,0,1}(1,1),} {\color[rgb]{0,0,1}(2,2),(3,1)}, which corresponds to the decomposition $\sigma=e^{y_{1}+y_{2}}+3e^{2y_{1}+2y_{2}}-e^{2y_{1}+y_{2}}$ associated to $h(u_{1},u_{2})=2+3\cdot 2^{u_{1}+u_{2}}-3^{u_{1}}$.

Example .3

We consider the following symmetric tensor or homogeneous polynomial:

[TABLE]

The associated series is

[TABLE]

To decompose it into a sum of powers of linear forms, we apply the border basis algorithm to the series $\sigma$. The algorithm projects successively the monomials $1,x_{1},x_{2},x_{1}^{2},x_{1}x_{2},x_{2}^{2},\ldots$ onto the family of polynomials $\boldsymbol{p}$, starting with $\boldsymbol{p}=[1]$. We obtain $\mathbf{x}^{\boldsymbol{b}}=\boldsymbol{c}=[1,x_{1},x_{2}]$, $\boldsymbol{p}=[1,x_{1}-6,x_{2}+\frac{1}{13}x_{1}-\frac{32}{13}]$ and the border basis is

[TABLE]

giving the projection of the border monomials \boldsymbol{d}=[{\color[rgb]{1,0,0}x_{1}^{2},x_{1}x_{2},x_{2}^{2}}] on the basis $\mathbf{x}^{\boldsymbol{b}}$. The decomposition of $\tau$ is deduced from the eigenvectors of the operator of multiplication by $x_{1}$:

[TABLE]

Its eigenvalues are [{\color[rgb]{1,0,0}-1,1,2]} and the eigenvectors correspond to the polynomials

[TABLE]

Computing $\omega_{i}=\langle\sigma\mid\boldsymbol{u}_{i}\rangle$ and $\xi_{i,j}=\frac{\langle\sigma\mid x_{j}\boldsymbol{u}_{i}\rangle}{\langle\sigma\mid\boldsymbol{u}_{i}\rangle}$ (see [22]), we obtain the decomposition:

[TABLE]

Example .4

We consider the algebraic code over $\mathbbm{K}=\mathbbm{Z}/32003\mathbbm{Z}$ defined by

[TABLE]

where

[TABLE]

and $\xi_{i}$ is the $i^{\mathrm{th}}$ column of $\Xi$. Suppose that we receive the word

[TABLE]

which is the sum $r=c+\omega$ of a code word $c\in C$ and an error vector $\omega\in\mathbbm{K}^{11}$. We want to correct it and find the corresponding word $c$ of the code $C$.

Computing the syndromes $\sigma_{\alpha}=\sum_{i=1}^{11}r_{i}\xi_{i}^{\alpha}=\sum_{i=1}^{11}\omega_{i}\xi_{i}^{\alpha}$ for $|\alpha|\leq 2$ and the corresponding (truncated) generating series, we get

[TABLE]

We apply the border basis algorithm to obtain error locator polynomials. The monomials are considered in the order $\mathbf{x}^{\mathbf{a}}=[1,x_{1},x_{2},x_{3},x_{1}^{2},x_{1}\,x_{2},\ldots,x_{3}^{2}]$. Here are the different steps, where $\mathfrak{n}$ denotes the new monomial introduced at each loop of the algorithm.

Step 1. $\mathfrak{n}=1$, $\mathbf{x}^{\boldsymbol{b}}=[1]$, $\mathbf{x}^{\boldsymbol{c}}=[x_{1}]$, $\boldsymbol{k}=[]$.

Step 2. $\mathfrak{n}=x_{1}$, $\mathbf{x}^{\boldsymbol{b}}=[1,x_{1}]$, $\mathbf{x}^{\boldsymbol{c}}=[x_{1},1]$, $\boldsymbol{k}=[]$.

Step 3. $\mathfrak{n}=x_{2}$, $\mathbf{x}^{\boldsymbol{b}}=[1,x_{1}]$, $\mathbf{x}^{\boldsymbol{c}}=[x_{1},1]$, $\boldsymbol{k}=[x_{2}+\frac{1}{2}x_{1}+\frac{3}{2}]$.

Step 4. $\mathfrak{n}=x_{3}$, $\mathbf{x}^{\boldsymbol{b}}=[1,x_{1}]$, $\mathbf{x}^{\boldsymbol{c}}=[x_{1},1]$, $\boldsymbol{k}=[x_{2}+\frac{1}{2}x_{1}+\frac{3}{2},x_{3}-1]$.

The algorithm stops at this step, since the new monomial $\mathfrak{n}=x_{1}^{2}$ is of degree $2$ and $\mathfrak{n}\cdot\mathbf{x}^{\boldsymbol{c}}\not\subset\mathbf{x}^{\mathbf{a}}$. It outputs two error locator polynomials: $x_{2}+\frac{1}{2}x_{1}+\frac{3}{2},x_{3}-1$.

We check that only $\xi_{5},\xi_{10}$ are roots of the error locator polynomials. We deduce the non-zero weights $\omega_{5},\omega_{10}$ by solving the system $\omega_{5}\xi_{5}^{\alpha}+\omega_{10}\xi_{10}^{\alpha}=\sigma_{\alpha}$ for $\alpha\in\{(0,0,0),(1,0,0)\}$. This yields $\omega_{5}=1,\omega_{10}=-1$, so that the code word is

[TABLE]