Mathieu Collowald

Evelyne Hubert
email: evelyne.hubert@inria.fr

Algorithms for Computing Cubatures Based on Moment Theory

Quadrature is an approximation of the definite integral of a function by a weighted sum of function values at specified points, or nodes, within the domain of integration. Gaussian quadratures are constructed to yield exact results for any polynomials of degree 2r -1 or less by a suitable choice of r nodes and weights. Cubature is a generalization of quadrature in higher dimension.

In this article we elaborate algorithms to compute all minimal cubatures for a given domain and a given degree. We propose first an algorithm in symbolic computation to characterize all cubatures of a given degree with a fixed number of nodes. The determination of the nodes and weights is then left to the computation of the eigenvectors of the matrix identified at the characterization stage and can be performed numerically.

The characterisation of cubatures on which our algorithms are based stems from moment theory. We formulate the results there in a basis independent way : Rather than considering the moment matrix, the central object in moment problems, we introduce the underlying linear map from the polynomial ring to its dual, the Hankel operator. This makes natural the use of bases of polynomials other than the monomial basis, and proves to be computationally relevant, either for numerical properties or to exploit symmetry.

Introduction

Quadrature and cubature generally refer to the approximation of integrals. Among the possible interpretations, this article looks at the approximation of the integral of a function f over a domain D in R n by a weighted sum of evaluations of f at well chosen points, known as nodes, in D. The nodes shall be independent of the function f and the weights are preferably positive. The criteria best suited to choose the nodes and weights can be open for debate and might depend on the intention or application. In this paper we adopt the degree criterion that presides over the well known Gaussian quadratures in the univariate case : the formula shall be exact for all polynomials up to a certain degree d. The algebraic techniques developed in this manuscript can nonetheless be adapted for the formula to be exact on alternative subspaces of polynomials, or even function spaces. Minimizing the number of nodes used for the criterion to be satisfied is the spirit in which the techniques are developed.

In dimension n = 1, Gaussian quadratures optimally achieve an exact formula for degree 2r -1 with r nodes. These nodes are the roots of the orthogonal polynomial of degree r associated with the considered integral operator. In higher dimension, the analogous problem does not enjoy such a clear-cut answer. Relevant integration schemes are still being investigated for standard domains like the simplex and the hypercube. Those are necessary for the discretization of partial differential equations, whether in high degree for isogeometric analysis [START_REF] Hughes | Efficient quadrature for nurbs-based isogeometric analysis[END_REF][START_REF] Mantzaflaris | Matrix generation in isogeometric analysis by low rank tensor approximation[END_REF] or in high dimensions [START_REF] Bungartz | Sparse grids[END_REF][START_REF] Le Maître | Spectral methods for uncertainty quantification[END_REF]. A simple product of Gaussian quadratures is "cursed by dimensionality". Our original motivation resides in the use of efficient integration schemes in geometric modeling [START_REF] Fuentes Suárez | Convolution surfaces with varying radius: Formulae for skeletons made of arcs of circles and line segments[END_REF][START_REF] Hubert | Convolution surfaces based on polygons for infinite and compact support kernels[END_REF][START_REF] Hubert | Convolution Surfaces based on Polygonal Curve Skeletons. (An Application of Symbolic Integration to Graphics)[END_REF].

Cubatures are the multidimensional analogues of quadratures. We mostly consider the linear forms

Ω ∶ R[x] → R p ↦ ∫ R n p(x)dµ(x)
where µ is a positive Borel measure on R n , the closed support supp µ of which is compact, and R[x] = R[x 1 , . . . , x n ] denotes the R-linear space of polynomials with n variables. A linear form

Λ ∶ R[x] → R p ↦ r j=1
a j p(ξ j )

with r ∈ N, a j > 0 and ξ j ∈ R n pairwise distinct, is called a cubature of degree d for Ω if it satisfies Λ(p) = Ω(p) for all p ∈ R[x] ≤d , where R[x] ≤d is the space of polynomials of degree d or less. The points ξ j are the nodes and the numbers a j are the weights of the cubature. Such a cubature is called inside if the nodes ξ j lie on supp µ and minimal if the number of nodes r is minimal.

When Ω is definite, i.e. when Ω(p 2 ) = 0 implies that p = 0, at least for p ∈ R[x] ≤d , it will appear clearly in our constructive approach that a cubature of degree d cannot have less than ⌊d 2⌋+n n nodes. This is known as Mysovskikh's theorem [START_REF] Dunkl | Orthogonal polynomials of several variables, volume 81 of Encyclopedia of Mathematics and its Applications[END_REF]Theorem 3.7.4]. A cubature that attains the lower bound is called a Gaussian cubature as it is then analogous to the univariate case. Gaussian cubatures are rather rare and the measures that admit Gaussian cubatures for any degree are intriguing [START_REF] Li | Discrete Fourier analysis, cubature, and interpolation on a hexagon and a triangle[END_REF][START_REF] Schmid | On bivariate Gaussian cubature formulae[END_REF]. A review of sharper lower bounds can be found in [START_REF] Möller | Lower bounds for the number of nodes in cubature formulae[END_REF], [START_REF] Cools | Cubature formulae and orthogonal polynomials[END_REF]Section 7].

The existence of inside cubatures of any degree is guaranteed by Tchakaloff's theorem [START_REF] Tchakaloff | Formules de cubatures mécaniques à coefficients non négatifs[END_REF] [70, Section 3.3].

It also provides a bound on the number of nodes of minimal cubatures, namely d+n n , where d is the degree of the cubature. Analogous results [START_REF] Bayer | The proof of Tchakaloff's theorem[END_REF][START_REF] Putinar | A note on Tchakaloff's theorem[END_REF] were later proved with less constraint on the measure µ and recently this upper bound was revisited for cubatures in the plane [START_REF] Riener | Optimization approaches to quadrature: new characterizations of Gaussian quadrature on the line and quadrature with few nodes on plane algebraic curves, on the plane and in higher dimensions[END_REF]. However, the proof of Tchakaloff's theorem is not constructive.

Radon's work [START_REF] Radon | Zur mechanischen Kubatur[END_REF] marked a starting point in the construction of cubatures thanks to orthogonal polynomials. As presented in reviews on the topic [START_REF] Cools | Constructing cubature formulae: the science behind the art[END_REF][START_REF] Stroud | Approximate calculation of multiple integrals[END_REF], cubatures have been constructed using several techniques based on zeros of multivariate orthogonal polynomials, on product of Gaussian quadratures, or on solutions of multivariate nonlinear systems. As a result, cubature formulae are available for a large variety of regions and weight functions [START_REF] Cools | Monomial cubature rules since "Stroud": a compilation[END_REF][START_REF] Stroud | Approximate calculation of multiple integrals[END_REF].

More recently, a characterization of cubatures based on moment matrices was initiated in [START_REF] Fialkow | A moment matrix approach to multivariable cubature[END_REF]. Additional contributions in this direction were made regarding Gaussian cubatures [START_REF] Lasserre | The existence of Gaussian cubature formulas[END_REF] or on a numerical method to compute some of them [START_REF] Abril Bucero | On the construction of general cubature formula by flat extensions[END_REF]. They are based on algebraic properties of moment matrices that emerged in the realm of global optimization [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF][START_REF] Laurent | A generalized flat extension theorem for moment matrices[END_REF].

The moment matrix relating to the linear form Λ can be interpreted as the matrix of a linear map from the polynomial ring to its dual in the monomial basis. In the univariate case, a moment matrix is Hankel. Alternative bases are relevant, for instance when symmetry is to be accounted for [START_REF] Collowald | Cubatures with symmetry[END_REF] or when better numerical stability is sought [START_REF] Beckermann | How to choose modified moments?[END_REF], [START_REF] Golub | Matrices, moments and quadrature with applications[END_REF]Section 5.3]. Rather than focus on the moment matrix, we introduce the Hankel operator as the central object. An algebraic characterisation of cubatures can then be expressed in terms of the rank and the positivity of the Hankel operator. This characterisation translates naturally into linear algebra for any choice of a basis for the polynomial ring.

We present an algorithm in symbolic computation that provides necessary and sufficient conditions for the existence of a cubature of degree d with r nodes. The algorithm improves on Gauss-Bareiss elimination by taking advantage of the positivity constraint. Beside being fraction free, the diagonalisation it performs avoids the complicated pivoting thanks to a pertinent observation about positive semidefinite matrices. The output of the algorithm consists of sets equalities and inequalities bearing on unknown generalized moments, up to a determined finite order. Each solution of these correspond in a unique way to a cubature. The nodes and the weights of the characterized cubature can then be computed as eigenvectors of a matrix identified at the first stage. This second step can safely be performed numerically.

The paper is organized as follows. Section 3 introduces the Hankel operator associated to a linear form and shows how to characterize cubatures with linear algebra. Section 4 describes the fraction and pivoting free algorithm to diagonalize parametric symmetric matrices over their locus of positivity. Section 5 makes precise how the tools exposed so far apply to algorithmically characterize the cubatures for a given integration form and compute its nodes and weights. Section 6 applies the algorithm to find all the minimal cubatures of degree 5 for the hexagon. We show that there is a one parameter family of such cubatures, when only two were known. We also demonstrate on this example how symmetry can be exploited to alleviate the computational cost.

In Section 2 we treat anew the well known univariate case of Gaussian quadratures with a view to its present generalization to cubatures. We then exhibit a number of related univariate problems that can be generalized with the multivariate techniques presented in this mansucript.

The univariate case : Gaussian quadratures

In this section, the classical quadrature problem is introduced (see e.g. [START_REF] Stoer | Introduction to numerical analysis[END_REF]Chapter 3] or [START_REF] Stroud | Approximate calculation of multiple integrals[END_REF]Chapter 1.3] and references herein) and methods to determine Gaussian quadratures are presented. Though the subject is tied in with orthogonal polynomials (see e.g. [21, Chapter 1]), we favor here a presentation that extends to the multidimensional case that we shall deal with. Here, R[x] denotes the R-vector space of polynomials in one variable and coefficients in R. Then R[x] ≤d is the space of polynomials of degree d and less. A linear form Λ on R[x] defined by

Λ ∶ R[x] → R, p ↦ r j=1
a j p(ξ j ), with r ∈ N, r ≥ 1, a j ∈ R, a j > 0, and ξ j ∈ R, pairwise distinct, is called a quadrature of degree d for Ω if it satisfies

Ω(p) = Λ(p), i.e. b a p(x)ω(x)dx = r j=1 a j p(ξ j ), ∀p ∈ R[x] ≤d
and if this equality does not hold for at least one polynomial p of degree d + 1. The points ξ j are the nodes and the coefficients a j are the weights of the quadrature. Such a quadrature is called inside if the nodes ξ j lie on the interval [a, b] and minimal if the number of nodes r is minimal for a fixed degree d.

In this univariate case, minimal inside quadratures with positive weights are known and called Gaussian quadratures (see e.g. [START_REF] Stoer | Introduction to numerical analysis[END_REF]Chapter 3.6] for more details). A Gaussian quadrature with r nodes is of degree 2r -1.

In next paragraph we show two related ways of computing the nodes. This is done thanks to the properties of the moment matrices:

Ω(x i+j-2 ) 1≤i,j≤r+1 = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ Ω(1) Ω(x) . . . Ω(x r ) Ω(x) . . . ⋮ ⋮ . . . ⋮ Ω(x r ) . . . . . . Ω(x 2r ) ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ m 0 m 1 . . . m r m 1 . . . ⋮ ⋮ . . . ⋮ m r . . . . . . m 2r ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = m i+j-2 1≤i,j≤r+1 .
These are Hankel matrices. But rather than their displacement structure, let us observe that these are the matrices of the scalar product Φ

∶ R[x] × R[x] → R, (p, q) ↦ Ω(pq)
restricted to R-vector space R[x] ≤r endowed with the monomial basis {1, x, . . . , x r }. The bilinear map defines a duality pairing to which is associated the injective map

Ĥ ∶ R[x] → R[x] * p ↦ Λ p , where Λ p ∶ R[x] → R q ↦ Λ(p q).
In the multivariate case we shall actually focus on this Hankel operator Ĥ rather than the bilinear map Φ.

Once the nodes of the cubature are found, the weights of a quadrature are generally obtained afterwards by solving the Vandermonde linear system

⎛ ⎜ ⎜ ⎜ ⎝ 1 1 . . . 1 ξ 1 ξ 2 . . . ξ r ⋮ ⋮ ⋮ ξ r-1 1 ξ r-1 2 . . . ξ r-1 r ⎞ ⎟ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎜ ⎝ a 1 a 2 ⋮ a r-1 ⎞ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎝ m 0 m 1 ⋮ m r-1 ⎞ ⎟ ⎟ ⎟ ⎠ (1) obtained from the equations Λ(x k ) = r j=1 a j ξ k j = m k = Ω(x k ) ∀k = 0, . . . , r -1. (2)

Nodes as roots of orthogonal polynomials

Assume that a Gaussian quadrature of degree 2r -1 exists for Ω, i.e. there exists a linear form

Λ ∶ R[x] → R, p ↦ r j=1 a j p(ξ j )
with a j > 0 and

ξ j ∈ R pairwise distinct such that Ω(x k ) = m k = Λ(x k ) for all k ≤ 2r -1.
Consider then the monic polynomial π r of degree r whose roots are the nodes ξ 1 , . . . , ξ r

π r (x) = r j=1 (x -ξ j ) = x r -τ r-1 x r-1 -. . . -τ 0 .
Since π r (ξ j ) = 0 for all j = 1, . . . , r, we have Λ(π r ) = 0 and even Λ(pπ r ) = 0 ∀p ∈ R[x]. Taking p(x) = x k for k = 0, . . . , r, we deduce by linearity of Λ that r-1 i=0 τ i Λ(x k+i ) -Λ(x k+r ) = 0 ∀k = 0, . . . , r.

Thus, the vector τ 0 . . . τ r-1 -1 t is in the kernel of the Hankel matrix

H (r) 1 = Λ(x i+j-2 ) 1≤i,j≤r+1 = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ Λ(1) Λ(x) . . . Λ(x r ) Λ(x) . . . ⋮ ⋮ . . . ⋮ Λ(x r ) . . . . . . Λ(x 2r ) ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ . By hypothesis Λ(x k ) = m k = Ω(x k ) for 0 ≤ k ≤ 2r -1.
To obtain the nodes of Λ one can thus:

• Determine the polynomial π r , through its coefficients τ 0 , . . . , τ r-1 in the monomial basis, by solving the Hankel linear system

⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ m 0 m 1 . . . m r-1 m 1 . . . ⋮ ⋮ . . . ⋮ m r-1 . . . . . . m 2r-2 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ τ 0 τ 1 ⋮ ⋮ τ r-1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ m r m r+1 ⋮ ⋮ m 2r-1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
• Determine the roots of π r .

To make the link with a classical approach to Gaussian quadratures, we remark here that the polynomial π r constructed above is the monic orthogonal polynomial of degree r for the scalar product Φ defined by Ω. Indeed, since Λ is a cubature of degree 2r -1, we have

Φ(p, π r ) = Ω(pπ r ) = Λ(pπ r ) = 0, ∀p ∈ R[x] ≤r-1 . (3)
Orthogonal polynomials have real roots, within the support [a, b] of the measure [START_REF] Szegö | Orthogonal polynomials[END_REF], as is expected for a quadrature.

Nodes as generalized eigenvalues

Adding some columns to the Hankel linear system above, which results in shifting left the columns of

H (r-1) 1 , one observes ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ m 0 m 1 . . . m r-1 m 1 . . . ⋮ ⋮ . . . ⋮ m r-1 . . . . . . m 2r-2 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ H (r-1) 1 ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0 . . . . . . 0 τ 0 1 ⋱ ⋮ τ 1 0 ⋱ ⋱ ⋮ ⋮ ⋮ ⋱ ⋱ 0 ⋮ 0 . . . 0 1 τ r-1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ M (r-1) x = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ m 1 m 2 . . . m r m 2 . . . ⋮ ⋮ . . . ⋮ m r . . . . . . m 2r-1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ H (r-1) x , that is H (r-1) 1 M (r-1) x = H (r-1)
x , where H (r-1)

q = Λ(q x i+j-2 ) 1≤i,j≤r for q ∈ R[x]. (4)
M (r-1) x
is the companion matrix of the polynomial π r . Its eigenvalues are thus the roots of the polynomial π r . Hence the generalized eigenvalues of the pair of Hankel matrices (H

(r-1) x , H (r-1) 1
) are the sought nodes ξ 1 , . . . , ξ r .

Both H (r-1) 1 and H (r-1) x are real and symmetric matrices. As the matrix of a scalar product, H (r-1) 1 is furthermore positive definite. Let v be a generalized eigenvector of (H

(r-1) x , H (r-1) 1 ) associated to the eigenvalue v, i.e. H (r-1) x v = ξH (r-1) 1 v. Then ξ = vT Hxv vT H1v = v T Hx v v T H1 v = ξ.
This proves that the generalized eigenvalues of (H (r-1) x , H (r-1) 1

) are real.

In the multivariate setting, the matrices of the Hankel operator will take the role of the above Hankel matrices. The way to generalize the companion matrix M (r-1) x is to see it as the matrix of the multiplication by x, modulo π r , expressed in the monomial basis {1, x, . . . , x r } of quotient space R[x] (π r ).

Jacobi matrix and modified moment matrices

In the above, the matrices that come in the determination of a quadrature were interpreted as matrices of symmetric bilinear forms or linear maps in the monomial basis. Alternative bases are of interest from different points of view.

On one hand we considered the bilinear forms Consider first the basis Π defined by the orthonormal polynomials πr , obtained by normalizing the monic orthogonal polynomials π r . Then

Φ ∶ R[x] × R[x] → R, (p, q) ↦ b a p(x)q(x)ω(x)dx and Φ x ∶ R[x] × R[x] → R, (p, q) ↦ b a xp(x)q(x)ω(x)dx
• The matrix of Φ in the basis Π is the identity.

• The matrix of Φ x in the basis Π is the tridiagonal Jacobi matrix

⎛ ⎜ ⎜ ⎜ ⎝ a 0 b 0 0 . . . b 0 a 1 b 1 ⋱ 0 b 1 ⋱ ⋱ ⋮ ⋱ ⋱ ⋱ ⎞ ⎟ ⎟ ⎟ ⎠
, where (a r ) r∈N and (b r ) r∈N are sequences of real numbers that reflect the recurrence equation of order

xπ r (x) = b r πr+1 (x) + a r πr (x) + b r-1 πr-1 (x) ∀r ≥ 1 xπ 0 (x) = b 0 π1 (x) + a 0 π0 (x)
The matrix of the multiplication by x in R[x] (π r ) is thus the truncated r × r Jacobi matrix. Its eigenvalues are the roots of the orthonormal polynomial πr , or equivalently the nodes ξ 1 , . . . , ξ r . See [START_REF] Stoer | Introduction to numerical analysis[END_REF]Theorem 3.6.20] for the classical link between eigenvalues of tridiagonal matrices and roots of orthogonal polynomials.

One may also consider the matrices of Φ and Φ x in a basis of orthogonal polynomials for another measure. This correspond to considering modified moments and provide some computational advantage in terms of numerical stability [START_REF] Beckermann | How to choose modified moments?[END_REF], [START_REF] Sack | An algorithm for Gaussian quadrature given modified moments[END_REF].

Another interesting choice of alternative basis arises in the presence of symmetry. Assume that the weight function ω ∶ [-a, a] → R is even. Then the odd moments are zero : Ω(x 2k+1 ) = 0. If we simply reorder the monomials in a basis by listing first the even degree monomials and then the odd degree monomials, we observe that the matrix H1 of Φ in this basis becomes block diagonal. The same is true of the matrix Hx 2 of Φ x 2 and the generalized eigenvalues of the pair of truncated matrices ( Hx 2 , H1 ) are ξ 1 1 , . . . , ξ 2 r . This all generalizes to the multivariate case that will be developed in [START_REF] Collowald | Cubatures with symmetry[END_REF].

Related problems

The techniques presented in this paper have other applications than cubatures. Indeed, a number of classical problems can be formulated as we did for quadratures, in the univariate case, and cubatures, in the multivariate case.

In the univariate case, one seeks to retrieve pairs {(a 1 , ξ 1 ), . . . , (a r , ξ r )} from the moments (µ k ) k∈N or (µ k ) 0≤k≤R with R an upper bound of r under the assumption that

µ k = r j=1 a j ξ k j .
In the case of quadratures, µ k ∈ R and we expect a j and ξ j to be in R. In general though µ k , a j and ξ j are in C. Furthermore, the number of terms r might be an additional unknown of the problem. As demonstrated above the problem can be solved (uniquely) when R = 2r -1 moments are available.

In the multivariate case, the input is indexed by k = (k 1 , . . . , k n ) ∈ N n and each ξ j belongs to C n . The relationships are given by

µ k = r j=1
a j ξ k1 j,1 . . . ξ kn j,n .

An additional problem in the multivariate case, even when r is known, is to find an appropriate K-vector space basis of the quotient algebra K[x 1 , . . . , x n ] I ξ , where I ξ is the ideal of all polynomials vanishing on {ξ 1 , . . . , ξ r } ⊂ C n .

Prony's method for exponential interpolation [START_REF] Baron De Prony | Essai expérimental et analytique sur les lois de la dilatabilité des fluides élastique et sur celles de la force expansive de la vapeur de l'eau et de la vapeur de l'alkool, à différentes températures[END_REF]. One wants to retrieve f ∶ C → C as

f (x) = r j=1
a j e λj x from the knowledge of µ k = f (k) with k ∈ N. One retrieves the above problem by considering ξ j = e λj .

In its multivariate version f ∶ C n → C and

f (x) = r j=1
a j e ⟨λj , x⟩ where x = (x 1 , . . . , x n ), λ j ∈ C n and ⟨y, x⟩ = y 1 x 1 + . . . + y n x n . With ξ j1 = e λj1 , . . . , ξ jn = e λjn the problem is then recast into

f (x) = r j=1 a j ξ x j .
Under some natural assumptions, solutions have been proposed for the univariate case [START_REF] Potts | Parameter estimation for exponential sums by approximate prony method[END_REF], for the multivariate case with a univariate resolution (projection method) [START_REF] Potts | Parameter estimation for multivariate exponential sums[END_REF] and with a multivariate approach [START_REF] Kunis | A multivariate generalization of prony's method[END_REF][START_REF] Sauer | Prony's method in several variables: symbolic solutions by universal interpolation[END_REF][START_REF] Mourrain | Polynomial exponential decomposition from moments[END_REF].

Sparse interpolation [START_REF] Ben-Or | A deterministic algorithm for sparse multivariate polynomial interpolation[END_REF]. Assuming a (multivariate) polynomial f has a small support

f (x) = r j=1 a j x αj ,
one wishes to retrieve the exponents α j ∈ N n from evaluations of the polynomial at chosen points. For k ∈ N, one chooses µ k = f (p k 1 , . . . , p k n ) where p j are distinct prime numbers. From µ 0 , . . . , µ 2r-1 one can retrieve ξ j = p αj = p αj1 1 . . . p αjn n so that the exponent can be found by factorization [START_REF] Ben-Or | A deterministic algorithm for sparse multivariate polynomial interpolation[END_REF][START_REF] Kaltofen | Improved sparse multivariate polynomial interpolation algorithms[END_REF][START_REF] Kaltofen | Early termination in sparse interpolation algorithms[END_REF]. Replacing the prime numbers by roots of unity, the authors in [START_REF] Giesbrecht | Symbolic-numeric sparse interpolation of multivariate polynomials[END_REF] proposed a symbolic-numeric solution.

For latest development in Prony's method and sparse interpolation, one can consult the material attached to the 2015 Dagstuhl seminar Sparse modeling and multi-exponential analysis (http://www.dagstuhl.de/ 15251).

Pole estimation, Padé approximant [START_REF] Henrici | Applied and computational complex analysis[END_REF]. The input data are the coefficients of the Taylor expansion at z = 0 of a function f ∶ C → C and one wishes to find the poles of the function.

f (x) = r j=1 a j 1 -ξ j x = k∈N µ k x k .
The multivariate version that can be approached is

f (x) = r j=1 a j ⟨ξ j , x⟩ = k∈N n µ k x k .
In the univariate case all polynomials can be factored over C into degree 1 polynomials. Therefore the problem covers all rational functions with distinct poles. The restriction to denominators that can be factored into linear form is more restrictive in the multivariate case.

Shape-from-moment problem [START_REF] Gravin | The Inverse Moment Problem for Convex Polytopes[END_REF][START_REF] Milanfar | Reconstructing polygons from moments with connections to array processing[END_REF]. The problem consists in recovering the vertices of a convex n-dimensional polytope V from its moments. Projection methods have been developed in [START_REF] Collowald | Numerical reconstruction of convex polytopes from directional moments[END_REF][START_REF] Gravin | The Inverse Moment Problem for Convex Polytopes[END_REF] based on moments in a direction δ ∈ R n and Brion's identities

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ (k+n)! k! V < x, δ > k dx = r j=1 a j < v j , δ > n+k k ≥ 0 0 = r j=1 a j < v j , δ > n-k 1 ≤ k ≤ n.
The coefficients a j are nonzero real numbers and < v j , δ > are the projections of the vertices v j on the direction δ. Taking ξ j =< v j , δ >, the formula can be recast into

µ k = r j=1 a j ξ k j k ∈ N,
where µ k is related to the left hand side of the system of equations above. The set of projected vertices can thus be recovered. Different projections are then required and matching processes are presented in [START_REF] Collowald | Numerical reconstruction of convex polytopes from directional moments[END_REF][START_REF] Gravin | The Inverse Moment Problem for Convex Polytopes[END_REF].

The case n = 2 treated with complex moments in [START_REF] Elad | Shape from moments-an estimation theory perspective[END_REF][START_REF] Golub | A stable numerical method for inverting shape from moments[END_REF][START_REF] Milanfar | Reconstructing polygons from moments with connections to array processing[END_REF] is linked to this general case in [START_REF] Collowald | Numerical reconstruction of convex polytopes from directional moments[END_REF].

A multidimensional treatment of the shape-from-moment problem is however impossible due to the lack, up to our knowledge, of an equivalent formula for the moments

∫ V x k dx with k ∈ N n .
Symmetric tensor decomposition [START_REF] Brachat | Symmetric tensor decomposition[END_REF]. The general (multivariate) problem is to find, for a homogeneous polynomial

f (z) = k =d d k µ k z k ,
the minimal rank r such that there exist (a 1 , ξ 1 ), . . . , (a r , ξ r ) ∈ C × C n such that f can be written as a linear combination of linear forms to the d th power:

f (z) = r i=1 a i ⟨ξ i , z⟩ d .
De-homogenizing the binary case (n = 2), we obtain an equivalent univariate problem

f (z) = µ d z d + . . . + d k µ k z k + . . . + µ 0 = r i=1 a i (1 -ξ i z) d
initially solved by Sylvester [START_REF] Sylvester | Sur une extension du théorème de clebsch relatif aux courbes du quatrième degré[END_REF].

Characterization of cubatures through Hankel operators

In the following, K denotes a subfield of C. At some points we shall need to make clear whether K is C, R or a computable extension of Q, in the sense that exact arithmetic can be performed in

K. K[x] = K[x 1 , . . . , x n ]
denotes the ring of polynomials in the variables x = (x 1 , . . . , x n ) with coefficients in K and K[x] ≤δ the K-vector space of polynomials of degree at most δ.

The dual of K[x], the set of K-linear forms on K[x], is denoted by K[x] * . A typical example of a linear form on K[x] is the evaluation e ξ at a point ξ of K n . It is defined by

e ξ ∶ K[x] → K, p ↦ p(ξ).
Other examples of linear forms on K[x] are given by linear combinations of evaluations

Λ ∶ K[x] → K, p ↦ r j=1
a j e ξj (p), with a j ∈ K ∖ {0} and ξ j ∈ K n , or by the integration over a domain

D ⊂ R n Ω ∶ R[x] → R, p ↦ D p(x)dx.
We shall say that Λ

∶ R[x] → R is a cubature if Λ = ∑ r j=1
a j e ξj , with the nodes ξ j ∈ R n , and the weights a j > 0.

This section is first aimed at characterizing such linear forms Λ through their associated Hankel operators : cubatures correspond to positive semi-definite Hankel operators of finite rank. Second, we present a criterion of existence for such a form Λ to exist from the knowledge of the values of Λ on K[x] ≤δ , for some δ < ∞ : Hankel operators of a given finite rank can be uniquely determined from their restriction to K[x] ≤δ thanks to the flat extension criterion. Last, we shall show how to compute the nodes ξ 1 , . . . , ξ r and the weights a 1 , . . . , a r from the values of Λ on K[x] ≤δ : The nodes of the cubature are the roots of the ideal formed by the kernel of the Hankel operator. They can be computed as generalized eigenvalues and eigenvectors.

When reviewing the above results as needed for cubatures, we strive to remain simple, providing the arguments and proofs that basically rely on linear algebra and matrix computations. A noteworthy feature of our review is that we free the results from the monomial basis by making the Hankel operator the central object rather than the moment matrix. Alternative bases are indeed useful, for instance when symmetry is to be accounted for [START_REF] Collowald | Cubatures with symmetry[END_REF] or when better numerical stability is sought [START_REF] Beckermann | How to choose modified moments?[END_REF], [START_REF] Golub | Matrices, moments and quadrature with applications[END_REF]Section 5.3]. Our formulations are thus readily amenable for application in a more varied context.

The concept of flat extension was introduced for the classical truncated moment problem [START_REF] Curto | Solution of the truncated complex moment problem for flat data[END_REF][START_REF] Curto | The truncated complex K-moment problem[END_REF]. It was incorporated in a successful line of developments in global optimisation and algebraic geometry [START_REF] Lasserre | Global optimization with polynomials and the problem of moments[END_REF][START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF]. Our reference text is the book chapter [START_REF] Laurent | Sums of squares, moment matrices and optimization over polynomials[END_REF]. The link between cubatures and flat extensions was first expanded on in [START_REF] Fialkow | A moment matrix approach to multivariable cubature[END_REF]; [START_REF] Lasserre | The existence of Gaussian cubature formulas[END_REF] then provided a criterion for the existence of a Gaussian cubature while [START_REF] Abril Bucero | On the construction of general cubature formula by flat extensions[END_REF] showed how to compute some cubatures as a global optimization problem. The connection between Hankel operators and multiplication maps, that allows one to compute the cubature nodes as an eigenvalue problem, was for instance used in [START_REF] Brachat | Symmetric tensor decomposition[END_REF][START_REF] Mourrain | Polynomial exponential decomposition from moments[END_REF]. in the context of symmetric tensor decomposition and multivariate Prony's method.

Hankel operator associated to a linear form

For a linear form Λ ∈ K[x] * , the associated Hankel operator Ĥ is the K-linear map

Ĥ ∶ K[x] → K[x] * p ↦ Λ p , where Λ p ∶ K[x] → K q ↦ Λ(p q).
The kernel

I Λ ∶= {p ∈ K[x] Λ p = 0} = {p ∈ K[x] Λ(p q) = 0, ∀q ∈ K[x]} of Ĥ is an ideal of K[x]. Let us observe that if Λ = ∑ r j=1
a j e ξj , with a j ∈ K and ξ j ∈ K n , then the ideal

I(ξ) ∶= {p ∈ K[x] p(ξ i ) = 0, 1 ≤ i ≤ r} is included in I Λ .
We shall see that the reciprocal is true as well.

The Hankel operator can be represented by an infinite matrix once we fix a basis for K[x] and K[x] * , but we shall examine the truncated operators to remain in the realm of finite dimension when it comes to matrices. If V is a K-linear subspace of K[x] we define Ĥ V thanks to the restrictions Λ p V of Λ p to V :

Ĥ V ∶ V → V * p ↦ Λ p V , Assume V is the K-linear span ⟨b 1 , . . . , b s ⟩ of a linearly independent set B = {b 1 , . . . , b s } in K[x].
Then the matrix of Ĥ V in the basis B and its dual basis is

H B 1 = (Λ(b i b j )) 1≤i,j≤s .
Note that if Λ = ∑ r i=1 a i e ξj we have the following matrix equality that generalizes the factorization of Hankel matrices in terms of Vandermonde matrices:

H B 1 = W B ξ T A W B ξ , where A = diag(a 1 , . . . , a r ) and W B ξ = (b j (ξ i )) 1≤i≤r,1≤j≤s .
Hence Ĥ has rank at most r when Λ = ∑ r j=1 a j e ξj . In all cases, the rank of Ĥ is the dimension of K[x] I Λ as a K-vector space.

Before further investigation of Hankel operators, we shall review some results about ideals such as I(ξ) and I Λ , when Ĥ is of finite rank.

Zero-dimensional ideals

In the univariate case, the nodes of a quadrature are characterized as the roots of a univariate polynomial p.

If p ∈ K[x] is of degree r then K[x] = ⟨1, x, . . . , x r-1 ⟩ ⊕ (p),
where (p) is the ideal generated by p. In other words the classes of 1, x, . . . , x r-1 in the quotient K[x] (p) form a basis of this quotient when considered as a K-linear space. The number of roots of p, counted with multiplicities, is equal to the dimension of K[x] (p). Observe furthermore that the companion matrix associated to p corresponds to the matrix of the multiplication by x in K[x] (p); Its eigenvalues are the roots of p. In this section, we generalize these notions to the multivariate case. The set of nodes of a cubature shall be characterized as the roots, or variety, of a zero dimensional ideal in the multivariate polynomial ring

K[x] = K[x 1 , . . . , x n ].
Given an ideal I in K[x], its variety is understood in C n and is denoted by

V C (I) = {ξ ∈ C n p(ξ) = 0, ∀p ∈ I} . Conversely, for V ⊂ C n , I(V ) = {p ∈ K[x] p(ξ) = 0, ∀ξ ∈ V } is a radical ideal.
In the univariate case, radical ideals are generated by square-free polynomials. In the multivariate case, an ideal

I is radical if p k ∈ I ⇒ p ∈ I, for any k > 1. Theorem 3.1 [17, Chapter 2.2] Let I ⊂ K[x] be an ideal. The variety V C (I) ⊂ C n is a finite set if and only if the algebra K[x]
I is finite-dimensional as a K-vector space. Then the dimension of K[x] I as a K-vector space is greater than or equal to the number of points in V C (I). Equality occurs if and only if I is a radical ideal.

When I is not a radical ideal, an appropriate notion of multiplicity for the roots of I allows to obtain a more general version of the above result [START_REF] Cox | Using algebraic geometry[END_REF]Chapter 4]. In either case, an ideal I is said to be zero-dimensional if K[x] I is finite-dimensional as a K-vector space. Contrary to the univariate case, there is no universal bases for the quotient spaces K[x] I of a given dimension.

The class of any polynomial

f ∈ K[x] modulo I is denoted by [f ].
To shorten some statements we shall say that

B = {b 1 , . . . , b r } ⊂ K[x] is a basis of I Λ when [B] = {[b 1 ], . . . , [b r ]} is. Given a polynomial p ∈ K[x]
, the map proof:

M p ∶ K[x] I → K[x] I [f ] ↦ [pf ]
If M B p = m ij 1≤i
,j≤r is the matrix of the multiplication operator M p , we have

pb j ≡ ∑ r i=1 m ij b i mod I, ∀j = 1, . . . , r. Evaluating this equality at a root ξ of I brings p(ξ)b j (ξ) = r k=1 m kj b k (ξ) ∀j = 1, . . . , r, which can be rewritten as p(ξ i ) b 1 (ξ i ) . . . b r (ξ i ) = b 1 (ξ i ) . . . b r (ξ i ) M B p . ∎ Assume V C (I) = {ξ 1 , . . . , ξ s }. Introducing the (generalized) Vandermonde matrix W B ξ = (b j (ξ i )) 1≤i≤s, 1≤j≤r (5)
the above result can be written as a matricial identity

W B ξ M B p = DW B ξ ,
where D is the diagonal matrix whose diagonal elements are p(ξ 1 ), . . . , p(ξ s ).

As in the univariate case, the Vandermonde matrix appears naturally in Lagrange interpolation. We shall call a set { 1 , . . . , r } ⊂ K[x] a Lagrange basis associated to ξ = {ξ 1 , . . . , ξ r } if it satisfies i (ξ j ) = δ ij . The above theorem shows their existence. It is unique only modulo the ideal I(ξ) of all the polynomials vanishing at ξ 1 , . . . , ξ r .

Theorem 3.3 Let I be a zero-dimensional radical ideal in K[x] with V C (I) = {ξ 1 , . . . , ξ r } ⊂ K n . Then B = {b 1 , . . . , b r } is a basis of K[x] I if and only if the matrix W B ξ = (b j (ξ i )) 1≤i,j≤r is invertible. In this case the r polynomials 1 , . . . , r in K[x] defined by 1 . . . r = b 1 . . . b r W B ξ -1 (6) form a basis of K[x] I that satisfies i (ξ j ) = δ ij . For any p ∈ K[x], [ i ] is an eigenvector of M p associated to the eigenvalue p(ξ i ). proof: Let α = (α 1 , . . . , α r ) ∈ K r such that α 1 b 1 + ⋅ ⋅ ⋅ + α r b r ≡ 0 mod I. Since V C (I) = {ξ 1 , . . . , ξ r }, we have α 1 b 1 (ξ i ) + ⋅ ⋅ ⋅ + α r b r (ξ i ) = 0 ∀i = 1, . . . , r.
With the help of the matrix W B ξ = (b j (ξ i )) 1≤i,j≤r , this equality can be rewritten as W B ξ α = 0. Thus, the matrix W B ξ is invertible if and only if (α 1 , . . . , α r ) = (0, . . . , 0). This means that

W B ξ is invertible if and only if b 1 , . . . , b r are linearly independent modulo I. Since Card B = dim K[x] I, we have that W B ξ is invertible if and only if B is a basis of K[x] I.
Evaluating the defining equality for at ξ i we obtain that v = 1 (ξ i ) . . . r (ξ r ) is the solution to the linear system vW B ξ = b 1 (ξ i ) . . . b r (ξ i ) . The unique solution to this system is the row vector that has zero entries except for a 1 at the i-th position.

The last point comes from the matricial equality W B ξ M B p = DW B ξ following Theorem 3.2, mutipliying it on the right and left by W B ξ -1 . ∎

Hankel operators of finite rank

In this section we examine the Hankel operators associated to linear forms ∑ r i=1 a i e ξj . Several results in this section require additionally that I Λ is radical. This property will appear in Section 3.4 as a consequence of the positivity of the Hankel operator associated to cubature. The goal here is to show how the weights a 1 , . . . , a r and the nodes ξ 1 , . . . , ξ r can be retrieved as generalized eigenvectors.

Theorem 3.4 A linear form Λ is equal to ∑ r i=1 a i e ξj for some ξ i ∈ K n and a i ∈ K ∖ {0}, if and only if its associated Hankel operator Ĥ is of finite rank r and I Λ radical. In this case, V C (I Λ ) = {ξ 1 , . . . , ξ r } and a i = Λ( i ), where { 1 , . . . , r } is a Lagrange basis associated to {ξ 1 , . . . , ξ r }.

proof: The rank of Ĥ is the dimension of K[x] I Λ . Assume that K[x] I Λ is of dimension r and that I Λ is radical. There are thus distinct ξ 1 , . . . , ξ r ∈ C n that form the variety of I Λ . We can call on Theorem 3.3 to exhibit a Lagrange basis for K Let f ∈ I(ξ), then, for all q ∈ K[x], Λ(qf ) = ∑ r j=1 a j q(ξ j )f (ξ j ) = 0. Thus, I(ξ) ⊂ I Λ . Conversely, consider q ∈ I Λ and take 1 , . . . , r to be a Lagrange basis associated to ξ 1 , . . . , ξ r . Then a i q(ξ i ) = Λ(q i ) = 0 for all 1 ≤ i ≤ r. Since a i ≠ 0, we conclude that q ∈ I(ξ) so that I Λ ⊂ I(ξ). ∎ One can formulate a similar result without assuming that I Λ is radical [55, Theorem 3.1], with some derivation operators coming into play in Λ. One then touches on Hermite interpolation instead of Lagrange interpolation.

[x] I Λ . Since K[x] = ⟨ 1 , . . . , r ⟩ ⊕ I Λ , Λ
In the spirit of Section 3.2, our next endeavor is to characterize the bases K[x] I Λ and determine the matrices of the multiplication maps therein. This shall be done in terms of the matrices of the Hankel operators and their inverses.

Bases of K[x] I Λ

A basis of K[x] I Λ as a K-vector space is the image by the natural projection K

[x] → K[x] I Λ of a linearly independent set B ⊂ K[x] s.t. K[x] = ⟨B⟩ ⊕ I Λ .
Recall from linear algebra, see for instance [30, Proposition V, Section 2.30], that this latter equality implies:

K[x] * = (I Λ ) ⊥ ⊕ B ⊥ and (I Λ ) ⊥ → ⟨B⟩ * Λ ↦ Λ ⟨B⟩ is an isomorphism,
where, for any

V ⊂ K[x], V ⊥ = {Γ ∈ K[x] * Γ(v) = 0, ∀v ∈ V } is a linear subspace of K[x] * .
Note that the image of Ĥ lies in (I Λ ) ⊥ . With the natural identification of ⟨B⟩ with K[x] I Λ , and ⟨B⟩ * with

(K[x] I Λ )
* , the factorization of Ĥ by the natural projection

K[x] → K[x] I Λ defines the injective morphism H ∶ K[x] I Λ → (K[x] I Λ )
* through the commuting diagram:

K[x] (I Λ ) ⊥ (K[x] I Λ ) * K[x] I Λ Ĥ ≅ H
If the rank of Ĥ is finite and equal to r, then the dimension of K[x] I Λ and (K[x] I Λ ) * , as K-vector spaces, is r and the injective linear operator H is then an isomorphism. When B ⊂ K[x] is a basis for K[x] I Λ , we can identify K[x] I Λ with ⟨B⟩ and thus H B 1 is the matrix of H in this basis and its dual basis. It follows that H B 1 is invertible. The converse assertion is also true. It gives us a criterion for

B ⊂ K[x] to define a basis of K[x] I Λ . Theorem 3.5 Assume that rank Ĥ = r < ∞ and consider B = {b 1 , . . . , b r } ⊂ K[x]. Then B is a basis of K[x] I Λ if and only if the matrix H B 1 is invertible. proof: If H B 1 is invertible then B is linearly independent modulo I Λ and thus B is a basis of the r- dimensional vector space K[x] I Λ . Indeed, assume c = (c 1 , . . . , c r ) ∈ K r is such that c 1 b 1 + ⋅ ⋅ ⋅ + c r b r ≡ 0 mod I Λ . Using the definition of I Λ , we get c 1 Λ(b 1 b i ) + ⋅ ⋅ ⋅ + c r Λ(b r b i ) = 0, ∀i = 1, . . . , r. These equalities amount to H B 1 c = 0. When K[x] = ⟨B⟩ ⊕ I Λ , i.e. when B is a basis of K[x] I Λ , then, with the natural identification of ⟨B⟩ with K[x] I Λ , H B 1 is the matrix of H. Since H is an isomorphism, H B 1 must be invertible. ∎ Multiplication maps in K[x] I Λ
We now assume that the Hankel operator Ĥ associated to Λ has finite rank r. Then K[x] I Λ is of dimension r when considered as a linear space over K. In this paragraph we make explicit how to obtain the matrices of the multiplication maps from the matrices of associated Hankel operators.

For p ∈ K[x], consider the multiplication map

M p ∶ K[x] → K[x] q ↦ p q and M p ∶ K[x] I Λ → K[x] I Λ [q] ↦ [p q].
M p is a well defined linear map respecting the following commuting diagram

K[x] K[x] K[x] I Λ K[x] I Λ M p M p
Let us temporarily introduce the Hankel operator Ĥp associated to Λ p . Note that Ĥp = Ĥ ○ M p . Therefore the image of Ĥp is included in the image of Ĥ and ker Ĥ ⊂ ker Ĥp . We can thus construct the map

H p ∶ K[x] I Λ → (K[x] I Λ )
* that satisfies H p = H ○ M p and the following commuting diagram.

K[x] (I Λ ) ⊥ (K[x] I Λ ) * K[x] I Λ K[x] K[x] I Λ Ĥp ≅ H p M p M p
H B p = H B 1 M B p
where

H B 1 = (Λ(b i b j )) 1≤i,j≤r and H B p = (Λ(p b i b j )) 1≤i,j≤r
proof: The matrix of H p in B and its dual basis is H B p while H p = H ○ M p . ∎ Combining the above with Theorem 3.2 and Theorem 3.3 leads then to the following result.

Corollary 3.7 Let H ∶ K[x] I Λ → K[x] I Λ be the Hankel operator associated to Λ = ∑ r i=1 a i e ξi . Consider { 1 , . . . , r } ⊂ K[x]
a Lagrange basis associated to ξ 1 , . . . , ξ r . For any p ∈ K[x], i is a generalized eigenvector of (H p , H) associated to p(ξ i ).

Positivity

Associating a Hankel operator Ĥ ∶ K[x] → K[x] * to a linear form Λ ∶ R[x] → R is somewhat equivalent to associating to Λ the bilinear form (p, q) ↦ Λ(pq) = Ĥ(p)(q). When this associated bilinear form is positive semi-definite, that is Λ(p 2 ) = Ĥ(p)(p) ≥ 0 for all p ∈ R[x], then Λ and Ĥ are also said to be positive semi-definite. It is easily seen to be the case when Λ is the Hankel operator associated the linear form Ω ∶ p ↦ ∫ D p(x)dx or Λ = ∑ r i=1 a i e ξj , with ξ i ∈ R n and a i > 0. Note that the matrices H B 1 introduced in Section 3.1 are positive semi-definite when Ĥ has this property.

Theorem 3.8 If Λ ∶ R[x] → R is a positive semi-definite linear form then I Λ = {p Λ(pq) = 0, ∀q ∈ K[x]} is a real radical ideal.
What real radical means in general is to be found in [START_REF] Bochnak | Real algebraic geometry[END_REF]. Of relevance in this article is the fact that a zero dimensional real radical ideal of R[x] has all its roots in R n , as opposed to C n , and they are simple roots. Based on fundamental results in real algebraic geometry, the above assertion is proved by an elegant and short argument in [START_REF] Laurent | A generalized flat extension theorem for moment matrices[END_REF]Lemma 5.2].

A cubature is thus well characterized by its Hankel operator by combining Theorem 3.8 and Theorem 3.4. Theorem 3.9 Let Ĥ be the Hankel operator associated to the linear form

Λ ∶ R[x] → R[x]. Then Λ = ∑ r i=1
a i e ξi for some a i > 0 and ξ i ∈ R n if and only if Ĥ is positive semi-definite of finite rank r. When this is the case, ξ 1 , . . . , ξ r are the roots of the kernel I Λ of Ĥ.

The algorithm we shall develop in Section 5 establishes the existence of a cubature thanks to the above criterion before it goes on to compute the nodes. The following property specific to positive semi-definite forms allows to restrict the search space of the algorithm. Proposition 3.10 Assume that the linear form Λ ∶ R[x] → R is positive semi-definite and let Ĥ be its associated Hankel operator. For each δ ∈ N consider furthermore

Ĥ(δ) ∶ R[x] ≤δ → R[x] * ≤δ p ↦ Λ (δ) p where Λ (δ) p ∶ R[x] ≤δ → R q ↦ Λ(pq)
.

Then ker Ĥ(δ) = ker Ĥ ∩ R[x] ≤δ .
Furthermore rank Ĥ(δ) = rank Ĥ(δ+1) implies that Ĥ is of finite rank r = rank Ĥ(δ) .

proof: The inclusion ker Ĥ ∩ R[x] ≤δ ⊂ ker Ĥ(δ) is easy. The converse inclusion follows from the semidefiniteness of the bilinear form (p, q) ↦ Λ(pq). If p ∈ ker Ĥ(δ) then in particular Λ(p 2 ) = 0. The Cauchy-Schwarz inequality (Λ(pq) 2 ≤ Λ(p 2 )Λ(q 2 )) implies then that Λ(pq) = 0 for all q ∈ R[x]. Hence p ∈ ker Ĥ.

It follows that ker

Ĥ(δ) = ker Ĥ(δ+1) ∩ R[x] ≤δ . Thus if rank Ĥ(δ) = rank Ĥ(δ+1) , all the element of R[x] ≤δ+1 ∖ R[x]
≤δ are linearly dependent on the elements of R[x] ≤δ modulo ker Ĥ(δ+1) . In particular for each x α with α = δ + 1 there exists p α ∈ K[x] ≤δ such that x α + p α ∈ ker Ĥ(δ+1) ⊂ ker Ĥ. The latter is an ideal hence, for all 1 ≤ i ≤ n, the polynomial x i x α +x i p α of degree δ+2 belongs to ker Ĥ, and thus to ker Ĥ(δ+2) = ker Ĥ∩R[x] ≤δ+2 . It follows that rank Ĥ(δ+2) = rank Ĥ(δ+1) = rank Ĥ(δ) . The claim follows by induction. ∎

Observe that in general rank Ĥ(0) ≤ rank Ĥ(1) ≤ . . . ≤ rank Ĥ(δ) ≤ rank Ĥ(δ+1) ≤ . . . and this sequence of ranks becomes stationary when Ĥ is of finite rank. Yet you may have Ĥ(δ) = rank Ĥ(δ+1) < rank Ĥ(δ+1) when Ĥ is not positive semi-definite.

The proof above actually points to the fact that a basis for ker Ĥ(δ+1) provides a set of generators for I Λ (cf. also [START_REF] Laurent | Sums of squares, moment matrices and optimization over polynomials[END_REF]Theorem 5.19]). We do not dwell on this as we wish to emphasize that the nodes of the cubature are more straightforwardly obtained as generalized eigenvalues, a problem for which that are well studied numerical algorithms.

Before computing the nodes, we shall first determine if the characterized cubature is inside. Assuming that the support K of the cubature is given as a semi algebraic set, i.e. K = {x ∈ R n g 1 (x) ≥ 0, . . . , g s (x) ≥ 0}, for some polynomials g 1 , . . . , g s ∈ R[x], the following criterion allows one to test if the cubature is inside.

Corollary 3.11 Consider Λ = ∑ r i=1 a i e ξi , with ξ i ∈ R n and a i > 0, and q ∈ R[x]. Then q(ξ i ) ≥ 0, for 1 ≤ i ≤ r, if and only if the Hankel operator Ĥq associated to Λ q is positive semi-definite.

proof: As Λ q = ∑ r i=1 a i q(ξ i )e ξi , by Theorem 3.9, Ĥq is positive semi-definite if and only if a i q(ξ i ) ≥ 0 for all 1 ≤ i ≤ r. ∎

Flat extensions

In determining a cubature of degree 2d for a linear form Ω ∶ p ↦ ∫ D p(x)dx one shall seek to extend the linear form Λ (d) = Ω R[x] ≤2d to a linear form Λ on R[x] such that the associated Hankel operator is positive semi-definite of finite rank. The existence of such a form needs to be verified in some finite way though.

A linear form Λ (δ+κ) on R[x] ≤2δ+2κ is an extension of a given linear form Λ

(δ) on R[x] ≤2δ if its restriction to R[x] ≤2δ is Λ (δ) , that is if Λ (δ+κ) (p) = Λ (δ) (p) ∀p ∈ R[x] ≤2δ .
Let Λ (δ) be a linear form on K[x] ≤2δ . Similarly to Section 3.1, we associate with Λ (δ) the Hankel operator

H (δ) ∶ K[x] ≤δ → K[x] * ≤δ p ↦ Λ (δ) p
, where Λ

(δ) p ∶ K[x] ≤δ → K q ↦ Λ(p q).
A linear form Λ (δ+κ) on R[x] ≤2δ+2κ is a flat extension of a given linear form Λ (δ) if furthermore the rank of its associated Hankel operator H (δ+κ) is the rank of the Hankel operator H (δ) associated with the linear form Λ (δ) . In this case, H (δ+κ) is positive semi-definite if and only if H (δ) is positive semi-definite.

The flat extension theorem below was proved in terms of moment matrices, i.e. matrices of Ĥ(δ) in the monomial basis of K[x] ≤δ [START_REF] Curto | Solution of the truncated complex moment problem for flat data[END_REF][START_REF] Curto | The truncated complex K-moment problem[END_REF][START_REF] Laurent | Revisiting two theorems of Curto and Fialkow on moment matrices[END_REF][START_REF] Laurent | Sums of squares, moment matrices and optimization over polynomials[END_REF][START_REF] Laurent | A generalized flat extension theorem for moment matrices[END_REF]. It can be interpreted as an involutivity or regularity criterion [START_REF] Malgrange | Cartan involutiveness = Mumford regularity[END_REF][START_REF] Malgrange | Systèmes différentiels involutifs, volume 19 of Panoramas et Synthèses [Panoramas and Syntheses[END_REF], and was generalized as such in [START_REF] Ma | A certificate for semidefinite relaxations in computing positive-dimensional real radical ideals[END_REF].

Theorem 3.12 (Flat extension theorem) Let Λ (δ) be a linear form on R[x] ≤2δ . Assume that Λ (δ) is a flat extension of its restriction to R[x] ≤2δ-2 . Then there exists a unique flat extension of

Λ (δ) to R[x] ≤2δ+2κ for all κ ≥ 1.
4 Fraction-free diagonalization of positive semidefinite matrices

In this section we consider a symmetric matrix whose entries are polynomials in some parameters h 1 , . . . , h t over a computable field extension K of Q is. We present an algorithm in symbolic computation that computes a partition of the locus of positive semi-definiteness of the matrix in R t . Each region in this partition is given by a set of equalities and inequalities together with a parametrization of the diagonalization valid on this region of the parameter space. Additional specifications for this algorithm provides the central computational ingredient for determining the existence of a cubature with given degree and number of nodes.

While the problem can be approached with Gaussian elimination it is preferable to provide a fraction free approach that limits the expression growth, a key difficulty in symbolic computation. Bareiss' algorithm for the triangularization of a given matrix [START_REF] Bareiss | Sylvester's identity and multistep integer-preserving Gaussian elimination[END_REF][START_REF] Bareiss | Computational solutions of matrix problems over an integral domain[END_REF][START_REF] Geddes | Algorithms for computer algebra[END_REF] is the model for such an approach, one that does not require the computation of greatest common divisors. It was revisited for the diagonalization of a given symmetric matrix [START_REF] Hong | Exact division applied to diagonalizing symmetric matrices[END_REF][START_REF] Rouillier | Formules de Bareiss et réduction de formes quadratiques[END_REF]. An algorithmic difficulty in the description and implementation of the algorithm is the necessity of pivoting the rows and columns when the next diagonal element is tested to be zero. We observe, and this appears to be new, that the pivoting is unnecessary when we restrict to positive semi-definite symmetric matrices.

Positive semi-definite matrices

A symmetric matrix A = (a ij ) 1≤i,j≤c , with entries in R, is called positive semi-definite if v t Av ≥ 0, ∀v ∈ R c . It is called positive definite if it satisfies furthermore v t Av = 0 ⇔ v = 0. Given any matrix A = (a ij ) 1≤i,j≤c , the determinant of its submatrix A J = (a ij ) (i,j)∈J 2 with J ⊂ {1, .
. . , c} is called a principal minor. These minors provide a characterization of positivity for symmetric matrices [74, Theorem 7.2].

Theorem 4.1 Let A be a symmetric matrix with entries in R. Then

• A is positive semi-definite if and only if every principal minor of A is nonnegative.

• A is positive definite if and only if every leading principal minor of A is positive.

Instead of requiring all minors to be nonnegative, we can consider c specific linear combinations of the minors, namely the coefficients of the characteristic polynomial [START_REF] Basu | Algorithms in real algebraic geometry[END_REF]Chapter 8.2.4].

When the matrix A depends polynomially on some parameters the locus of positivity is thus a semi algebraic set. Describing semi-algebraic sets is a central question in computational real algebraic geometry [START_REF] Basu | Algorithms in real algebraic geometry[END_REF]. When one is interested in finding a (single) point in this semi-algebraic set, SemiDefinite Program solvers can provide a numerical solution. Exact algorithms were also developed to decide whether this particular semialgebraic set is empty or not and, in the negative case, exhibits an algebraic representation of a point in this semi-algebraic set. We refer to [START_REF] Henrion | Exact algorithms for linear matrix inequalities[END_REF] and the reference therein. We shall develop an approach that provides the basis in which the associated bilinear form is diagonal in addition to a description of the semi-algebraic description of the locus of positivity. As a side note, the Maple implementation of our approach actually succeeded with Section 6 when the Maple implementation of [START_REF] Henrion | Exact algorithms for linear matrix inequalities[END_REF] did not.

We think it worth sharing the original observation on Gauss algorithm that prompted our approach. Gaussian elimination is a staple algorithm to compute, among other things, the rank of a matrix. Proposition 4.2 below basically offers a simplification of this algorithm that decides if a symmetric matrix is positive semi-definite. It is based on an observation that shows that there is no need to permute the columns when a pivot is zero. Proposition 4.2 Let A = (a ij ) 1≤i,j≤c be a symmetric matrix with entries in R.

1. If A is a 1 × 1 matrix, then A is positive semi-definite if and only if a 11 ≥ 0.

Otherwise:

(a) If a 11 < 0, then A is not positive semi-definite.

(b) If a 11 = 0, then A is positive semi-definite if and only if a 1j = 0 for every j = 1, . . . , c and the submatrix obtained by deleting the first row and the first column is positive semi-definite.

(c) If a 11 > 0, then for each i > 1 subtract ai1 a11 times row 1 from row i and delete the first row and the first column. Then A is positive semi-definite if and only if the resulting matrix is positive semi-definite.

proof: Since the coefficient a 11 is a principal minor, a 11 ≥ 0 is a necessary condition for the positive semi-definiteness of the matrix A. In the particular case of a 1 × 1 matrix A, this is also a sufficient condition by Theorem 4.1. It remains now to study the cases 2.(b) and 2.(c).

Assuming that a 11 = 0, then we distinguish two cases:

• If a 1j = 0 for all j = 1, . . . , c, then every principal minor det A J with 1 ∈ J is zero since the coefficients of a row of the matrix A J are zero. As a consequence, A is positive semi-definite if and only if the matrix (a ij ) 2≤i,j≤c is positive semi-definite.

• If there is a nonzero coefficient in the first row of A, that is if there exists j such that a 1j ≠ 0, then the principal minor det a 11 a 1j a j1 a jj is negative since det a 11 a 1j a j1 a jj = a 11 a jja 1j a j1 = -a 2 1j < 0.

Then A is not positive semi-definite.

Assuming that a 11 > 0, the row operations describe the matrix equality A = LU , that is

⎛ ⎜ ⎝ a 11 . . . a 1c ⋮ ⋮ a c1 . . . a cc ⎞ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 1 0 . . . . . . 0 a21 a11 1 ⋱ ⋮ ⋮ 0 ⋱ ⋱ ⋮ ⋮ ⋮ ⋱ ⋱ 0 ac1 a11 0 . . . 0 1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ a 11 a 12 . . . . . . a 1c 0 ã22 . . . . . . ã2c ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ 0 ãc2 . . . . . . ãcc ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ with ãij = a ij -ai1 a11 a 1j for all (i, j) ∈ N 2 such that 2 ≤ i, j ≤ c. Since we have ãij = a ij - a i1 a 11 a 1j = a ji - a 1i a 11 a j1 = ãji ∀2 ≤ i, j ≤ c,
the matrix à = (ã ij ) 2≤i,j≤c is a symmetric matrix.

For every J ⊂ {1, . . . , c} with 1 ∈ J, we also have det A J = det L J det U J and det L J = 1 so that the principal minor det A J of A is the principal minor det U J of U . We have furthermore that det U J = a 11 det U J∖{1} . Since a 11 is positive, the principal minors of U are nonnegative if and only if the ones of the matrix U without the first row and the first column are nonnegative. ∎ When the matrix A depends polynomially on some parameters, one can add a case distinction at each step to provide the conditions on the parameters for the matrix to be positive semidefinite and compute the rank. However the algorithm suggested by Proposition 4.2 runs then with rational functions whose numerators and denominators grow fast in degree. We shall thus introduce a fraction-free variant to this algorithm, elaborating on Bareiss elimination algorithm.

Fraction-free triangularization and diagonalization

We seek to perform elimination in matrices the entries of which lie in a polynomial ring R = K[h 1 , . . . , h t ] where Q ⊂ K ⊂ R is a field that allows to perform exact arithmetic operations. Typically K is an algebraic extension of Q as for instance

K = Q[ √ 3].
The fraction-free approach of Bareiss [START_REF] Bareiss | Sylvester's identity and multistep integer-preserving Gaussian elimination[END_REF][START_REF] Bareiss | Computational solutions of matrix problems over an integral domain[END_REF][START_REF] Geddes | Algorithms for computer algebra[END_REF] to Gaussian elimination limits the growth of the entries by dividing a known extraneous factor at each step. The intermediate results can furthermore be expressed in terms of minors of the original matrix. This gives a clear idea on the growth and the specialization property of the elimination obtained. The scheme was revisited in [START_REF] Hong | Exact division applied to diagonalizing symmetric matrices[END_REF][START_REF] Rouillier | Formules de Bareiss et réduction de formes quadratiques[END_REF] for the diagonalization of symmetric matrices. The essence of these algorithms lie in the following results. The first statement is the heart of Bareiss's algorithm [START_REF] Bareiss | Sylvester's identity and multistep integer-preserving Gaussian elimination[END_REF][START_REF] Bareiss | Computational solutions of matrix problems over an integral domain[END_REF][START_REF] Geddes | Algorithms for computer algebra[END_REF]. The second observation was presented in [START_REF] Hong | Exact division applied to diagonalizing symmetric matrices[END_REF][START_REF] Rouillier | Formules de Bareiss et réduction de formes quadratiques[END_REF]. 0,0 = 1 and for an integer such that 1 ≤ ≤ c -1, let

A ( ) = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ a (0) 1,1 a (0) 1,2 ⋯ ⋯ ⋯ ⋯ a (0) 1,c 0 a (1) 2,2 a (1) 2,c ⋮ ⋱ ⋱ ⋮ ⋮ ⋱ a ( -1) , ⋯ ⋯ a ( -1) ,c ⋮ 0 a ( ) +1, +1 ⋯ a ( ) +1,c ⋮ ⋮ ⋮ ⋮ 0 ⋯ ⋯ 0 a ( ) c, +1 ⋯ a ( ) c,c ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ with a ( ) i,j = det ⎛ ⎜ ⎜ ⎜ ⎝ a 1,1 ⋯ a 1, a 1,j ⋮ ⋮ ⋮ a ,1 ⋯ a , a ,j a i,1 ⋯ a i, a i,j ⎞ ⎟ ⎟ ⎟ ⎠ for < i, j ≤ c.
Let be an integer such that 1 ≤ ≤ c -1.

1. If a (k-2) k-1,k-1 ≠ 0 for all k = 1, . . . , , then A ( ) = L . . . L 1 A with L k = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 1 0 . . . . . . . . . . . . . . . 0 0 ⋱ ⋱ ⋮ ⋮ ⋱ ⋱ ⋱ ⋮ ⋮ 0 1 ⋱ ⋮ ⋮ ⋮ -l k+1 l k ⋱ ⋮ ⋮ ⋮ ⋮ 0 ⋱ ⋱ ⋮ ⋮ ⋮ ⋮ ⋱ ⋱ 0 0 . . . 0 -l c 0 . . . 0 l k ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ with l k = a (k-1) k,k a (k-2) k-1,k-1 and l j = a (k-1) j,k a (k-2) k-1,k-1 for k + 1 ≤ j ≤ c.

If

A is furthermore symmetric, then the matrix Ã( ) = L ( ) AL ( ) t , with L ( ) = L . . . L 1 , is also symmetric. This matrix Ã( ) is given by

Ã( ) = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ã(0) 11 0 ⋯ ⋯ ⋯ 0 0 ⋱ ⋮ ⋮ ã( -1) , 0 ⋯ 0 0 ã( ) +1, +1 ⋯ ã( ) +1,c ⋮ ⋮ ⋮ ⋮ 0 0 ã( ) c, +1 ⋯ ã( ) c,c ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ , with ã(k) i,j = a (k) i,j a (k-1) k,k . (7)
The multiplication on the left by the matrix L correspond to an elimination step : a multiple of the -th row is subtracted to a non zero multiple of each of the row below so as to have zeros appear under the diagonal in the -th column. The simplest choice of multiples would be a

( -1) i,
and a ( -1) ,

. But one proves [START_REF] Bareiss | Sylvester's identity and multistep integer-preserving Gaussian elimination[END_REF][START_REF] Bareiss | Computational solutions of matrix problems over an integral domain[END_REF][START_REF] Geddes | Algorithms for computer algebra[END_REF] that the so obtained entries are all divisible by the previous pivot a

( -2)
-1, -1 . It can thus be divided out. If we consider a symmetric c × c matrix A with entries in R that is positive definite then, according to Theorem 4.1, a ( -1) ,

> 0 for all 1 ≤ ≤ c. We can thus proceed with elimination steps without any pivoting and obtain a diagonal matrix Ã(c-1) . In general though one might encounter a pivot a

(k-1) k,k
that is zero. To go further in the elimination one then has to introduce a pivoting strategy, which differs for diagonalization [START_REF] Hong | Exact division applied to diagonalizing symmetric matrices[END_REF] compared to triagularization [START_REF] Bareiss | Computational solutions of matrix problems over an integral domain[END_REF][START_REF] Geddes | Algorithms for computer algebra[END_REF]. We are here concerned with positive semi-definite matrices. A property of the matrices A (l) similar to the one disclosed in Corollary 4.2 allows us to forgo the pivoting.

Locus of semi-positivity by diagonalization

Consider now the integral domain R to be the polynomial ring

K[h] = K[h 1 , . . . , h t ], where Q ⊂ K ⊂ R.
Our goal is to determine the locus, in R t , such that the symmetric matrix A with entries in K[h] becomes semi-definite under substitution. This locus is a semi-algebraic set : it is given as a union of sets, each of which is described by a finite set of polynomial equalities and inequalities.

To formalize our statements let us introduce a specialization φ ∶ K[h] → R to be a K-morphism. A typical example is given by

φ̵ h ∶ K[h] → R, p ↦ p( ̵ h) with ̵ h = ( ̵ h 1 , . . . , ̵ h t ) a point in R t . Given a matrix A = (a ij ) 1≤i,j≤c with a ij ∈ K[h],
we denote by φ(A) the matrix obtained by applying φ to the coefficients of A, that is

φ(A) = (φ(a ij )) 1≤i,j≤c = (ā ij ) 1≤i,j≤c . Likewise ā( ) ij = det ⎛ ⎜ ⎜ ⎜ ⎝ ā1,1 ⋯ ā1, ā1,j ⋮ ⋮ ⋮ ā ,1 ⋯ ā , ā ,j āi,1 ⋯ āi, āi,j ⎞ ⎟ ⎟ ⎟ ⎠ = φ a ( ) ij ∀1 ≤ ≤ c -1, ∀ < i, j ≤ c.
= K[h], Q ⊂ K ⊂ R, let φ ∶ R → R be a specialisation. Assume that φ a (k-1) k,k > 0 for 1 ≤ k < .
For φ(A) to be semi-definite positive it is necessary that φ a , +j = 0 for all 1 ≤ j ≤ cand the principal submatrix φ A ˆ , where A ˆ is obtained from A by removing the -th row and -th column, is positive semi-definite of rank r. proof: Let A J be the ( +2)×( +2) leading principal submatrix of A. If we apply the adequately truncated L ( ) to A J , we obtain the ( + 2) × ( + 2) leading principal submatrix of A ( ) . We have then a

( -1) , 2 -1 k=1 a (k-1) k,k det(A J ) = k=1 a (k-1) k,k a ( ) +1, +1 a ( ) +2, +2 -a ( ) +1, +2 2
.

If a

(k-1) kk > 0 for all 1 ≤ k ≤ and a ( ) +1, +1 = 0 while a

( -1)
+1, +2 ≠ 0 then det A J < 0. According to Theorem 4.1 this contradicts A being positive semi-definite.

Notice that we used the fact that

det L k = ⎛ ⎝ a (k-1) k,k a (k-2) k-1,k-1 ⎞ ⎠ c-k so that det L ( ) = a ( -1) , c- -1 k=1 a (k-1) k,k
and replaced c by + 2 to reflect the truncation.

The same argument, changing A J to be another principal matrix, shows that for A to be positive semi-definite when a

( ) +1, +1 = 0, it is required that a ( ) +1, +j = 0 for all 2 ≤ j ≤ c -. ∎ Notice that for k ≤ we have A ˆ (k) = A (k)
ˆ . So when we encounter a pivot that is zero, we ignore that row and column and continue on. The entries of the successive matrices we encounter are thus determinants of a principal submatrix of A, whose size is the rank of A.

We give a recursive description of our algorithm. We will use the following notation. For a matrix A = (a ij ) 1≤i,j≤c with entries in R and p ∈ R, we define the first elimination matrix with the last non-zero pivot p as Z,E], where

L(A, p) = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 1 0 . . . . . . 0 - a 2,1 p a 1,1 p ⋱ ⋮ ⋮ 0 ⋱ ⋱ ⋮ ⋮ ⋮ ⋱ ⋱ 0 - a c,1 p 0 . . . 0 a 1,1 p ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ Algorithm 4.5 Diagonalization & Positivity Input : ▷ A = (a ij ) 1≤i,
• P is a list of polynomials in K[h], • Z is a set of polynomials in K[h],
• E is a list of elimination matrices in K(h) whose denominators are power products of elements in P .

The set S satisfies the following property: For a specialization φ ∶ K[h] → R, φ(A) is positive semi-definite if and only if there is a triplet [P , Z, E] in S such that φ(p) > 0 ∀p ∈ P and φ(q) = 0 ∀q ∈ Z.

In this case Card(P ) -1 is the rank of φ(A) and, letting L be the product of the elements of E, φ(LAL t ) is a diagonal matrix whose non zero entries are the product of two consecutive elements in φ(P ) = {φ(p), p ∈ P }.

If c = 0 then return [P , Z, E]; Otherwise:

If a 11 ∉ K or a 11 > 0 then

• Let p be the last element of P .

• Let A 1 be the submatrix of L(A, p) A obtained by removing the first row and first column.

• Append a 11 to P to give P 1 .

• Append L(A, p) to E to give E 1 .

• S 1 := Diagonalization & Positivity ( A 1 , P 1 , Z, E 1 , T 1 ).
If a 11 ∉ K or a 11 = 0 then

• Z 2 := {a 1,j 1 ≤ j ≤ n} ∪ Z.
• A 2 is obtained from A by removing the first row and first column.

• S 2 := Diagonalization & Positivity ( A 2 , P , Z 2 , E, T 2 ). return S ∶= S 1 ∪ S 2 .
The elements in P and Z are principal minors of A of size up to order Card(P ) -1. This determines their degrees in terms of the degrees of the entries of A. The output set S can have as many as 2 c triplets. One can lower this number if we can dismiss certain branches by checking if the semi-algebraic set defined by a pair (P, Z) is consistent. Though algorithmic [START_REF] Basu | Algorithms in real algebraic geometry[END_REF], this is no light task. In its application to cubature, the algorithm will be made more specific, which shall result in a restricted output.

Diagonalization for determining the existence of cubatures

Taking into account the specificity of the cubature problem, we revisit the procedure in the previous section.

There are less cases to be distinguished so that the size of the output is smaller.

The cubature problem comes with not only a c × c symmetric matrix A whose entries belong to K[h] = K[h 1 , . . . , h t ] but also three input indices c ′ , c ′′ and r, which satisfy c ′ ≤ r ≤ c ′′ ≤ c. With these indices are considered two leading principal submatrices of A, A ′ and A ′′ , respectively of size c ′ and c ′′ . The entries of A are of degree 1 or less in general, but the entries of A ′ more specifically belong to K, and A ′ is positive definite. We are now looking for the locus in R t where the respective specializations Ā′′ and Ā of A ′′ and A are semi-definite positive of rank r -extending the terminology introduced in Section 3.5, Ā is a flat extension of Ā′′ . Algorithm 4.5 can then be appropriately modified to determine this locus. The output size is then reduced to c ′′ -c ′ r-c ′ triplets :

• Since A ′ is positive definite, the algorithm runs without splitting as long as the matrix is of size bigger than cc ′ . The first c ′ pivots are added to P .

• Since A ′′ is required to have the same rank as A, the algorithm is stopped when the matrix is of size cc ′′ . All its entries are added to Z.

• Only the branches that can terminate with Card P = r + 1 (as we included 1 as first pivot in P ) need to be continued.

The changes are straightforward. We simply provide the formal description of the input and output so as to serve as reference when this algorithm is called upon.

Algorithm 4.6 Diagonalization & Positivity with Rank Constraints

Input :

▷ Integers c ′ ≤ r ≤ c ′′ ≤ c. ▷ A c × c symmetric matrix A with entries in K[h] ≤1 , whose c ′ × c ′ leading principal submatrix is positive definite with entries in K. ▷ P a list of polynomials in K[h]. ▷ Z a set of polynomials in K[h]. ▷ E a list of elimination matrices. Output: ▷ A set S of c ′′ -c ′ r-c ′ triplets [P , Z, E], where • P is a list of polynomials in K[h] with Card P = r + 1, • Z is a set of polynomials in K[h],
• E is a list of elimination matrices in K(h) whose denominators are power products of elements in P .

The set S satisfies the following property: For a specialization φ ∶ K[h] → R, φ(A) and φ(A ′′ ), where A ′′ is the leading principal submatrix of size c ′′ , are positive semi-definite of rank r if and only if there is a triplet [P, Z, E] such that φ(p) > 0 ∀p ∈ P and φ(q) = 0 ∀q ∈ Z.

In this case, letting L be the product of the elements of E, φ(LAL t ) is a diagonal matrix whose non zero entries are the product of two consecutive elements in φ(P ) = {φ(p), p ∈ P }.

Algorithms to determine cubatures

This section offer the synthesis of our approach to determine cubatures of degree d for an integral operator.

Our approach is two fold. We first parametrize the cubatures with a given number of nodes, hence establishing their existence, or inexistence. We then provide the means to compute their nodes and weights. A salient feature in this approach is that the cubatures of a given degree d and with a given number of nodes are described exhaustively.

All along this section, the polynomial ring K[x] is endowed with a basis M such that, for any δ ∈ N,

M (δ) = M∩K[x]
≤δ is a basis of K[x] ≤δ . Any basis consisting of homogeneous polynomials fits this description, as for instance the monomial basis. But so does a set of orthogonal polynomials for a given measure. We assume that M = {m 1 , m 2 , . . .} with M (δ) = {m 1 , m 2 , . . . , m r δ }, where r δ = dim K[x] ≤δ .

We seek to determine the cubatures of degree d for a positive linear form

Ω ∶ R[x] → R
i.e. we look for the linear forms

Λ ∶ R[x] → R p ↦ ∑ r i=1 a i p(ξ i ) , for some r ∈ N, a i > 0 and ξ i ∈ R n , such that Λ(m i ) = Ω(m i ), for all 1 ≤ i ≤ dim R[x] ≤d .
One is either given r or one might be seeking the minimal possible r. Our emphasis is to describe the algorithmic aspects of determining cubatures and will not dwell on the finest hypotheses under which the existence is secured, or on the sharper bound known on the number of nodes. We are mostly considering the forms Ω defined as Ω(p) = ∫ R n pdµ where µ is a positive Borel measure with compact support. It is also reasonable to assume that Ω is definite, i.e. Ω(p 2 ) > 0 ⇒ p = 0, at least when deg p ≤ d.

In Section 3 we reviewed the conditions for a linear form to be a cubature. These conditions translate as rank and positivity constraints on the associated Hankel operator. For the specific problem of determining an interpolatory cubature of degree d for a integral operator Ω, these conditions can be checked in a finite dimension associated with an integer δ ≤ d + 1. This is explained in Section 5.1. In Section 5.2 we thus parametrize the space of linear forms on R[x] ≤2δ whose restrictions to R[x] ≤d coincide with Ω and apply the algorithm established in Section 4 to determine the locus of cubatures with r nodes among them. Following additional results in Section 3 we then show how, for each characterized cubature, we can retrieve its nodes as eigenvectors in Section 5.4. The weights are then obtained as the solution of a nonsingular linear system. The whole process is demonstrated for the minimal cubatures of degree 5 of the hexagon in Section 6.

Upper bound on the minimal number of nodes

The question of existence of a cubature of any given degree and a bound for the minimal number of nodes was provided by Tchakaloff [START_REF] Tchakaloff | Formules de cubatures mécaniques à coefficients non négatifs[END_REF] [70, Section 3.3].

Theorem 5.1 Let µ be a positive measure with compact support K in R n and d a fixed integer. For some r ≤ dim R[x] ≤d , there exist r points ξ 1 , . . . , ξ r ∈ K, and r positive real numbers a 1 , . . . , a r such that

Ω(p) = r i=1 a i p(ξ i ), ∀p ∈ R[x] ≤d .
The above theorem was generalized in [START_REF] Putinar | A note on Tchakaloff's theorem[END_REF][START_REF] Bayer | The proof of Tchakaloff's theorem[END_REF] by reducing the hypothesis on the measure µ. The theorem remains true for any finite dimensional spaces of polynomials. It actually generalizes to spaces of functions other than polynomials [START_REF] Richter | Parameterfreie Abschätzung und Realisierung von Erwartungswerten[END_REF]. See [START_REF] Schmüdgen | The moment problem[END_REF]Theorem 1.2.4]. These non contructive proofs call on Carathéodory's theorem in convex geometry.

Cubatures whose weights are obtained uniquely from the nodes are called interpolatory. Since the weights are solution of the linear system

W (d) ξ ⎛ ⎜ ⎝ a 1 ⋮ a r ⎞ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎝ Ω(m 1 ) ⋮ ⋮ Ω(m r d ) ⎞ ⎟ ⎟ ⎟ ⎠ , where W (d) ξ = (m j (ξ i )) 1≤i≤r1≤j≤r δ a cubature is interpolatory if the rank w (d) ξ of W (d) ξ
is equal to the number r of nodes. The following result shows that a subset of the nodes of a non-interpolatory cubature can be used as the set of nodes of an interpolatory one [14, Section 6.1]. In the litterature on computing cubatures, this was observed in [START_REF] Davis | A construction of nonnegative approximate quadratures[END_REF], [START_REF] Stroud | Approximate calculation of multiple integrals[END_REF]. It is also a consequence of Carathéodory's theorem [START_REF] Rockafellar | Convex analysis[END_REF]Corollary 17.1.2].

Lemma 5.2 Assume Λ = ∑ s i=1 a i e ξi is a cubature of degree d for Ω ∶ R[x] → R. Then there is a cubature with w (d) ξ nodes in {ξ 1 , . . . , ξ s }.
A corollary to Theorem 5.1 is thus the existence of an inside interpolatory cubature of degree d with at most dim R[x] ≤d nodes. Furthermore any minimal cubature needs to be interpolatory.

Algorithm to determine the existence of a cubature with r nodes

According to Theorem 3.9 and Theorem 3.12, a cubature for Ω of degree d with r nodes exists if we can find δ ∈ N and a linear form

Λ (δ) ∶ R[x] ≤2δ → R such that (1) Λ (δ) (p) = Ω(p), for all p ∈ R[x] ≤d ;
(2) Λ (δ) is positive semi-definite and its associated Hankel operator

Ĥ(δ) ∶ R[x] ≤2δ → R[x] *
≤2δ is of rank r;

(3) Λ (δ) is a flat extension of its resriction Λ (δ-1) to R[x] ≤2δ-2 . In other words Ĥ(δ) has the same rank as the Hankel operator

Ĥ(δ-1) ∶ R[x] ≤δ-1 → (R[x] ≤δ-1 ) * associated to Λ (δ-1) .
The cubature is then the unique flat extension Λ ∶ R[x] → R of Λ (δ) predicted by Theorem 3.12. Theorem 3.9 ensures then that Λ = ∑ r i=1 a i e ξi for

• ξ 1 , . . . , ξ r ∈ R n , the roots of the ideal

I Λ = {p Λ(pq) = 0, ∀q ∈ K[x]};
• a i = Λ( i ) > 0, where { 1 , . . . , r } is a Lagrange basis for {ξ 1 , . . . , ξ r }.

Some inequalities need to be satisfied between d, δ and r for the above to make sense. First 2) and ( 3) can be checked computationally on the matrix

⌊ d 2 ⌋ < δ and r ≤ dim R[x] ≤δ-1 . Note also that the restriction to R[x] ≤⌊ d 2 ⌋ of the Hankel operator associated to Ω agrees with Ĥ(⌊ d 2 ⌋) . If Ω is definite positive then we must have dim R[x] ≤⌊ d 2 ⌋ ≤ r. Given Λ (δ) ∶ R[x] ≤2δ → R such that Property (1) stands, Properties (
H (δ) 1 = (Λ(m i m j )) 1≤i,j≤r δ of Ĥ(δ) in the basis M (δ) of R[x]
≤δ , and its dual. Namely:

H (δ) 1
is positive semi-definite of rank r and its leading principal submatrix of size r δ-1 is already of rank r. One can then see how Algorithm 4.5 applies in determining the existence of cubatures.

Let us introduce Λ(δ) ∶ K[x] ≤2δ → K[h r d +1 , . . . , h r 2δ ] the linear map such that Λ(m i ) = Ω(m i ) if i ≤ r d , i.e. deg(m i ) ≤ d, h i if i > r d , i.e. deg(m i ) > d.
and build the matrix H(δ)

1 = Λ(m i m j ) 1≤i,j≤r δ whose entries are polynomials in K[h] = K[h r d +1 , . . . , h r 2δ
] of degree at most one. The existence of the desired cubature boils down to determining values

̵ h = ( ̵ h r d +1 , . . . , ̵ h r 2δ ) ∈ R r 2δ -r d such that under the specialization φ ∶ K[h] → R f ↦ f ( ̵ h) the matrix φ H(δ)
1 is positive semi-definite of rank r and has its principal leading submatrix of size r δ-1 already of rank r. The polynomial equalities and inequalities on the unknowns h for these conditions to be fulfilled are given by the output of Algorithm 4.6 applied to A ∶=

H(δ) 1 , c ′ ∶= dim R[x] ⌊ d 2 ⌋ , c ′′ ∶= dim R[x] δ-1 , and c ∶= dim R[x] δ .
Note though that some branches of the algorithm can be dismissed at an early stage by taking into account Proposition 3.10.

The remaining issue is how to choose δ to be certain that any cubature Λ with r nodes is such that Ĥ(δ) is of rank r. Were we to look for all possible cubatures with r nodes, we would need to take δ = r. If we restrict to searching for interpolatory cubatures we can limit ourselves to δ = d + 1 thanks to Lemma 5.2 : With Theorem 3.3 and Theorem 3.9, the latter indeed implies that for an interpolatory cubature Λ = ∑ r i=1 a i e ξi a subset of M (d) needs to be a basis for R[x] I Λ , where I Λ is the ideal of {ξ 1 , . . . , ξ r }. It follows that Ĥ(d) is of rank r from Theorem 3.5. One can furthermore take into account Proposition 3.10: when looking for all interpolatory cubatures with r nodes, it is sufficient to take

δ = 1 + min dim K[x] ≤d , ⌊ d 2 ⌋ + r -dim K[x] ⌊ d 2 ⌋ . Algorithm 5.3 Existence-of-a-cubature Input : ▷ The degree d of the sought cubature. ▷ A number r of nodes, dim R[x] ≤⌊ d 2 ⌋ ≤ r ≤ dim R[x] ≤d . ▷ The values of the linear form Ω on the basis {m 1 , . . . , m r d } of R[x] ≤d . Output : A set S of triplets [P, Z, J] where ▷ P is a list of r polynomials in K[h r d +1 , . . . , h r 2δ ]; ▷ Z is a set of polynomials in K[h r d +1 , . . . , h r 2δ ]; ▷ J is an increasing sequence of r integers between 1 and dim R[x] δ-1 .
The set S satisfies the following two converse properties:

• For any ̵ h = ( ̵ h r d +1 , . . . , ̵ h r 2δ ) ∈ R r 2δ -r d such that there is a triplet [P, Z, J] in S with p( ̵ h) > 0 ∀p ∈ P and q( ̵ h) = 0 ∀q ∈ Z,
there exists a unique cubature Λ of degree d for Ω with r nodes satisfying

Λ(m j ) = Ω(m j ) if j ≤ r d ̵ h j if r d < j ≤ r 2δ
This cubature is interpolatory and Λ(m i m j ) i,j∈J is non singular.

• For a linear form Λ ∶ R[x] → R to be an interpolatory cubature of degree d for Ω with r nodes there must be a triplet [P, Z, J] in S such that the specialization

φ ∶ K[h r d +1 , . . . , h r 2δ ] → R h j ↦ Λ(m j )
satisfies φ(p) > 0 ∀p ∈ P and φ(q) = 0 ∀q ∈ Z.

1. Let δ ∶= 1 + min dim K[x] d , ⌊ d 2 ⌋ + r -dim K[x] ⌊ d 2 ⌋ .
2. Construct the matrix H(δ)

1 = Λ(m i m j ) 1≤i,j≤c where Λ ∶ K[x] ≤2δ → K[h r d +1 , . . . , h r 2δ ] is the linear map such that Λ(m i ) = Ω(m i ) if i ≤ r d , i.e. deg(m i ) ≤ d, h i if i > r d , i.e. deg(m i ) > d. 3. Call Algorithm 4.6 [Diagonalization & Positivity with Rank Constraints] with A ∶= H(δ) 1 , c ∶= dim R[x] ≤δ , c ′ ∶= dim K[x] ⌊ d 2 ⌋ , and c ′′ ∶= dim R[x] ≤δ-1 ;
4. For each triplets [P, Z, E] in the output of the above, consider J to be sequence of the indices of the nonzero diagonal elements appearing in P .

This algorithm gives us a modus operandi to determine the minimal number of nodes in the cubatures of degree d for Ω. As seen in Lemma 5.2 these need to be interpolatory cubatures. The above algorithm is run successively on an increasing r, starting with a lower bound. When Ω is positive definite, the most general lower bound is dim R .

[x] ≤⌊ d 2 ⌋ . Yet,
Hence H B p H B 1 -1 W B ξ T = W B ξ T D.
Most polynomials p of degree 1 separate the roots of I Λ , i.e. are such that p(ξ 1 ), . . . , p(ξ r ) are pairwise distinct. It is thus fair to assume that we chose p as such, by taking a random combination of the polynomials in M (1) = {m 1 , m 2 , . . . , m n+1 }. The associated eigenvectors are thus unique, up to a scalar factor.

For any computed eigenvector v of With the hypotheses on d and Ω, we have M (1) = {m 1 , . . . , m n+1 } ⊂ B. With the previous computation we have thus determined the vector m 2 (ξ i ) . . . m n+1 (ξ i ) T for each 1 ≤ i ≤ r. This latter vector provides us with the coordinates of each node ξ i , for 1 ≤ i ≤ r.

H B p H B 1 -1 there is thus a λ ∈ R ∖ {0} and a node ξ such that v = λ ⎛ ⎜ ⎝ b 1 (ξ) ⋮ b r (ξ) ⎞ ⎟ ⎠ . Since b 1 ∈ R,
Computing the weights by solving a non singular linear system

According to Theorem 3.3, the matrix W B ξ that we determined at the previous stage is invertible. The unique solution a ∈ R r to the linear system

W B ξ ⎛ ⎜ ⎝ a 1 ⋮ a r ⎞ ⎟ ⎠ = ⎛ ⎜ ⎝ Ω(b 1 ) ⋮ Ω(b r ) ⎞ ⎟ ⎠ (8)
thus provides us with the weights of the cubature Λ.

Minimal cubatures of degree 5 for the regular hexagon

A minimal cubature of degree 5 for the regular hexagon is described in [START_REF] Stroud | Approximate calculation of multiple integrals[END_REF] under the name H 2 ∶ 5 -1. Only two such cubatures were known [START_REF] Lyness | Quadrature rules for regions having regular hexagonal symmetry[END_REF][START_REF] Stroud | Approximate calculation of multiple integrals[END_REF], and both are D 6 -invariant, where D 6 is the group of symmetry of the hexagon; It consists of reflections and rotations. As a mean to illustrate the algorithms presented in this article we show that there is actually a one parameter family of minimal cubatures. They are C 6 invariant, where C 6 is the subgroup of D 6 consisting of rotations.

We first use the monomial basis and straightforwardly apply the algorithm presented in Section 5 to determine all the cubatures of degree 5 with 7 nodes for the hexagon. The C 6 -invariance reveals itself ther experimentally. We then illustrate the use of a symmetry adapted basis to exploit and ascertain this invariance. Computing symmetric cubatures is the subject of a forthcoming paper. We give here a foretaste of the results on this example. We exploit the obtained block diagonal structure of the matrices of the Hankel operators to alleviate the computational cost the algorithms and deduce properties of the cubature from the size of the blocks.

Let H be the regular hexagon whose vertices are given by (±1, 0), (± 1 2 , ± √

2 ) (see Figure 1(a)). Its moments of order less than or equal to 5 are

H 1dx = 3 √ 3 2 , H x 2 1 dx = H x 2 2 dx = 5 √ 3 16 H x 4 1 dx = H x 4 2 dx = 21 √ 3 160 , H x 2 1 x 2 2 dx = 7 √ 3 160
and zero otherwise.

(1, 0) 1 2 , √ 3 2 -1 2 , √ 3 2
(-1, 0) Note that in virtue of [54] [14,Theorem 8.3] there cannot be cubature of degree 5 with dim R[x] ≤2 = 6 nodes, the hexagon being centrally symmetric. This could also be uncovered computationally following the process described in Section 5. We thus look for a cubature with 7 nodes.

-1 2 , - √ 3 2 1 2 , - √ 3 2 • • • • • • (a) Vertices • (0, 0) √ 14

Existence

We consider first the monomial basis M = {1, x 1 , x 2 , x 2 1 , x 1 x 2 , x 2 2 , x 3 1 , . . .} ordered with the graded reverse lexicographic order. We proceed as described in Section 5.2 to determine the sets of equalities and inequalities that characterize the cubatures of degree 5 for H with 7 nodes.

1. δ ∶= 1 + min dim K[x, y] ≤5 , ⌊ 5 2 ⌋ + 7 -dim K[x, y] ≤⌊ 5 2 ⌋ = 4;
We thus consider the monomial basis M (4) of R[x] ≤4 ordered according to the graded reverse lexicographic order

M = {1, x 1 , x 2 , x 2 1 , x 1 x 2 , x 2 2 , x 3 1 , x 2 1 x 2 , x 1 x 2 2 , x 3 2 , x 4 1 , x 3 1 x 2 , x 2 1 x 2 2 , x 1 x 3 2 , x 4 2 }.
2. The parameterized moment matrix, with Λ(x i 1 x j 2 ) = √ 3 160 h ij for i + j > 5, is We consider now the triplet [P, Z, J] with J = {1, 2, 3, 4, 5, 6, 7}. The only non constant polynomial in P is a positive multiple of 50h 60 -441. By computing a Gröbner basis we determine that the solutions of the equalities in Z for which 50h 60 -441 ≠ 0 satisfy :

HM 1 = √ 3 160 ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝
⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ There are dim R[x] ≤8 -dim R[x] ≤5 = 24
h 24 = - 196

Weights and nodes

After examining the output of Algorithm 5.3 we know that the cubatures of degree 5 with 7 nodes for the hexagon are in one to one correspondance with the real numbers θ in ]π, π].

We now proceed to compute the nodes and weights of the associated cubatures for a selection of values of the parameter θ following Section 5.4. We first consider the two cases θ = 0 and θ = π where the computations can be done in exact arithmetic. We shall retrieve the two D 6 -invariant cubatures that already exist in the literature. We then compute numerically the nodes and weights for a set of of other values of θ.

The symmetric cubatures previously known: The case θ = 0 belongs to the triplet with J = {1, 2, 3, 4, 5, 6, 7}.

Hence B = 1, x 1 , x 2 , x 2 1 , x 1 x 2 , x 2 2 , x 3 1 is a basis of R[x] I Λ . For p = x 1 + 5x 2 we have H B 1 = √ 3 160 ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝
⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ and H B p = √ 3 160 ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ . A matrix of normalized eigenvectors for (M B p ) T = H B p H B 1
-1 is then given by the transpose of

W B ξ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 1 0 0 0 0 0 0 1 √ 14
⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ .
Since x 1 and x 2 are the second and third polynomials in B, the second and third columns of W B ξ are the coordinates of the nodes. The nodes of the cubature are thus: We thus recognize the cubature given in [START_REF] Stroud | Approximate calculation of multiple integrals[END_REF]. It is represented in golden yellow in Figure 2.

ξ 1 = (0, 0), ξ 2 = √ 14 5 , 0 , ξ 3 = - √ 14
The case θ = π belongs to the triplet with J = {1, 2, 3, 4, 5, 6, 8} so that C = {1,

x 1 , x 2 , x 2 1 , x 1 x 2 , x 2 2 , x 2 1 x 2 } is a basis of R[x] I Λ . For p = x 1 + 5x 2 we have H C 1 = √ 3 160 ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 240
⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ A matrix of normalized eigenvectors for (M C p ) T = H C p H C 1
-1 is given by the transpose of

W C ξ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 1 0 0 0 0 0 0 1 0 √ 14
⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ .
Since x 1 and x 2 are the second and third polynomials in C, the second and third columns of W C ξ are the coordinates of the nodes and thus a second D 6 -invariant cubature [START_REF] Lyness | Quadrature rules for regions having regular hexagonal symmetry[END_REF]. It is represented in purple in Figure 2. Table 1: Nodes of the cubatures (excepted (0, 0)) for selected values of θ.

ξ 1 = (0, 0), ξ 2 = 0, √ 14 5 , ξ 3 = 0, - √ 14
New cubatures: In the case θ ∈] -π, π[, B = 1, x 1 , x 2 , x 2 1 , x 1 x 2 , x 2 2 , x 3 1 is a basis of R[x] I Λ . For p = x 1 + 5x 2 we have H B 1 = √ 3 160 ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 240
0 0 h 60 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ and H B p = √ 3 160 ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0
For each selected value of the parameter θ ∈]-π, π[, we compute (numerically) the eigenvectors of the matrix

M B p T = H B p H B 1
-1 and deduce the coordinates of the nodes (see Table 1 and Figure 2). We note that (0, 0) is a node in all the so obtained cubatures. We can also observe that the other 6 nodes all lie on a circle and form a C 6 orbit. Furthermore solving the Vandermonde-like linear system (8), we get the same weights for all selected values of the parameter θ

⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ a 1 a 2 a 3 a 4 a 5 a 6 a 7 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0.6649837921916226 0.3221820698602823 0.3221820698602823 0.3221820698602823 0.3221820698602823 0.3221820698602823 0.3221820698602823 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ,
where the first weight is associated with the node (0, 0). It thus looks like the computed cubatures are C 6 -invariant. Two cubatures are represented in different shades of brown in Figure 2.

Figure 2: Four distinct minimal cubatures : the origin is a node for all the cubatures. The other nodes are represented with a distinct color for each distinct cubautre.

Alternative basis and C 6 -invariant cubatures

The observed symmetry in the minimal cubatures computed above implies a certain structure. It manifested itself with the presence of many zeros in the matrices of the Hankel operators. We shall present here how to exploit the symmetry in the computation by working in an alternative basis. First it bring out a block diagonal structure in the matrices of the Hankel operators. Furthermore we can collect qualitative information at early stages and restrict our attention to a single sub-block in the multiplication and Vandermonde matrices. The results we apply results will actually be presented and proved in a forthcoming article [START_REF] Collowald | Cubatures with symmetry[END_REF]. This section is to be understood as a teaser.

Each block pertains to an irreducible representation of D 6 . The first block correspond to the trivial representation and consists of invariants. The 5-th and 6-th blocks, as well as the 7-th and 8-th blocks, work in pairs as they relate to irreducible representations of degree 2. Note that the first and second blocks are C 6 -invariants. Pairing all the blocks above actually provide a symmetry adapted basis according to C 6 that consists of 4 blocks. A crucial property is that if a linear form Λ is D 6 -invariant then Λ(d ij ) = 0 when i ≥ 2. If it is C 6 -invariant then Λ(d ij ) = 0 when i ≥ 3. Here we shall look for a cubature that is C 6 invariant and thus assume that it is the case.

A basis of R[x] ≤4 is given by the union of the following four blocks, where we ordered the elements by degree: To ease following the argument we provide the degrees of these elements as follows: deg B = {0, 2, 4} ∪ {3, 3} ∪ {1, 1, 3, 3} ∪ {2, 2, 4, 4, 4, 4} .

B = {d
The computation gain appears as follow : with Λ(d ij ) = h ij for i = 1, 2 and deg d ij > 5 and Λ(d ij ) = 0 for 3 ≤ i ≤ 8, the parameterized matrix of the Hankel operator is block diagonal : The triplets [P, Z, J] now combine four blocks [P (i) , Z (i) , J (i) ] for 1 ≤ i ≤ 4. Two index sets J (i) are predetermined owing to the degrees appearing in the B (i) : J (1) = {1, 2} and J (4) = {1, 1}. There is a total of four possibilities for the triplets J. They are given by the pairs J (2) , J (3) taken to be either ({1}, {1, 2}), ({2}, {1, 2}), ({}, {1, 2, 3}), or ({}, {1, 2, 4}). The last two possibilities bring inconsistent systems of polynomial equations Z.

H1 = (3) 2 diag H(1) 1 , H (2) 1 , H(3) 1 , H(4)
0 h 14 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , H(4) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 7
The case J (2) , J (3) = ({2}, {1, 2})

The set of equations and inequations admit a single solution : among the blocks of the symmetry adapted basis already bears some qualititative information. First, as there are two polynomials in the first block we know that there are two orbits of nodes. As there is a total of 7 nodes, we deduce that it can only be an orbit of 6 nodes in addition to the origin. We can furthermore deduce from the number of basis elements in each block the types of the orbits of nodes. In this case the orbit with 6 nodes has 2 nodes on the x-axis. The nodes are thus completely determined if we know the radius of the circle on which they lie.

h 14 =
H B p is also block diagonal for any invariant polynomial p. We have The eigenvalues of M B p T are the evaluations of p at the nodes. We hence see that the origin is one of the nodes, while the other six nodes lie on a circle of radius The case J (2) , J (3) = ({1}, {1, 2})

H B 1 ∶= (3) diag
The set of equations and inequations admit a one dimensional family of solutions. They can be parametrized by θ ∈]0, 2 π[ and are then given by h 14 = 49 250 , h 15 = 49 500 (cos θ + 1), h 16 = 343 3125 , h 17 = 343 6250 (cos θ + 1), h 21 = 49 500 sin θ, h 22 = 343 6250 sin θ.

We can infer that there are two orbits, one being the origin and and the other an orbit of 6 nodes.

For θ = π we have h 21 = 0 and h 22 = 0 so that it determines a D 6 invariant cubature. From the distribution of the block elements in the basis, we deduce the type of the 6 node orbit : it has two nodes on the y-axis. Computing the nodes of the cubature determined by the parameter θ. We shall first determine the radius of the orbits of nodes and then one node per orbit. The eigenvalues of M p are thus 0 and 14 5 , the latter with a multiplicity of 6. This implies that the nodes other than the origin are all on a circle of radius where α is to be determined.

H C 1 = √ 3 diag
As can be observed above, all the values of an invariant p at the different orbits of nodes already appear as eigenvalues of the first block M To find the value of q at the nodes in terms of the parameter θ describing the cubature we consider The value of q at the nodes is thus either 0 or Determining the weights. Since we have a C 6 -invariant cubature, the nodes on the same orbit have equal weight. Here we have two orbits of nodes and we show that we only need to solve a 2 × 2 linear system to retrieve the two distinct weights.

H (1) q = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣
The vector of weights a = a 1 a 2 a 2 a 2 a 2 a 2 a 2 T is solution of W B ξ a = m where m = 3 After the existence of the invariant cubature was obtained, we have thus determined the nodes and weights by dealing only with c × c matrices, where c is the number of orbits of nodes.

  Consider the linear form Ω on R[x] defined byΩ ∶ R[x] → R, p ↦ b a p(x)ω(x)dx,where [a, b] ⊂ R is an interval (finite or infinite) and ω is a nonnegative weight function defined on this interval. It is typically assumed that this function is measurable and and non negative on [a, b] and even that b a p 2 (x)ω(x)dx = 0 ⇒ p(x) = 0. Furthermore, the moments are given by m k = b a x k ω(x)dx, for k ∈ N.

.

  The matrices, in the monomial basis, of their restrictions to R[x] ≤r-1 are H

Theorem 3 . 2

 32 is a well defined linear map, simply called the multiplication by p. Its eigenvalues consist of the values of the function p on V C (I)[START_REF] Cox | Using algebraic geometry[END_REF] Theorem 4.5]. We shall give a proof of this result by exhibiting the left eigenvectors of the matrix of M p . This is a hint that duality can shed further light on the structure of K[x] I in general[START_REF] Mourrain | Polynomial exponential decomposition from moments[END_REF]. Let I be a zero-dimensional ideal in K[x] and assume B = {b 1 , . . . , b r } is a basis of K[x] I. For p a polynomial in K[x] consider M B p the matrix of the multiplication operator M p in the basis B. Then for any ξ ∈ V C (I) the row vector b 1 (ξ) . . . b r (ξ) is a left eigenvector of M B p associated to the eigenvalue p(ξ).

  is entirely determined by its values on 1 , . . . , r . Hence Λ = ∑ r i=1 a i e ξj . Assume now that Λ = ∑ r j=1 a j e ξj with a j ∈ C ∖ {0} and ξ j ∈ C n pairwise distinct. Let I(ξ) be the ideal of polynomials vanishing at the points ξ 1 , . . . , ξ r . By the Strong Nullstellensatz [17, Chapter 1.4], the ideal I(ξ) is radical. Furthermore K[x] I(ξ) is of dimension r. We prove next that I Λ = I(ξ).

Proposition 3 . 6

 36 Assume the Hankel operator Ĥ associated to the linear form Λ has finite rank r. Let B = {b 1 , . . . , b r } be a basis of K[x] I Λ . Then the matrix M B p of the multiplication by p in K[x] I Λ satisfies

Theorem 4 . 3

 43 Consider A = (a ij ) 1≤i,j≤c a c × c matrix with entries in an integral domain R. Let A(0)

Theorem 4 . 4

 44 With the hypotheses and notations of Theorem 4.3 with A symmetric and R

( - 1

 1 ) , ≥ 0. If φ a ( -1) , = 0 then φ(A) is positive semidefinite of rank + r if and only if φ a ( -1)

  j≤c a c × c symmetric matrix with entries in K[h]. ▷ P a list of polynomials in K[h] -P stands for pivots or positive. The default value is [1]. ▷ Z a set of polynomials in K[h] -Z stands for zero. The default value is ∅. ▷ E a list of elimination matrices. The default value is the empty list []. Output: ▷ A set S of triplets [P ,

1 - 1

 11 under certain symmetry properties of Ω, or dimension condition, sharper lower bounds are known [14, Section 8.3] for a polynomial p ∈ R[x] ≤1 ∖ R. According to Theorem 3.5 and Theorem 3.6 H B 1 is invertible and M B p = H B H B p is the matrix of the multiplication by p in K[x] I Λ expressed in the basis B. By Theorem 3.2, the eigenvalues of M B p are the values of p at the nodes {ξ 1 , . . . , ξ r } of Λ and W B ξ M B p = DW B ξ where D = diag (p(ξ 1 ), . . . , p(ξ r ))

  we can determine λ and hence the whole vector b 1 (ξ) . . . b r (ξ) T . By computing the r eigenvectors we can thus determine the Vandermonde matrix W B ξ .

Figure 1 :

 1 Figure 1: The regular hexagon H and the cubature H 2 ∶ 5 -1

3 H( 4 ) 1 and 4 . 50 h

 341450 parameters. The 6 × 6 leading principal submatrix has entries in K = Q[ √ 3] and is positive definite. 3. We call Algorithm 4.6 with the matrix A = 160 √ c ′ = 6, r = 7, c ′′ = 10 and c = 15. There are then 4 1 = 4 triplets [P, Z, J], where J = {1, 2, 3, 4, 5, 6, k} with k ∈ {7, 8, 9, 10}. For the two triplets [P, Z, J] where J = {1, 2, 3, 4, 5, 6, k} with k ∈ {9, 10}, the sets of equalities Z are inconsistent. These triplets thus bring no cubature. For the triplet [P, Z, J] where J = {1, 2, 3, 4, 5, 6, 8}, the set of equalities Z has a single solution h 24 = 49 50 , h 42 = 147 50 , h 15 = 0, h 60 = 441 50 , h 33 = 0, h 51 = 0, h 06 = 539 07 = 0, h 16 = 0, h 25 = 0, h 34 = 0, h 43 = 0, h 52 = 0, h 61 = 0, h 70 = 0, h 08 = 14749 2500 , h 17 = 0, h 26 = 343 2500 , h 35 = 0, h 44 = 1029 2500 , h 53 = 0, h 62 = 3087 2500 , h 71 = 0, h 80 = 9261 2500 while the set P has a single non constant polynomial, which is a positive multiple of 50h 42 -49. This triplet corresponds thus to a single cubature.

,

  and for p = x 2 + y 2

√ 14 5 . 14 5 cos k π 3 , √ 14 5 sin k π 3 0

 514143 We can thus conclude that the nodes are :{(0, 0)} ∪ √ ≤ k ≤ 5 .

For all 0

 0 < θ < 2π, a basis of R[x] I Λ is given by C = {d11 , d 12 ; d 31 ; d 51 , d 61 ; d 71 , d 81 }

2 For p = x 2 + y 2

 222

√ 14 25 ,

 25 independently of the parameter θ. The set of nodes is thus

3 1

 3 matrix M C p .Consider now the second generating invariant q = d 15 = y 2 3 x 2y 2 2 . Observe that q cos(α) 2 4 cos(α) 2 -1 2

2 5 4 ( 1 +

 241 cos (θ)) 4⋅7 3 5 6 (1 + cos (θ))

47 3 5 6 ( 1 +≤ k ≤ 5 where α is such that 1 -= 1 2 ( 1 +

 615121 cos (θ)).For the cubature determined by the parameter θ, the nodes are thus cos(α)2 4 cos(α) 2 -1 2 cos (θ)) .

  that we actually have to solve the linear system 1

  2 satisfied by the orthonormal polynomials πr (see e.g. [21, Corollary 1.3.10],[2, Chapter 1]):

  51 h 42 h 33 h 70 h 61 h 52 h 43 h 34 42 h 33 h 24 h 61 h 52 h 43 h 34 h 25 33 h 24 h 15 h 52 h 43 h 34 h 25 h 16 24 h 15 h 06 h 43 h 34 h 25 h 16 h 07 21 0 0 h 60 h 51 h 42 h 70 h 61 h 52 h 43 h 80 h 71 h 62 h 53 h 44 0 0 0 h 51 h 42 h 33 h 61 h 52 h 43 h 34 h 71 h 62 h 53 h 44 h 35 7 0 0 h 42 h 33 h 24 h 52 h 43 h 34 h 25 h 62 h 53 h 44 h 35 h 26 0 0 0 h 33 h 24 h 15 h 43 h 34 h 25 h 16 h 53 h 44 h 35 h 26 h 17 21 0 0 h 24 h 15 h 06 h 34 h 25 h 16 h 07 h 44 h 35 h 26 h 17 h 08

	240	0	0	50	0	50	0	0	0	0	21	0	7	0	21
	0	50	0	0	0	0	21	0	7	0	0	0	0	0	0
	0	0	50	0	0	0	0	7	0	21	0	0	0	0	0
	50	0	0	21	0	7	0	0	0	0	h 60 h 51 h 42 h 33 h 24
	0	0	0	0	7	0	0	0	0	0	h 51 h 42 h 33 h 24 h 15
	50	0	0	7	0	21	0	0	0	0	h 42 h 33 h 24 h 15 h 06
	0 h 60 h 0 21 0 0 0 0 0 7 0 0 0 h 51 h 0 7 0 0 0 0 h 42 h 0 0 21 0 0 0 h 33 h

  [START_REF] Baron De Prony | Essai expérimental et analytique sur les lois de la dilatabilité des fluides élastique et sur celles de la force expansive de la vapeur de l'eau et de la vapeur de l'alkool, à différentes températures[END_REF] , h 33 = -h 15 , h 51 = h 15 , h 07 = 0, h 16 = 0, h 25 = 0, h 34 = 0, h 43 = 0, h 52 = 0, h 61 = 0, h 70 = 0, All the unknowns can be expressed explicitly in terms of either h 60 or h 15 while

				25	+ h 60 , h 42 =	294 25	-h 60 , h 06 = -h h 44 = 98 5 1029 2500 ,
	h 08 =	7889 500	-	28 25	h 60 , h 26 = -	2401 500	+	14 25	h 60 , h 62 =	3087 500	-	14 25	h 60 , h 80 = -	3087 500	+	28 25	h 60 ,
				h 35 = -	7 25	h 15 , h 53 = -	7 25	h 15 , h 71 =	21 25	h 15 , h 17 =	21 25	h 15 .
								50 2 h 60 -	49 5	2	+ 50 2 h 15	2 = 49 2 .

  With this exact computation, we also get the exact values of the weights (a 1 , a 2 , a 3 , a 4 , a 5 , a 6 , a 7 ) =

																				5	, 0 ,
	ξ 4 =	√ 14 10	,	√ 10 42	, ξ 5 = -	√ 10 14	,	√ 10 42	, ξ 6 = -	√ 10 14	, -	√ 10 42	, ξ 7 =	√ 10 14	, -	√ 10 42	.
					43 112 √ 3	,	125 672 √	3	,	125 672 √	3	,	125 672 √	3	,	125 672 √	3	,	125 672 √	3	,	125 672 √ 3	.

  11 , d 12 , d 13 } ∪ {d 31 , d 41 } ∪ {d 51 , d 61 , d 52 , d 62 } ∪ {d 71 , d 81 , d 72 , d 73 , d 82 , d 83 } .

  To determine the existence conditions of a cubature we can thus apply Algorithm 4.6 separately to H(1) , H(2) , H(3) and H(4) . We thus deal with matrices with sizes 3 × 3, 2 × 2, 4 × 4, 6 × 6 instead of a single matrix of size 15 × 15.

	20 0 h 14 h 15 -h 21 h 17 h 17 -h 22 0 h 14 h 15 0 7 20 0 -h 21 h 14 h 15 h 21 0 h 16 h 17 0 h 22 0 0 h 14 0 -h 22 h 16 h 17 h 21 h 15 h 22 0 h 17 h 17	⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

  49 250 , h 15 = 49 250 , h 16 = 343 3125 , h 17 = 343 3125 , h 21 = 0, h 22 = 0. Since h 21 = 0, h 22 = 0, we know that it determines a D 6 -invariant cubature. A basis of R[x] I Λ is given by B = {d 11 , d 12 ; d 41 ; d 51 , d 61 ; d 71 , d 81 } The distribution of the elements of B

September 2018

To appear in Studies in Applied Mathematics, OPSFA special issue.

Extension and inside-ness

When looking at the cubatures for Ω ∶ R[x] → R given as the integration according to a measure µ, i.e. Ω(p) = ∫ R n pdµ, it might be necessary to ensure that the nodes lie on the support K of this measure, i.e. that the cubature is inside. When K is a semi-algebraic set it is possible to test this property of a cubature before actually computing its nodes. The key property is Corollary 3.11. Assume

and consider κ ∈ N be such that 2κ ≥ deg g k for all k = 1, . . . , s.

In previous section we have characterized a cubature Λ of degree d for Ω with r nodes thanks to its restriction Λ (δ) to R[x] ≤2δ , for some δ ∈ N. Λ (δ) being a flat extension of Λ (δ-1) , it determines uniquely Λ ∶ R[x] → R. There are several reasons for wanting to know an intermediate flat extension Λ (δ+κ) of Λ (δ) , as for instance testing if the cubature is inside.

To directly obtain this extension Algorithm 5.3 can be slightly modified as follow: In Step 2. we introduce

In Step 3. we then call Algorithm 4.6 with A ∶= H

The triplets produced by this modification are in one to one correspondence with the triplets produced by the original algorithm, only their component Z has more elements. As guaranteed by Theorem 3.12 one can combine the elements of Z so as to express h r 2δ +1 , . . . , h r 2δ+κ explicitly in terms of h r d +1 , . . . , h r 2δ .

Computation of the weights and the coordinates of the nodes

Assume that a cubature Λ of degree d for Ω with r nodes has been identified by a solution ( ̵ h r d +1 , . . . , ̵ h r 2δ ) of one triplet [P, Z, J] output by Algorithm 5.3. This means that Λ is the unique flat extension of Λ

and that B = {m

Λ (δ) and hence Λ are positive semi-definite. Their Hankel operators are of rank r. By Theorem 3.9 Λ = ∑ r i=1 a i e ξi for some a i > 0 and ξ i ∈ R n . The nodes {ξ 1 , . . . , ξ r } are the roots of the ideal I Λ , which is generated by ker Ĥ(δ) .

Following Section 3 we can compute the nodes and the weights of the cubature Λ with linear algebra operations. We describe the process with the light assumption that d ≥ 2 and that Ω(p 2 ) > 0 for all p ∈ R[x] ≤1 ∖ {0} as it avoids some cumbersome case distinctions. With this assumption, B = {b 1 , . . . , b r } contains M (1) , the basis of R[x] ≤1 .

Computing the nodes as eigenvectors

First we form the symmetric matrices

All the real solutions of this equation are on a circle and thus can be parameterized with a single parameter θ as Taking into account the positivity constraint 50h 60 -441 = 49(1 + cos θ) > 0, all the solutions of the triplet are parametrized explicitly by θ, for θ ∈]π, π[. Note though that the solution of Z for θ = π is the solution of the previous triplet.

Inside-ness

Proposition 3.11 gives us a method to show that the cubatures characterized above are inside for every value of θ ∈]π, π].

The regular hexagon H is defined as a semi-algebraic set by

Taking the basis M (3) = {1, x 1 , x 2 , x 2 1 , x 1 x 2 , x 2 2 , x 3 1 , x 2 1 x 2 , x 1 x 2 2 , x 3 2 } of R[x] ≤3 , we construct the localizing matrices H