Angelos Mantzaflaris
email: angelos.mantzaflaris@inria.fr

Bernard Mourrain
email: bernard.mourrain@inria.fr

Agnes Szanto
email: aszanto@ncsu.edu

Punctual Hilbert Schemes and Certified Approximate Singularities

Keywords: certification, singularity, multiplicity structure, Newton's method, inverse system, multiplication matrix

In this paper we provide a new method to certify that a nearby polynomial system has a singular isolated root and we compute its multiplicity structure. More precisely, given a polynomial system f = (f 1 , . . . , f N ) ∈ C[x 1 , . . . , xn] N , we present a Newton iteration on an extended deflated system that locally converges, under regularity conditions, to a small deformation of f such that this deformed system has an exact singular root. The iteration simultaneously converges to the coordinates of the singular root and the coefficients of the so-called inverse system that describes the multiplicity structure at the root. We use α-theory test to certify the quadratic convergence, and to give bounds on the size of the deformation and on the approximation error. The approach relies on an analysis of the punctual Hilbert scheme, for which we provide a new description. We show in particular that some of its strata can be rationally parametrized and exploit these parametrizations in the certification. We show in numerical experimentation how the approximate inverse system can be computed as a starting point of the Newton iterations and the fast numerical convergence to the singular root with its multiplicity structure, certified by our criteria.

Introduction

Local numerical methods such as Newton iterations have proved their efficiency to approximate and certify the existence of simple roots. However for multiple roots they dramatically fail to provide fast numerical convergence and certification. The motivation for this work is to find a method with fast convergence to an exact singular point and its multiplicity structure for a small perturbation of the input polynomials, and to give numerical tests that can certify it. The knowledge of the multiplicity structure together with a high precision numerical approximation of a singular solution can be valuable information in many problems.

In [START_REF] Mourrain | Isolated points, duality and residues[END_REF] a method called later integration method is devised to compute the so-called inverse system or multiplicity structure at a multiple root. It is used in [START_REF] Mantzaflaris | Deflation and certified isolation of singular zeros of polynomial systems[END_REF] to compute an approximation of the inverse system, given an approximation of that root and to obtain a perturbed system that satisfies the duality property. However, this method did not give a way to improve the accuracy of the initial approximation of the root and the corresponding inverse system. In [START_REF] Hauenstein | On deflation and multiplicity structure[END_REF] a new one-step deflation method is presented that gives an overdetermined polynomial system in the coordinates of the roots and the corresponding inverse system, serving as a starting point for the present paper. However, for certification, [START_REF] Hauenstein | On deflation and multiplicity structure[END_REF] refers to the symbolic-numeric method in [START_REF] Ayyildiz Akoglu | Certifying solutions to overdetermined and singular polynomial systems over Q[END_REF] that only works if the input system is given exactly with rational coefficients and have a multiple root with the prescribed multiplicity structure.

In the present paper we give a solution for the following problem:

Problem 1.1. Given a polynomial system f = (f 1 , . . . , f N ) ∈ C[x] N and a point ξ ∈ C n , deduce an iterative method that converges quadratically to the triple (ξ * , µ * , ǫ * ) such that ξ * ∈ C n , µ * defines the coefficients of a basis Λ * = {Λ * 1 , . . . , Λ * r } ⊂ C[d ξ * ] dual to the set B ξ * = {(xξ * ) β1 , . . . , (xξ * ) βr } ⊂ C[x] and ǫ * defines a perturbed polynomial system f ǫ * := f + ǫ * B ξ * with the property that ξ * is an exact multiple root of f ǫ * with inverse system Λ * . Furthermore, certify this property and give an upper bound on the size of the perturbation ǫ * .

The difficulty in solving Problem 1.1 is that known polynomial systems defining the coordinates of the roots and the inverse system are overdetermined, and we need a square subsystem of it in the Newton iterations to guarantee the existence of a root together with the quadratic convergence. Thus, roots of this square subsystem may not be exact roots of the complete polynomial system, and we cannot certify numerically that they are approximations of a root of the complete system. This is the reason why we introduce the variables ǫ that allow perturbation of the input system. One of the goals of the present paper is to understand what kind of perturbations are needed and to bound their magnitude.

Certifying the correctness of the multiplicity structure that the numerical iterations converge to poses a more significant challenge: the set of parameter values describing an affine point with multiplicity r forms a projective variety called the punctual Hilbert scheme. The goal is to certify that we converge to a point on this variety. We study an affine subset of the punctual Hilbert scheme and give a new description using multilinear quadratic equations that have a triangular structure. These equations appear in our deflated polynomial system, have integer coefficients, and have to be satisfied exactly without perturbation, otherwise the solution does not define a proper inverse system, closed under derivation. Fortunately, the structure allowed us to define a rational parametrization of a strata of the punctual Hilbert scheme, called the regular strata. In turn, this rational parametrization allows certification when converging to a point on this regular strata.

Our method comprises three parts: first, we apply the Integration Method (Algorithm 1) with input f and ξ to compute an approximation of the multiplicity structure, second, an analysis and certification part (see Section 6 and Algorithm 2), and third, a numerical iteration part converging to the exact multiple root with its multiplicity structure for an explicit perturbation of the input system (see Section 5).

Related Work. There are many works in the literature studying the certification of isolated singular roots of polynomial systems. One approach is to give separation bounds for isolated roots, i.e. a bound that guarantees that there is exactly one root within a neighborhood of a given point. Worst case separation bounds for square polynomial systems with support in given polytopes and rational coefficients are presented in [START_REF] Emiris | The DMM bound: multivariate (aggregate) separation bounds[END_REF]. In the presence of singular roots, turned into root clusters after perturbations, these separation bounds separate the clusters from each other and bound the cluster size. [START_REF] Yakoubsohn | Finding a cluster of zeros of univariate polynomials[END_REF][START_REF] Yakoubsohn | Simultaneous computation of all the zero-clusters of a univariate polynomial[END_REF][START_REF] Giusti | On location and approximation of clusters of zeros of analytic functions[END_REF] give separation bounds and numerical algorithms to compute clusters of zeroes of univariate polynomials. [START_REF] Dedieu | On simple double zeros and badly conditioned zeros of analytic functions of n variables[END_REF] extends α-theory and gives separation bounds for simple double zeroes of polynomial systems, [START_REF] Giusti | On location and approximation of clusters of zeros: Case of embedding dimension one[END_REF] extend these results to zeroes of embedding dimension one.

Another approach, called deflation, comprises of transforming the singular root into a regular root of a new system and to apply certification techniques on the new system. [START_REF] Kanzawa | Approximate singular solutions of nonlinear equations and a numerical method of proving their existence[END_REF] uses a square deflated system to prove the existence of singular solutions. [START_REF] Leykin | Newton's method with deflation for isolated singularities of polynomial systems[END_REF] devises a deflation technique that adds new variables to the systems for isolated singular roots that accelerates Newton's method and [START_REF] Leykin | Higher-order deflation for polynomial systems with isolated singular solutions[END_REF] modifies this to compute the multiplicity structure. [START_REF] Rump | Verified error bounds for multiple roots of systems of nonlinear equations[END_REF] computes error bounds that guarantee the existence of a simple double root within that error bound from the input, [START_REF] Li | Verified error bounds for isolated singular solutions of polynomial systems: case of breadth one[END_REF][START_REF] Li | Verified error bounds for isolated singular solutions of polynomial systems[END_REF] generalizes [START_REF] Rump | Verified error bounds for multiple roots of systems of nonlinear equations[END_REF] to the breadth one case and give an algorithm to compute such error bound. [START_REF] Li | Verified error bounds for singular solutions of nonlinear systems[END_REF] gives verified error bounds for isolated and some nonisolated singular roots using higher order deflations. [START_REF] Dayton | Computing the multiplicity structure in solving polynomial systems[END_REF][START_REF] Wu | Computing the multiplicity structure from geometric involutive form[END_REF][START_REF] Zeng | The closedness subspace method for computing the multiplicity structure of a polynomial system[END_REF][START_REF] Wu | Determining singular solutions of polynomial systems via symbolic-numeric reduction to geometric involutive form[END_REF][START_REF] Dayton | Multiple zeros of nonlinear systems[END_REF][START_REF] Hao | Algorithm 931: an algorithm and software for computing multiplicity structures at zeros of nonlinear systems[END_REF] give deflation techniques based on numerical linear algebra on the Macaulay matrices that compute the coefficients of the inverse system, with improvements using the closedness property of the dual space. [START_REF] Giusti | Multiplicity hunting and approximating multiple roots of polynomial systems[END_REF][START_REF] Giusti | Approximation numérique de racines isolées multiples de systèmes analytiques[END_REF] give a new deflation method that does not introduce new variables and extends α-theory to general isolated multiple roots for the certification to a simple root of a subsystem of the overdetermined deflated system. In [START_REF] Hauenstein | On deflation and multiplicity structure[END_REF] a new deflated system is presented, its simple roots correspond to the isolated singular points with their multiplicity structure. A somewhat different approach is given in [START_REF] Ayyildiz Akoglu | Certifying solutions to overdetermined and singular polynomial systems over Q[END_REF], where they use a symbolic-numeric certification techniques that certify that polynomial systems with rational coefficients have exact isolated singular roots. More recently, [START_REF] Lee | On isolation of singular zeros of multivariate analytic systems[END_REF] design a square Newton iteration and provide separation bounds for roots when the deflation method of [START_REF] Leykin | Newton's method with deflation for isolated singularities of polynomial systems[END_REF] terminates in one iteration, and give bounds for the size of the clusters.

The certification approach that we propose is based on an algebraic analysis of some strata of the punctual Hilbert scheme. Some of its geometric properties have been investigated long time ago, for instance in [START_REF] Brianc ¸on | Description de Hilb n C{x, y}[END_REF][START_REF] Iarrobino | Punctual Hilbert Schemes[END_REF][START_REF] Brianc ¸on | Dimension of the punctual Hilbert scheme[END_REF] or more recently in the plane [START_REF] Bejleri | The tangent space of the punctual Hilbert scheme[END_REF]. However, as far as we know, the effective description that we use and the rational parametrization of the regular strata that we compute have not been developed previously.

Preliminaries

Let f := (f 1 , . . . , f N ) ∈ C[x] N with x = (x 1 , . . . , x n ). Let ξ = (ξ 1 , . . . , ξ n ) ∈ C n be
an isolated multiple root of f . Let I = f 1 , . . . , f N , m ξ be the maximal ideal at ξ and Q be the primary component of I at ξ so that √ Q = m ξ . The shifted monomials at ξ will be denoted for α = (α 1 , . . . , α n ) ∈ N n by

x α ξ := (x 1 -ξ 1 ) α1 • • • (x 1 -ξ n ) αn .

Consider the ring of power series C[[d

ξ ]] := C[[d 1,ξ , . . . , d n,ξ ]] and we denote d β ξ := d β1 1,ξ • • • d βn n,ξ , with β = (β 1 , . . . , β n ) ∈ N n . We identify C[[d ξ ]
] with the dual space C[x] * by considering the action of d β ξ on polynomials as derivations and evaluations at ξ, defined as (1)

d β ξ (p) := ∂ β (p) ξ = ∂ |β| p ∂x β1 1 • • •∂x βn n (ξ) for p ∈ C[x].
Hereafter, we reserve the notation d and d i for the dual variables while ∂ and ∂ xi for derivation. We indicate the evaluation at ξ ∈ C n by writing d i,ξ and d ξ , and for ξ = 0 it will be denoted by d. The derivation with respect to the variable

d i,ξ in C[[d ξ ]] is denoted ∂ d i,ξ (i = 1, . . . , n). Observe that 1 β! d β ξ ((x -ξ) α ) = 1 if α = β, 0 otherwise, where β! = β 1 ! • • • β n !. For p ∈ C[x] and Λ ∈ C[[d ξ ]] = C[x] * , let p•Λ : q → Λ(p q). We check that p = (x i -ξ i ) acts as a derivation on C[[d ξ ]]: (x i -ξ i ) • d β ξ = ∂ d i,ξ (d β ξ ) = β i d β-ei ξ
. Throughout the paper we use the notation e 1 , . . . , e n for the standard basis of C n or for a canonical basis of any vector space V of dimension n. We will also use

integrals of polynomials in C[[d ξ ]] as follows: for Λ ∈ C[[d ξ ]] and k = 1, . . . , n, ∫ k Λ denotes the polynomial Λ * ∈ C[[d ξ ]] such that ∂ d k,ξ (Λ * ) = Λ
and Λ * has no constant term. We introduce the following shorthand notation For t ∈ N, let D t be the elements of D of order ≤ t. As D is of dimension r, there exists a smallest t ≥ 0 s.t. D t+1 = D t . Let us call this smallest t, the nil-index of D and denote it by δ ξ (I), or simply by δ. As D is stable by the derivations ∂ d i,ξ , we easily check that for t ≥ δ ξ (I), D t = D and that δ ξ (I) is the maximal degree of elements of D.

∫ k Λ := ∫ k Λ(d 1,ξ , . . . , d k,ξ , 0, . . . , 0). (2) For an ideal I ⊂ C[x], let I ⊥ = {Λ ∈ C[[d ξ ]] | ∀p ∈ I, Λ(p) = 0}. The vector space I ⊥ is naturally identified with the dual space of C[x]/I. We check that I ⊥ is a vector subspace of C[[d ξ ]] which is closed under the derivations ∂ d i,ξ for i = 1, . . . , n. Lemma 2.1. If Q is a m ξ -primary isolated component of I, then Q ⊥ = I ⊥ ∩ C[d ξ ].
Let B = {x β1 ξ , . . . , x βr ξ } be a basis of C[x]/Q. We can identify the elements of C[x]/Q with the elements of the vector space span C (B). We define the normal form N (p) of a polynomial p in C[x] as the unique element b of span C (B) such that pb ∈ Q. Hereafter, we are going to identify the elements of C[x]/Q with their normal form in span C (B). For α ∈ N n , we will write the normal form of x α ξ as

(3) N (x α ξ ) = r i=1 µ βi,α x βi ξ .
2.1. The multiplicity structure. We start this subsection by recalling the definition of graded primal-dual pairs of bases for the space C[x]/Q and its dual. The following lemma defines the same dual space as in e.g. [START_REF] Dayton | Computing the multiplicity structure in solving polynomial systems[END_REF][START_REF] Dayton | Multiple zeros of nonlinear systems[END_REF][START_REF] Li | Verified error bounds for isolated singular solutions of polynomial systems[END_REF]], but we emphasize on a primal-dual basis pair to obtain a concrete isomorphism between the coordinate ring and the dual space. We can assume that β 1 = 0 and that the ordering of the elements in B by increasing degree. Define the set of exponents in B as

E := {β 1 , . . . , β r } ⊂ N n . (2) The unique dual basis Λ = {Λ 1 , Λ 2 , . . ., Λ r } of D ⊂ C[d ξ ] dual to B has the form Λ i = 1 βi! d βi ξ + |α|≤|β i | α ∈E µ βi,α 1
β! d α ξ . (3) We have 0 = ord(Λ 1 ) ≤ • • • ≤ ord(Λ r )
, and for all 0 ≤ t ≤ δ we have D t = span {Λ j : ord(Λ j ) ≤ t} , where D t denotes the elements of D of order ≤ t, as above.

A graded primal-dual basis pair (B, Λ) of D as described in Lemma 2.2 can be obtained from any basis Λ of D by first choosing pivot elements that are the leading monomials with respect to a graded monomial ordering on C[d], these leading monomials define B, then transforming the coefficient matrix of Λ into row echelon form using the pivot leading coefficients, defining Λ.

A monomial set B is called a graded primal basis of f at ξ if there exists Λ ⊂ C[d ξ ] such that (B, Λ) is a graded primal-dual basis pair and Λ is complete for f at ξ.

Next we describe the so-called integration method introduced in [START_REF] Mourrain | Isolated points, duality and residues[END_REF][START_REF] Mantzaflaris | Deflation and certified isolation of singular zeros of polynomial systems[END_REF] that computes a graded pair of primal-dual bases as in Lemma 2.2 if the root ξ is given. The integration method performs the computation of a basis order by order. We need the following proposition, a new version of [27, Theorem 4.2]:

Proposition 2.3. Let Λ 1 , . . . , Λ s ∈ C[d ξ ]
and assume that ord(Λ i ) ≤ t for some t ∈ N. Suppose that the subspace

D := span(Λ 1 , . . . , Λ s ) ⊂ C[d ξ ] is closed under derivation. Then ∆ ∈ C[d ξ ] with no constant term satisfies ∂ d k (∆) ∈ D for all k = 1, . . . , n if and only if ∆ is of the form ∆ = s i=1 n k=1 ν k i ∫ k Λ i (5) for some ν k i ∈ C satisfying s i=1 ν k i ∂ d l (Λ i ) -ν l i ∂ d k (Λ i ) = 0 for 1 ≤ k < l ≤ n. (6)
Furthermore, ( 5) and [START_REF] Dayton | Multiple zeros of nonlinear systems[END_REF] implies that

∂ d k (∆) = s i=1 ν k i Λ i for k = 1, . . . , n. (7)
Proof. Suppose Λ ∈ C[d] with no constant term satisfies ∂ d k (Λ) ∈ D for all k = 1, . . . , n. To prove [START_REF] Brianc ¸on | Dimension of the punctual Hilbert scheme[END_REF], we can proceed exactly as in the proof of [START_REF] Mourrain | Isolated points, duality and residues[END_REF]Theorem 4.2]: we write ∆ uniquely as

∆ = ∆ 1 (d 1 , . . . , d n ) + ∆ 2 (d 2 , . . . , d n ) + • • • , +∆ n (d n ) with ∆ i ∈ C[d i , . . . , d n ] \ C[d i+1 , . . . , d n ]. Then ∫ i ∂ di ∆ i = ∆ i .
Then we prove that by induction on k that if

σ k := ∆ 1 + • • • + ∆ k then ∆ k = s j=1 ν k j ∫ k Λ j -(σ k-1 -σ k-1 | d k =0 )
and

σ k = ∆ k + σ k-1 = s j=1 ν k j ∫ k Λ j + σ k-1 | d k =0 = s j=1 ν k j ∫ k Λ j + s j=1 ν k-1 j ∫ k Λ j | d k =0 + • • • + s j=1 ν 1 j ∫ k Λ j | d k =0,•••d2=0 .
Conversely, suppose that Λ ∈ C[d] with no constant term is of the form (5) satisfying [START_REF] Dayton | Multiple zeros of nonlinear systems[END_REF]. Define ∆1 = σ1 :=

s j=1 ν 1 j ∫ 1 Λ j and for k = 2, . . . n define ∆k := s j=1 ν k j ∫ k Λ j -(σ k-1 -σ k-1 | d k =0 )
and σk := ∆1 + • • • + ∆k . Then in the proof of [START_REF] Mourrain | Isolated points, duality and residues[END_REF]Theorem 4.2]

it is shown that ∆k ∈ C[d k , . . . , d n ] \ C[d k+1 , . . . , d n ] and σk = s j=1 ν k j ∫ k Λ j + s j=1 ν k-1 j ∫ k Λ j | d k =0 + • • • + s j=1 ν 1 j ∫ k Λ j | d k =0,•••d2=0 so we get that ∂ d k (Λ) = ∂ d k (σ k ) = s j=1 ν k j Λ j ∈ D t as claimed. Let Q be a m ξ -primary ideal. Proposition 2.3 implies that if Λ = {Λ 1 , . . . , Λ r } ⊂ C[d ξ ] with Λ 1 = 1 ξ is a basis of Q ⊥ , dual to the basis B = {x β1 ξ , . . . , x βr ξ } ⊂ C[x] of C[x]/Q with ord(Λ i ) = |β i |, then there exist ν k i,j ∈ C such that ∂ d k (Λ i ) = |βj|<|βi| ν k i,j Λ j .
Therefore, the matrix M k of the multiplication map

M k by x k -ξ k in the basis B of C[x]/Q is M k = [ν k j,i ] T 1≤i,j≤r = [µ βi,βj +e k ] 1≤i,j≤r using the notation (3) and the convention that ν k i,j = µ βi,βj +e k = 0 if |β i | ≥ |β j |. Consequently, ν k i,j = µ βi,βj+e k i, j, = 1, . . . , r, k = 1, . . . ,

n, and we have

Λ i = |βj|<|βi| n k=1 µ βi,βj+e k ∫ k Λ j
where µ βi,βj+e k is the coefficient of x βi ξ in the normal form of x βj +e k ξ in the basis B of C[x]/Q. Next we give a result that allows to simplify the linear systems involved in the integration method. We first need a definition:

Definition 2.4. Let E ⊂ N n be a set of exponents. We say that E is closed under division if β = (β 1 , . . . , β n ) ∈ E implies that β -e k ∈ E
as long as β k > 0 for all k = 1, . . . , n. We also call the corresponding primal basis B = {x β1 ξ , . . . , x βr ξ } closed under division. The following lemma provides a simple characterization of dual bases of inverse systems closed under derivation, that we will use in the integration algorithm.

Λ i = |βj |<|βi| n k=1 µ βi,βj+e k ∫ k Λ j . ( 8
)
Then D = span{Λ 1 , . . . , Λ r } is closed under derivation iff for all i, s = 1, . . . , r, |β s | < |β i | and k = l ∈ {1, . . . , n} we have j:|βs|<|βj|<|βi| µ βi,βj +e k µ βj,βs+e l -µ βi,βj+e l µ βj,βs+e k = 0. (9)
Furthermore, (B, Λ) is a graded primal-dual basis pair iff they satisfy (9) and

µ βi,βj+e k = 1 for β i = β j + e k 0 for β j + e k ∈ E, β i = β j + e k , (10)
Proof of Lemma 2.5. Assume Λ = {Λ 1 , . . . , Λ r } is linearly independent and D = span(Λ) is closed under derivation. For t ∈ {0, . . . , δ} denote by {Λ 1 , . . . , Λ rt } = Λ ∩ C[d ξ ] t and D t = span(Λ 1 , . . . , Λ rt ). Then by Proposition 2.3, Λ satisfy equations ( 7) for t = 0, . . . , δ and for j = 1, . . . , r, k = 1, . . . , n, we have

∂ d k (Λ j ) = |βs|<|βj| µ βj,βs+e k Λ s . Substituting this to (6) we get for i = 1, . . . , r |βj|<|βi| µ βi,βj+e k |βs|<|βj| µ βj,βs+e l Λ s -µ βi,βj+e l |βs|<|βj| µ βj ,βs+e k Λ s = 0. ( 11
)
Then using linear independence and collecting the coefficients of Λ s we get [START_REF] Doubilet | On the Foundations of Combinatorial Theory[END_REF]. Conversely, assume that ( 9) is satisfied. Then ( 11) is also satisfied. We use induction on t to prove that D t is closed under derivation. For t = 0 there is nothing to prove. Assume D t-1 is closed under derivation. Then by Proposition 2.3 if

|β j | < t then ∂ d k (Λ j ) = |βs|<|βj| µ βj,βs+e k Λ s for k = 1, . . . , n. Thus for |β i | = t, (11) implies that |βj|<|βi| µ βi,βj+e k ∂ d l (Λ j ) -µ βi,βj +e l ∂ d k (Λ j ) = 0.
Again, by Proposition 2.3 we get that D t is closed under derivation. Next, assume first that (B, Λ) is a graded primal-dual basis pair. This means that for i = 1, . . . , r and for l such that

|β l | ≤ |β i | δ i,l = Λ i x β l ξ = n k=1 |βj |<|βi| µ βi,βj +e k ∫ k Λ j x β l ξ = n k=1 |βj |<|βi| µ βi,βj +e k coeff( d β l β l ! , ∫ k Λ j )
Fix k to be the index of the last non-zero entry of β l . For all other k's d β l becomes zero when we substitute

0 into d k+1 , . . . , d n in ∫ k Λ j . Thus, Λ i x β l ξ = |βj|<|βi| µ βi,βj+e k coeff( d β l β l ! , ∫ k Λ j ) = |βj|<|βi| µ βi,βj+e k coeff( d β l -e k (β l -e k )! , Λ j ).
Since E is closed under division, β l -e k = β m ∈ E for some m < l. By duality, we have that coeff(

d βm (βm)! , Λ j ) = δ m,j , so Λ i x β l ξ = µ βi,βm+e k = µ βi,β l . To satisfy Λ i x β l ξ = δ i,l we must have µ βi,βm+e k = 1 if β i = β m + e k 0 if β m + e k = β l ∈ E but i = l.

Conversely, by induction on

t = |β i | we have that deg(Λ i ) ≤ |β i |. Then Λ i x β l ξ = 0 when |β l | > |β i |.
For

|β l | ≤ |β i |, Equations (10) imply that the coefficient of d β l β l ! in Λ i is 0 if i = l and 1 if i = l. Therefore (B, Λ) is a graded primal-dual basis pair.
To compute the inverse system D of f at a point ξ, we will consider the additional systems of equations in ξ and µ = {µ βi,α }:

(12) Λ i (f j ) = 0 for 1 ≤ i ≤ r, 1 ≤ j ≤ N.
Throughout the paper we use the following notation:

Notation 2.6. Let f 1 , . . . , f N ∈ C[x], ξ ∈ C n and fix t ∈ N. Let B t-1 = {x β1 ξ , . . . , x βr t-1 ξ } ⊂ C[x ξ ] t-1 be closed under division and Λ t-1 = {Λ 1 , . . . , Λ rt-1 } ⊂ C[d ξ ] t-1 dual to B t-1 with ∂ d k (Λ j ) = |βs|<|βj| µ βj ,βs+e k Λ s j = 1, . . . , r t-1 , k = 1, . . . , n.
Consider the following homogeneous linear system of equations in the variables {ν k j : j = 1, . . . , r t-1 , k = 1, . . . , n}:

j:|βs|<|β j |<t ν k j µ β j ,βs+e l -ν l j µ β j ,βs+e k = 0, 1 ≤ k < l ≤ n (13) ν k j = 0 if βj + e k = β l for some 1 ≤ l ≤ rt-1 (14) r t-1 j=1 n k=1 ν k j ∫ k Λj (f l ) = 0 l = 1, . . . , N. ( 15
)
We will denote by H t the coefficient matrix of the equations in ( 13) and ( 14) and by K t the coefficient matrix of the equations in ( 13)- [START_REF] Hao | Algorithm 931: an algorithm and software for computing multiplicity structures at zeros of nonlinear systems[END_REF]. By Proposition 2.3 and Lemma 2.5, if

K t ν = 0 where ν = [ν k j : j = 1, . . . , s, k = 1, . . . , n], then Λ = s j=1 n k=1 ν k j ∫ k Λ j ∈ (f ) ⊥ ∩ C[d ξ ] t = D t .
The main loop of the integration method described in Algorithm 1 consists of computing the new basis elements in D t and the new basis monomials in B t of degree t from the primal-dual basis pair (B t-1 , Λ t-1 ) in degree t -1.

Algorithm 1 produces incrementally a basis of D, similarly to Macaulay's method. The algorithmic advantage is the smaller matrix size in

O(r n 2 + N ) instead of N n+δ-1 δ
, where δ is the maximal degree (depth) in the dual, cf. [START_REF] Mantzaflaris | Deflation and certified isolation of singular zeros of polynomial systems[END_REF][START_REF] Hauenstein | On deflation and multiplicity structure[END_REF]. [START_REF] Dedieu | On simple double zeros and badly conditioned zeros of analytic functions of n variables[END_REF]. Output: Either "D t = D t-1 " or B t = {x β1 ξ , . . . , x βr t ξ } for some r t > r t-1 closed under division and Λ t = {Λ 1 , . . . , Λ rt } with Λ i of the form (8), satisfying (9), [START_REF] Emiris | The DMM bound: multivariate (aggregate) separation bounds[END_REF] and [START_REF] Giusti | On location and approximation of clusters of zeros: Case of embedding dimension one[END_REF].

Algorithm 1 Integration Method -Iteration t Input: t > 0, f = (f 1 , . . . , f N ) ∈ C[x] N , ξ ∈ C n , B t-1 = {x β1 ξ , . . . , x βr t-1 ξ } ⊂ C[x] closed under division and Λ t-1 = {Λ 1 , . . . , Λ rt-1 } ⊂ C[d ξ ] a basis for D t-1 dual to B t-1 , of the form
(1) Set up the coefficient matrix K t of the homogeneous linear system (13)-( 15) in Notation 2.6 in the variables {ν k j } j=1,...,rt-1, k=1,...,n associated to an element of the form Λ =

rt-1 j=1 n k=1 ν k j ∫ k Λ j . Let h t := dim ker K t . (2) If h t = 0 then return "D t = D t-1 ". If h t > 0 define r t := r t-1 + h t .
Perform a triangulation of H t by row reductions with row permutations and column pivoting so that the non-pivoting columns correspond to exponents β rt-1+1 , . . . ,

β rt with strict divisors in B t-1 . Let B t = B t-1 ∪ {x βr t-1 +1 ξ , . . . , x βr t ξ }. (3) Compute a basis Λ rt-1+1 , . . . , Λ rt ∈ C[d ξ ] of ker K t from
the triangular reduction of H t by setting the coefficients of the non-pivoting columns to 0 or 1. This yields a basis 8) so that Eq. ( 12) are satisfied. Eq. ( 10) are satisfied, since Λ t is dual to B t .

Λ t = Λ t-1 ∪ {Λ rt-1+1 , . . . , Λ rt } dual to B t . The coefficients ν k i,j of Λ i are µ βi,βj +e k in (
The full Integration Method consists of taking Λ 1 := 1 ξ for t = 0, a basis of D 0 and then iterating algorithm Integration Method -Iteration t until we find a value of t when D t = D t-1 . This implies that the order δ = δ ξ (f ) = t -1. This leads to the following definition. Definition 2.7. We say that Λ ⊂ C[d ξ ] is complete for f at ξ if the linear system K t of the equations ( 13)- [START_REF] Hao | Algorithm 931: an algorithm and software for computing multiplicity structures at zeros of nonlinear systems[END_REF]

in degree t = δ + 1 = ord(Λ) + 1 is such that ker K δ+1 = {0}.
Notice that the full Integration Method constructs a graded primal-dual basis pair (B, Λ). The basis Λ ⊂ (f ) ⊥ spans a space stable by derivation and is complete for f , so that we have span

(Λ) = (f ) ⊥ ∩ C[d ξ ] = Q ⊥ where Q is the primary component of (f ) at ξ.
To guarantee that B t is closed under division, one could choose a graded monomial ordering ≺ of C[d ξ ] and compute an auto-reduced basis of ker H t such that the initial terms for ≺ are d βi ξ . The set B t constructed in this way would be closed under division, since D t is stable under derivation. In the approach we use in practice, we choose the column pivot taking into account the numerical values of the coefficients and not according to a monomial ordering and we check a posteriori that the set of exponents is closed under division (See Example 7.1).

The main property that we will use for the certification of multiplicities is given in the next theorem.

Theorem 2.8. If ξ * is an isolated solution of the system f (x) = 0 and B is a graded primal basis at ξ * closed under division, then the system F (ξ, µ) = 0 of all equations (9), ( 10) and ( 12) admits (ξ * , µ * ) as an isolated simple root, where µ * defines the basis Λ * of the inverse system of (f ) at ξ dual to B, due to [START_REF] Dedieu | On simple double zeros and badly conditioned zeros of analytic functions of n variables[END_REF].

Proof. This is a direct consequence of [START_REF] Hauenstein | On deflation and multiplicity structure[END_REF]Theorem 4.11], since the system of equations ( 9)-( 12) is equivalent to the system ( 14) in [START_REF] Hauenstein | On deflation and multiplicity structure[END_REF]Theorem 4.11]. The equations ( 9) express the commutation of the transposed of the parametric operator of multiplication in B, which are the same as the equations of commutation of the operators. By Lemma 2.5, the equations ( 10) are equivalent to the fact that (B, Λ * ) is a graded primal-dual basis pair. Finally, the equations ( 12) are the same as N (f i ) = 0, i = 1, . . . , s where N is the parametric normal form defined in [START_REF] Hauenstein | On deflation and multiplicity structure[END_REF][see Definition 4.7 and following remark]. Therefore the two systems are equivalent. By [START_REF] Hauenstein | On deflation and multiplicity structure[END_REF]Theorem 4.11], they define the simple isolated solution (ξ * , µ * ), where µ * defines the basis Λ * dual to B due to (8).

Punctual Hilbert scheme

The results in Sections 3 and 4 do not depend on the point ξ ∈ C n , so to simplify the notation, we assume in these sections that ξ = 0. Let m = (x 1 , . . . , x n ) be the maximal ideal defining

ξ = 0 ∈ C n . Let C[d] be the space of polynomials in the variables d = (d 1 , . . . , d n ) and C[d] t ⊂ C[d] the subspace of polynomials in d of degree ≤ t.
For a vector space V , let G r (V ) be the projective variety of the r dimensional linear subspaces of V , also known as the Grassmannian of r-spaces of V . The points in G r (V ) are the projective points of

P(∧ r V ) of the form v = v 1 ∧ • • • ∧ v r for v i ∈ V .
α i ∈ N n , |α i | < r. If Λ = {Λ 1 , . . . , Λ r } is a basis of a r-dimensional space D in C[d] r-1 with Λ i = |α|<r µ i,α d α
α! , the Plücker coordinates of D are, up to a scalar, of the form ∆ α1,...,αr = det µ i,αj 1≤i,j≤r . In particular, a monomial set B = {x β1 , . . . ,

x βr } ⊂ C[x] r-1 has a dual basis in D iff ∆ β1,...,βr (D) = 0. If (B = {x βi } r i=1 , Λ = {Λ i } r i=1
) is a graded primal-dual basis pair, then µ i,βj = δ i,j . To keep our notation consistent with the previous sections, the coordinates of Λ i ∈ Λ when Λ is dual to B will be denoted by µ βi,α instead of µ i,α . By properties of the determinant, the Plücker coordinates of D are such that [START_REF] Hauenstein | On deflation and multiplicity structure[END_REF] µ βi,α = ∆ β1,...,βi-1,α,βi+1,...,βr ∆ β1,...,βr i = 1, . . . , r. For D ∈ H r , for t ≥ 0 we denote by D t the vector space of elements of D of order ≤ t. We verify that (2) The dual basis Λ = {Λ 1 , . . . , Λ r } of B satisfies Λ 1 = 1 and the equations (8), ( 9) and (10).

If D is the dual of an ideal Q = D ⊥ ⊂ C[x] and B = {x β1 , . . . , x βr } is a basis of C[x]/Q so that ∆ β1,...,βr (D) = 0, the normal form of x α ∈ C[x] r-1 modulo Q = D ⊥ in the basis B is N (x α ) =
3.2. D ∈ H r iff D ⊥ = Q is an m-primary ideal such that dim C[x]/Q = r. Proof. Let D ∈ H r . We prove that D ⊥ = Q is an m-primary ideal. As D is stable by derivation, Q = D ⊥ is an ideal of C[x]. This also implies that 1 ∈ D, so that Q ⊂ m. As dim D = dim C[x]/Q = r, δ = ord(D) is finite and m δ+1 ⊂ D ⊥ = Q. Therefore, Q is m-primary,
D ⊥ t = D ⊥ + m t+1
Proof

∂ d k (Λ j ) = rt-1 j=1 ν i,k Λ i for t = 1, . . . , o
. By Proposition 2.3, Equations ( 8) and ( 9) are satisfied. As B t is dual to Λ 1 , . . . , Λ rt , Equation ( 10) are satisfied.

(

) ⇒ (1) Let Λ i ∈ C[d] r-1 for i = 1, . . . , r be elements of C[d] 2
r-1 dual to B, which satisfies Equations ( 8), ( 9) and [START_REF] Emiris | The DMM bound: multivariate (aggregate) separation bounds[END_REF]. By induction on t = 0, . . . , δ = deg(B), we prove that if

Λ t = {Λ 1 , . . . , Λ rt } is dual to B t , then Λ 1 , . . . , Λ rt ∈ C[d] t .
The property is true for t = 0 since Λ 1 = 1. If it is true for t -1, for Λ j with j = r t-1 + 1, . . . , r t we have by ( 8 For a sequence h = (h 0 , h 1 , . . . , h δ ) ∈ N δ+1 + and 0 ≤ t ≤ δ, let h t = (h 0 , . . . , h t ), r t = t i=0 h i . For r ≥ 1 we denote by S r the set of sequences h of some length δ < r with h i = 0, h 0 = 1 and r δ = r. For h ∈ S r , we consider the following subvarieties of H rt :

H ht = {D ∈ H rt | dim D i = dim D ∩ C[d] i ≤ r i , i = 0, . . . , t}.
These are projective varieties in H rt defined by rank conditions on the linear spaces D ∩ C[d] i for D ∈ H rt , that can be expressed in terms of homogeneous polynomials in the Plücker coordinates of D. In particular, the varieties H h := H h δ are projective subvarieties of H r . They may not be irreducible or irreducible components of H r , but we have H r = ∪ h∈S r H h .

We will study a particular component of H h , that we call the regular component of

H h , denoted H reg h . It is characterized as follows. Let H reg h0 = { 1 } = {C[d] 0 } = G 1 (C[d] 0
) and assume that H reg ht-1 has been defined as an irreducible component of H ht-1 . Let

Wt = {(Dt-1, Et) | Dt-1 ∈ H h t-1 , Et ∈ Gr t (C[d]t), Dt-1 ⊂ Et, ∀i ∂ d i Et ⊂ Dt-1}
The constraints D t-1 ⊂ E t and ∂ di E t ⊂ D t-1 for i = 1, . . . , n define a linear system of equations in the Plücker coordinates of E t (see e.g. [START_REF] Doubilet | On the Foundations of Combinatorial Theory[END_REF]), corresponding to the equations ( 5), [START_REF] Dayton | Multiple zeros of nonlinear systems[END_REF]. By construction, the projection of

W t ⊂ H ht-1 × G rt (C[d] t ) on the second factor G rt (C[d] t ) is π 2 (W t ) = H ht and the projection on the first factor is π 1 (W t ) = H ht-1 .
There exists a dense subset U t-1 of the irreducible variety H reg ht-1 (with U t-1 = H reg ht-1 ) such that the rank of the linear system corresponding to ( 5) and ( 6) defining E t is maximal. Since π -1 1 (D

Rational parametrization

Let B = {x β1 , . . . , x βr } ⊂ C[x] r-1 be a monomial set. In this section we assume that B is closed under division and its monomials are ordered by increasing degree. For t ∈ N, we denote by

B t = B ∩ C[x] t , by B [t] the subset of its monomials of degree t. Let h t = |B [t] |, r t = 0≤i≤t h t = |B t | and δ = deg(B).
Let Since for D ∈ H B we have ∆ B (D) = 0, we can define the affine coordinates of H B using the coordinates of the elements of the basis Λ = {Λ 1 , . . . , Λ r } dual to B: µ βj ,α = ∆ β1,...,βj-1,α,βj+1,...,βr ∆ B : j = 1, . . . , r, |α| < r .

H B := {D ∈ H r | B t is a basis of C[x]/(D ⊥ + m t+1 ), t = 0, . . . ,
The following lemma shows that the values of the coordinates {µ βi,βj+e k : i, j = 1, . . . r, Proof. Since D is closed under derivation, by Proposition 2.3 there exist

|β j | < |β i |, k = 1, . . . , n} uniquely define Λ.
c i,s,k ∈ C such that ∂ d k (Λ i ) = |βs|<|βi| c i,s,k Λ s . Then µ β i ,β j +e k = Λi(x β j +e k ) = ∂ d k (Λi)(x β j ) = |βs|<|β i | c i,s,k Λs(x β j ) = c i,j,k .
The second claim follows from obtaining the coefficients in Λ recursively from Λ 1 = 1 and

Λ i = |βj |<|βi| n k=1 µ βi,βj +e k ∫ k Λ j for i = 2, . . . , r.
We define µ := {µ βi,βj+e k } i,j=1,...r,|βj|<|βi|,k=1,...,n , µ t := {µ βi,βj+e k ∈ µ :

|β i | ≤ t} ⊂ µ and µ [t] := {µ βi,βj +e k ∈ µ : |β j | = t} ⊂ µ t .
The next definition uses the fact that Equations ( 13) and ( 14) are linear in ν k j with coefficients depending on µ t-1 :

Definition 4.2. Given D t-1 ∈ H Bt-1 with a unique basis Λ t-1 = {Λ 1 , . . . , Λ rt-1 } with Λ i = |α|<t µ βi,α d α α!
for j = 1, . . . , r t-1 that is dual to B t-1 , uniquely determined by µ t-1 = {µ βi,βj+e k :

|β i | ≤ t -1, |β j | < |β i |}
as above. Recall from Notation 2.6 that H t is the coefficient matrix of the homogeneous linear system ( 13) and ( 14) in the variables {ν k j : j = 1, . . . , r t-1 , k = 1, . . . , n}. To emphasize the dependence of its coefficients on D t-1 or µ t-1 we use the notation H t (D t-1 ) or H t (µ t-1 ). For D ∈ H reg h in an open subset, the rank ρ t of H t (D t-1 ) is maximal.

The next definition describes a property of a monomial set B such that it will allow us to give a rational parametrization of H B . Definition 4.3. For t = 1, . . . , δ = deg(B) we say that can be parametrized by rational functions of free parameters µ. We present hereafter Algorithm 2 to compute such a parametrization iteratively.

D t ∈ G rt (C[d] t ) is regular for B t if, • dim(D t ) = r t = |B t |, • rank H t (D t-1 ) = ρ t the generic rank of H t on H reg ht , • ∆ B [t] (D [t] ) = 0 where ∆ B [t] (D [t] ) is the Plücker coordinate of D [t] ∈ G ht (C[d] r ) corresponding

Algorithm 2 Rational Parametrization -Iteration t

Input: t > 0, B t = {x β1 , . . . , x βr t } ⊂ C[x] t closed under division and regular, µ t-1 ⊂ µ t-1 and Φ t-1 : µ t-1 → q βj,α (µ t-1 ) |βj |≤t-1,|α|<r with q βj ,α ∈ Q(µ t-1 ) parametrizing a dense subset of H reg ht-1 . Output: µ t ⊂ µ t and Φ t : µ t → q βj,α |βj |≤t,|α|<r , q βj ,α ∈Q(µ t ) extending Φ t-1 and parametrizing a dense subset of

H reg ht . (1) Let H t be as in Notation 2.6, ν = [ν k j : j = 1, . . . , r t-1 , k = 1, . . . , n] T . Decompose H t (Φ t-1 (µ t-1 )) • ν = 0 as (17) A(µ t-1 ) B(µ t-1 ) C(µ t-1 )   ν ′ ν ′′ ν   = 0,
where ν ′ is associated to a maximal set of independent columns of H t (Φ t-1 (µ t-1 )), ν ′′ = {ν k j : x βj +e k ∈ B [t] } and ν refers to the rest of the columns. If no such decomposition exists, return "B t is not regular".

(2) For ν k j ∈ ν ′ express ν k j = ϕ k j (ν, ν ′′ ) ∈ Q(µ t-1 )[ν, ν ′′ ] 1 as the generic solution of the system H t (Φ t-1 (µ t-1 )) • ν = 0. (3) For i = r t-1 + 1, . . . , r t do: (3.1) Define µ [t],i := µ βi,βj+e k : ν j,k ∈ ν , µ ′ [t],i = {µ βi,βj+e k : ν k j ∈ ν ′ }, µ ′′ [t],i = {µ βi,βj+e k : ν k j ∈ ν ′′ }, and µ t := µ t-1 ∪ rt i=rt-1+1 µ [t],i . (3.2) For µ βi,βj+e k ∈ µ ′′ [t],i set q βi,βj+e k = µ βi,βj+e k = 1 if β i = β j + e k and 0 otherwise. (3.3) For µ βi,βj+e k ∈ µ ′ [t],i define q βi,βj+e k := ϕ k j (µ [t],i , µ ′′ [t],i ) ∈ Q(µ t ) (3.4)
For |α| < r and µ βi,α ∈ µ t find q βi,α using Lemma 4.1.

Proposition 4.4. Let B = {x β1 , . . . , x βr } ⊂ C[x] r-1 be closed under division and assume that B is a regular basis. There exist a subset µ ⊂ µ with |µ| = γ and rational functions q βj,α (µ) ∈ Q(µ) for j = 1, . . . , r and |α| < r, such that the map Φ : C γ → H B defined by , and a rational parametrization of a basis Λ 1 (µ t ), . . . , Λ rt (µ t ) of a generic element of H reg Bt . For t = 0, we define Λ 1 = 1 and µ 0 = ∅. Assume that there exist µ t-1 ⊂ µ t-1 and a rational parametrization Λ 1 (µ t-1 ), . . . , Λ rt-1 (µ t-1 ) of a basis dual to B t-1 for a generic element H Bt-1 defined by the map

Φ : µ → q βj,α (µ) j=1,...,
Φ t-1 : µ t-1 → q βj ,α (µ t-1 ) |βj|≤t-1,|α|<r .
This means that im Φ t-1 = H Bt-1 . Denote by

D t-1 (µ t-1 ) ∈ G rt-1 (Q(µ t-1 )[d] t-1
) the space spanned by {Λ 1 (µ t-1 ), . . . , Λ rt-1 (µ t-1 )} over the fraction field Q(µ t-1 ). By Theorem 3.3 and Lemma 2.5, to define µ t and to extend D t-1 (µ t-1 ) to D t (µ t ), we need to find Λ rt-1+1 , . . . , Λ rt of the form

Λ i = rt-1 j=1 n k=1 µ βi,βj +e k ∫ k Λ j (µ t-1 ) i = r t-1 + 1, . . . , r t ,
satisfying the system of equations ( 13) and ( 14), i.e. such that

Λ i ∈ ker H t (µ t-1 ) for i = r t-1 + 1, . . . , r t ,
where H t (µ t-1 ) = H t Φ t-1 (µ t-1 ) and Equations ( 12) are satisfied. Since B is a regular basis, the kernel of H t (µ t-1 ) over Q(µ t-1 ) contains a subspace

D [t] of dimension h t = |B [t] | with ∆ B [t] (D [t] ) = 0. Therefore, the systems H t (µ t-1 ) ν = 0 with ν = [ν k j : j = 1, . . . , r t-1 , k = 1, . .

. , n]

T can be decomposed as ( 18)

A(µ t-1 ) B(µ t-1 ) C(µ t-1 )   ν ′ ν ′′ ν   = 0,
where ν ′ is associated to a maximal set of independent columns of H t (µ t-1 ), ν ′′ = {ν k j :

x βj +e k ∈ B [t]
} and ν is associated to the remaining set of columns. Note that |ν| = dim(ker H t (µ t-1 ))h t . Thus, ν ′′ ∪ ν is the set of free variables of the homogeneous system H t (µ t-1 ) ν = 0 and a general solution is such that the variables in ν ′ are linear functions of the variables in ν ′′ and ν, with rational coefficients in µ t-1 .

We obtain the coefficients of Λ rt-1+1 , . . . , Λ rt that satisfy equations ( 13) and ( 14) and ( 12) from the general solutions of H t (µ t-1 ) ν = 0 by further specializing the variables in ν ′′ to 0's and 1's, according the duality conditions. Define

µ [t],i := µ βi,βj +e k : ν j,k ∈ ν ⊂ µ [t]
. Thus, the parameters in µ [t] are linear functions of µ [t],i with rational coefficients in µ t-1 . The denominator in these coefficients is a factor of the numerator of a maximal non-zero minor of A(µ t-1 ). Note that the rest of the coefficients of Λ i are polynomial functions of the parameters µ t-1 ∪ µ [t] by Lemma 4.1. Define

µ t := µ t-1 ∪ rt i=rt-1+1 µ [t],i .
Thus, we get a parametrization of the coefficients of Λ rt-1+1 (µ t ), . . . , Λ rt (µ t ) in terms of µ t , which defines the degree t part of the map Φ t : µ t → (q βj ,α (µ t )) |βj|≤t,|α|<r . For D t ∈ H Bt , the coefficients of its basis dual to B t can be parametrized by Φ t parameter values µ t such that a maximal non-zero minor of A(µ t-1 ) in Q(µ t-1 ) does not vanish.

that the number of new parameters introduced is

|µ t \ µ t-1 | = (r t -r t-1 ) • |µ [t],i | = h t dim ker H t (µ t-1 ) -h t which is equal to γ [t] = dim G ht (ker H t (µ t-1 )) = dim G ht (ker H t (D t-1 )) for D t-1 generic in U t-1 as claimed.
To prove that Φ t parametrizes a dense subset of the projective variety H reg ht , note that the image im(Φ t ) of Φ t is a subset of H ht , the Zariski closure V t of im(Φ t ) is an irreducible subvariety of H ht . Furthermore, its projection

π t-1 (V t ) ⊂ H ht-1 is the closure of the image of im(Φ t-1 ) since if D t = imΦ t (µ * t ) then D t-1 = D t ∩ C[d] t-1 = Φ t-1 (µ * t-1 ). By induction hypothesis, π t-1 (V t ) = imΦ t-1 = H reg ht-1 . Thus, V t is the irreducible component of H ht which projects onto H reg ht-1
, that is H reg ht . Definition 4.5. We denote by H t (µ) a maximal square submatrix of A in [START_REF] Iarrobino | Punctual Hilbert Schemes[END_REF] such that det(H t (µ t-1 )) = 0.

The size of H t (µ) is the size of ν ′ in [START_REF] Iarrobino | Punctual Hilbert Schemes[END_REF], that is the maximal number of independent columns in H t (µ t-1 ).

Given an element

D = Λ 1 ∧ • • • ∧ Λ r ∈ G r (C[d] r-1
), in order to check that D is regular for B, it is sufficient to check first that ∆ B (D) = 0 and secondly that |H t (µ)| = 0 for all t = 0, . . . , δ, where µ = (µ β,α ) is the ratio of Plücker coordinates of D defined by the formula (16).

Newton's iterations

In this section we describe the extraction of a square, deflated system that allows for a Newton's method with quadratic convergence. We assume that the sole input is the equations f = (f 1 , . . . , f N ) ∈ C[x] N , an approximate point ξ ∈ C n and a tolerance ε > 0.

Using this input we first compute an approximate primal-dual pair (B, Λ) by applying the iterative Algorithm 1. The rank and kernel vectors of the matrices K t (see Algorithm 1) are computed numerically within tolerance ε, using SVD. Note that here and in Section 6 we do not need to certify the SVD computation but we are only using SVD to certify that some matrices are full rank by checking that the distance to the variety of singular matrices is bigger than the perturbation of the matrix. Thus we need a weaker test, which relies only on a lower bound of the smallest singular value.

The algorithm returns a basis B = {x β1 ξ , . . . , x βr ξ } with exponent vectors E = {β 1 , . . . , β r }, as well as approximate values for the parameters µ = {µ βi,βj+e k :

|β j | < |β i | ∈ E, k = 1, .
. . , n}. These parameters will be used as a starting point for Newton's iteration. Note that, by looking at B, we can also deduce the multiplicity r, the maximal order δ of dual differentials, the sequences r t = |B t |, and h t = |B [t] | for t = 0, . . . , δ.

Let F be the deflated system with variables (x, µ) defined by the relations (8) and Equations ( 9), ( 10) and ( 12) i.e.

F (x,µ)=                |βs|<|β j |<|β i | µ β i ,β j +e k µ β j ,βs+e l -µ β i ,β j +e l µ β j ,βs+e k =0 (a) for all i = 1, . . . , r, |βs| < |β i |, k = l ∈ {1, . . . , n} µ β i ,β j +e k = 1 for β i = β j + e k 0 for β j + e k ∈ E, β i = β j + e k , (b)
Λ i (f j ) = 0, i = 1, . . . , r, j = 1, . . . , N.

(c)

Here Λ 1 = 1 x and Λ i = |βj |<|βi|

n k=1 µ βi,βj+e k ∫ k Λ j ∈ C[µ][d x ]
denote dual elements with parametric coefficients defined recursively. Also, if Λ i = |α|≤|βi| µ βi,α

d α x α! then Λ i (f j ) = |α|≤|βi| µ βi,α ∂ α (f j )(x) α!
which is in C[x, µ] by Lemma 4.1. Note, however, that (a) and (b) are polynomials in C[µ], only (c) depends on x and µ. Equations (b) define a simple substitution into some of the parameters µ. Hereafter, we explicitly substitute them and eliminate this part (b) from the equations we consider and reducing the parameter vector µ. By Theorem 2.8, if B is a graded primal basis for f at the root ξ * then the above overdetermined system has a simple root at a point (ξ * , µ * ).

To extract a square subsystem defining the simple root (ξ * , µ * ) in order to certify the convergence, we choose a maximal set of equations whose corresponding rows in the Jacobian are linearly independent. This is done by extracting first a maximal set of equations in (a) with linearly independent rows in the Jacobian. For that purpose, we use the rows associated to the maximal invertible matrix H t (Definition 4.5) for each new basis element Λ i ∈ D [t] and t = 1, . . . , r. We denote by G 0 the subsystem of (a) that correspond to rows of H t .

We complete the system of independent equations G 0 with equations from (c), using a QR decomposition and thresholding on the transposed Jacobian matrix of G 0 and (c) at the approximate root. Let us denote by F 0 the resulting square system, whose Jacobian, denoted by J 0 , is invertible.

For the remaining equations F 1 of (c), not used to construct the square system F 0 , define Ω = {(i, j) : Λ i (f j ) ∈ F 1 }. We introduce new parameters ǫ i,j for (i, j) ∈ Ω and we consider the perturbed system

f i,ǫ = f i - j|(i,j)∈Ω ǫ i,j x βj ξ .
The perturbed system is f ǫ = fǫ B, where ǫ is the N × r matrix with [ǫ] i,j = ǫ i,j if (i, j) ∈ Ω and [ǫ] i,j = 0 otherwise. Denote by F (x, µ, ǫ) obtained from F (x, µ) by replacing Λ j (f i ) by Λ j (f i,ǫ ) for j = 1, . . . , r, i = 1, . . . , N . Then the equations used to construct the square Jacobian J 0 are unchanged. The remaining equations are of the form Λ j (f i,ǫ ) = Λ j (f i )ǫ i,j = 0 (i, j) ∈ Ω. Therefore the Jacobian of the complete system F (x, µ, ǫ) is a square invertible matrix of the form

J ǫ := J 0 0 J 1 Id
where J 1 is the Jacobian of the system F 1 of polynomials Λ j (f i ) ∈ C[x, µ] with (i, j) ∈ Ω. Since J ǫ is invertible, the square extended system F (x, µ, ǫ) has an isolated root (ξ * , µ * , ǫ * ) corresponding to the isolated root (ξ * , µ * ) of the square system F 0 . Furthermore, Λ * j (f i ) = ǫ * i,j = 0 for (i, j) ∈ Ω. Here Λ * 1 , . . . , Λ * r ∈ C[d ξ * ] are defined from (ξ * , µ * ) recursively by

Λ * 1 = 1 ξ * and Λ * i = |βj|<|βi| n k=1 µ * βi,βj+e k ∫ k Λ * j . (19)
We have the following property: Theorem 5.1. If the Newton iteration

(ξ k+1 , µ k+1 ) = (ξ k , µ k ) -J0(ξ k , µ k ) -1 F0(ξ k , µ k ),
starting from a point (ξ 0 , µ 0 ) converges when k → ∞, to a point (ξ * , µ * ) such that B is a regular basis for the inverse system D * associated to (ξ * , µ * ) and D * is complete for f , then there exists a perturbed system f i,ǫ * = f i -j|(i,j)∈Ω ǫ * i,j x βj ξ * with ǫ * i,j = Λ * j (f i ) such that ξ * is a multiple root of f i,ǫ * with the multiplicity structure defined by µ * .

Proof. If the sequence (ξ k , µ k ) converges to the fixed point (ξ * , µ * ), then we have F 0 (ξ * , µ * ) = 0 and in particular, G 0 (ξ * , µ * ) = 0 where G 0 (ξ * , µ * ) = 0 is the subset of equations selected from (a).

As µ * is regular for B, if it satisfies G 0 (ξ * , µ * ) = 0, it must satisfy all equations (a). Therefore µ * defines a point

D * = Λ * 1 ∧ • • • ∧ Λ * r ∈ H reg B . As (Λ * i
) is a basis of D * dual to B and f i,ǫ * = f ij:(i,j)∈Ω ǫ * i,j x βj ξ * with ǫ * i,j = Λ * j (f i ) for (i, j) ∈ Ω, we have that if (i, j) ∈ Ω then Λ * j (f i,ǫ * ) = Λ * j (f i )ǫ * i,j = 0. Otherwise Λ * j (f i,ǫ * ) = Λ * j (f i ), since it is one of the equations selected in (c) to construct the system F 0 and F 0 (ξ * , µ * ) = 0. This shows that

f ǫ * = (f i,ǫ * ) N i=1 ⊂ (D * ) ⊥ .
Since f ǫ * is obtained from f by adding elements in B, the system (c), at order δ+1 for f ǫ * and f are equivalent. Thus D * is complete for f and f ǫ and D * = (f ǫ * ) ⊥ ∩ C[d ξ * ] is the inverse system at ξ * of the system f ǫ * .

Certification

In this section we describe how to certify that the Newton iteration defined in Section 5 quadratically converges to a point that defines an exact root with an exact multiplicity structure of a perturbation of the input polynomial system f . More precisely, we are given f

= (f 1 , . . . , f N ) ∈ C[x] N , B = {x β1 , . . . , x βr } ⊂ C[x]
in increasing order of degrees and closed under division, δ := |β r |. We are also given the deflated systems F (x, µ), its square subsystem F 0 (x, µ) defined in Section 5 and F 1 (x, µ) the remaining equations in F (x, µ). Finally, we are given ξ 0 ∈ C n and µ 0 = {µ

(0) βi,βj+e k ∈ C : i, j = 1, . . . , r, |β j | < |β i |, k = 1, . . . , n}.
Our certification will consist of a symbolic and a numeric part: Regularity certification. We certify that B is regular (see Definition 4.3). This part of the certification is purely symbolic and inductive on t. Suppose for some t -1 < δ we certified that B t-1 is regular and computed the parameters µ t-1 and the parametrization Φ t-1 : µ t-1 → q βi,α (µ t-1 ) |βi|≤t-1,|α|≤t-1 (Algorithm 2). Then to prove that B t is regular, we consider the coefficient matrix H t of equations ( 13) and ( 14). We substitute the parametrization Φ t-1 to get the matrices H t (µ t-1 ). We symbolically prove that the rows of H t (µ t-1 ) (Definition 4.5) are linearly independent and span all rows of H t (µ t-1 ) over Q(µ t-1 ). If that is certified, we compute the parameters µ t and the parametrization Φ t : µ t → (q βi,α (µ t )) |βi|≤t,|α|≤t as in Algorithm 2 inverting the square submatrix H t of H t such that the denominators of q βi,α for |β i | = t divide det(H t (µ t-1 )) = 0.

Singularity certification.

(C1) We certify that the Newton iteration for the square system F 0 starting from (ξ 0 , µ 0 ) quadratically converges to some root (ξ * , µ * ) of F 0 , such that (ξ 0 , µ 0 ) -(ξ * , µ * ) 2 ≤ β, using α-theory. (C2) We certify that D * = span(Λ * ) is regular for B (see Definition 4.3), by checking that |H t (µ * )| = 0 for t = 1, . . . , δ (See Definition 4.5), using the Singular Value Decomposition of H t (µ 0 ) and the distance bound β between µ * and µ 0 . (C3) We certify that Λ * is complete for f at ξ * (see Definition 2.7), where Λ * ⊂ C[d ξ * ] is the dual systems defined from (ξ * , µ * ) recursively as in [START_REF] Lee | On isolation of singular zeros of multivariate analytic systems[END_REF]. This is done by checking that ker K δ+1 (ξ * , µ * ) = {0} (See Definition 2.7), using the Singular Value Decomposition of K δ+1 (ξ 0 , µ 0 ) and the distance bound β between (ξ * , µ * ) and (ξ 0 , µ 0 ).

Let us now consider for a point-multplicity structure pair (ξ 0 , µ 0 ) γ :

= sup k≥2 DF -1 0 (ξ 0 , µ 0 ) D k F0(ξ0,µ0) k! 1 k-1 , β := 2 DF -1
0 (ξ 0 , µ 0 ) F 0 (ξ 0 , µ 0 ) , α := β γ and for a matrix function A(ξ, µ), let L 1 (A; ξ 0 , µ 0 ; b) be a bound on its Lipschitz constant in the ball B b (ξ 0 , µ 0 ) of radius b around (ξ 0 , µ 0 ) such that A(ξ, µ) -A(ξ 0 , µ 0 ) ≤ L 1 (A; ξ 0 , µ 0 ; b) (ξ, µ)-(ξ 0 , µ 0 ) for (ξ, µ) ∈ B b (ξ 0 , µ 0 ). For a matrix M , let σ min (M ) be its smallest singular value. We have the following result: Theorem 6.1. Let B = {x β1 , . . . , x βr } ⊂ C[x] be closed under division and suppose B is regular. Suppose that α < α0 := 0.26141, L 1 (K δ+1 ; ξ 0 , µ 0 ; β) β < σ min (K δ+1 (ξ 0 , µ 0 )) and for t = 1, . . . , δ it holds that L 1 (H t ; µ 0 ; β) β < σ min (H t (µ 0 )). Then the Newton iteration on the square system F 0 starting from (ξ 0 , µ 0 ) converges quadratically to a point (ξ * , µ * ) corresponding to a multiple point ξ * with multiplicity structure µ * of the perturbed system f ǫ * = fǫ * B ξ * such that ǫ * ≤ F 1 (ξ 0 , µ 0 ) + L 1 (F 1 ; ξ 0 , µ 0 ; β) β, where B ξ * = {x β1 ξ * , . . . , x βr ξ * }. Proof. By the α-theorem [START_REF] Blum | Complexity and Real Computation[END_REF][Chap. 8, Thm. 1], the Newton iteration on F 0 starting from (ξ 0 , µ 0 ) converges quadratically to a point (ξ * , µ * ) such that (ξ * , µ * ) -(ξ 0 , µ 0 ) < β .

We deduce that

K δ+1 (ξ * , µ * ) -K δ+1 (ξ0, µ0) ≤ L1(K δ+1 ; ξ0, µ0; β) (ξ * , µ * ) -(ξ0, µ0) < σmin(K δ+1 (ξ0, µ0)).
Therefore K δ+1 (ξ * , µ * ) is within a ball around K δ+1 (ξ 0 , µ 0 ) of matrices of maximal rank, since σ min (K δ+1 (ξ 0 , µ 0 )) is the distance between K δ+1 (ξ 0 , µ 0 ) and the set of matrices not of maximal rank. Thus ker K δ+1 (ξ * , µ * ) = {0}. A similar argument shows that |H t (µ * )| = 0 for t = 1, . . . , δ. By Theorem 5.1, (ξ * , µ * ) defines a multiple root ξ * with multiplicity structure µ * for the perturbed system

f ǫ * = f -ǫ * B ξ * with ǫ * = F1(ξ * , µ * ) ≤ F1(ξ0, µ0) + F1(ξ * , µ * ) -F1(ξ * , µ * ) ≤ F1(ξ0, µ0) + L1(F1; ξ0, µ0; β) (ξ * , µ * ) -(ξ * 0 , µ * 0 ) ≤ F1(ξ0, µ0) + L1(F1; ξ0, µ0; β) β.

Experimentation

In this section we work out some examples with (approximate) singularities. The experiments are carried out using Maple, and our code is publicly available at https://github.com/filiatra/polyonimo. Example 7.1. We consider the equations

f 1 = x 3 1 + x 2 2 + x 2 3 -1, f 2 = x 3 2 + x 2 1 + x 2 3 -1, f 3 = x 3 3 + x 2 1 + x 2
2 -1, the approximate root ξ 0 = (0.002, 1.003, 0.004) and threshold ε = 0.01. In the following we use 32-digit arithmetic for all computations. both of them, which is sufficient, since the involved Lipschitz and γ constants are of the order of 1 for the input polynomials. In the first iteration we obtain β ≈ 0.00011 which clearly indicates that we are in the region of convergence. Indeed, the successive residuals for 4 iterations are 0.00603, 4.0•10 -5 , 2.07•10 -9 , 8.6•10 -18 , 3.55• 10 -35 . Clearly, the residual shrinks with a quadratic rate1 . We obtain ξ 4 = (1.8 • 10 -37 , 1.0, 2.8 • 10 -36 ) and the overdetermined system is satisfied by this point: F (ξ 4 , µ 4 ) ∞ = 8 • 10 -35 ; the resulting dual structure is D * 2 = {1, d 1 , d 3 , d 1 d 3 }. Example 7.2. We demonstrate how our method handles inaccuracies in the input, and recovers a nearby system with a true multiple point. Let

f 1 = x 1 2 + x 1 -x 2 + 0.003 , f 2 = x 2 2 + 1.004x 1 -x 2 .
There is a cluster of three roots around ξ 0 = (0.001, -0.002). Our goal is to squeeze the cluster down to a three-fold real root. We use 32 digits for the computation. Starting with ξ 0 , and a tolerance equal to 10 -2 Algorithm 1 produces an approximate dual 1, d 1 + 1.00099651d 2 , d 2 1 + 1.00099651d 1 d 2 + 1.00266222d 2 2 + 0.99933134d 2 and identifies the primal basis B = {1, x 1 , x 2 1 } using pivoting on the integration matrix. The sole stability condition reads µ 1,1µ 2,2 = 0, and

Λ 1 = 1, Λ 2 = d 1 + µ 1,1 d 2 , Λ 3 = d 2 1 + µ 1,1 d 1 d 2 + µ 2,1 d 2 + µ 2,2 µ 1,1 d 2
2 . The nearby system that we shall obtain is deduced by the residue in Newton's method. In particular, starting from ξ 0 , we consider the square system given by removing the equations Λ 1 (f 1 ) = 0 and Λ 2 (f 2 ) = 0. The rank of the corresponding Jacobian matrix remains maximal, therefore such a choice is valid. Newton's iterations converge quadratically to the point (ξ 5 , µ 5 ) = (1.1 • 10 -33 , 1.2 • 10 -33 , 1, 1, 1). The full residual is now F (ξ 5 , µ 5 ) = (0, 0.003, -10 -32 , 10 -32 , 0.004, 0, 0) .

This yields a perturbation f1 ≈ f 1 -0.003 and f2 ≈ f 2 -0.004(x 1ξ * 1 ) to obtain a system with an exact multiple root at the origin (cf. Th. 6.1). Of course, this choice of the square sub-system is not unique. By selecting to remove equations Λ 1 (f 1 ) = 0 and Λ 1 (f 2 ) = 0 instead, we obtain (ξ 5 , µ 5 ) = (0.00066578, -0.00133245, 1.001, 1.0, 1.001) and the residual F (ξ 5 , µ 5 ) = (0, 0.005, 0.002, 0, 0, 0, 0), so that the nearby system f * 1 ≈ x 1 2 + x 1x 2 + 0.008, f * 2 ≈ x 2 2 + 1.004x 1x 2 + 0.002 has a singularity at the limit point ξ * ≈ (0.00066578, -0.00133245) described locally by the coefficients µ * ≈ (1.001, 1.0, 1.001). Finally, consider the two square sub-systems as above, after changing f 1 , f 2 to define an exact three-fold root at the origin (i.e. f 1 = x 1 2 + x 1x 2 , f 2 = x 2 2 + x 1x 2 ). Newton's iteration with initial point ξ 0 on either deflated system converges quadratically to (ξ, µ) = (0, 1). This is a general property of the method: exact multiple roots and their structure are recovered by this process if ξ 0 is a sufficiently good initial approximation (cf. Section 5). We plan to develop this aspect further in the future.

Example 7.3. We show some execution details on a set of benchmark examples in taken from [START_REF] Dayton | Computing the multiplicity structure in solving polynomial systems[END_REF], see also [START_REF] Mantzaflaris | Singular zeros of polynomial systems[END_REF]. For this benchmark, we are given systems and points with multiplicities. We perturb the given points with a numerical perturbation of order 10 -2 . We use double precision arithmetic and setup Newton's iteration; with less than 10 iterations, the root was approximated within the chosen accuracy.

In Table 1, "IM" is the maximal size of the (numeric) integration matrix that is computed to obtain the multiplicity, "#µ" is the number of new parameters that are needed for certified deflation, "SC" is the number of stability constraints that were computed and "OS" stands for the size of the overdetermined system (equations × variables). This is the size of the Jacobian matrix that must be computed and inverted in each Newton's iteration. We can observe that the number of parameters required can grow significantly. Moreover, these parameters induce non-trivial denominators in the rational functions q βj ,α (µ) of Prop. [START_REF] Brianc ¸on | Description de Hilb n C{x, y}[END_REF]

  This lemma shows that to compute Q ⊥ , it suffices to compute all polynomials of C[d ξ ] which are in I ⊥ . Let us denote this set D = I ⊥ ∩ C[d ξ ]. It is a vector space stable under the derivations ∂ d i,ξ . Its dimension is the dimension of Q ⊥ or C[x]/Q, that is the multiplicity of ξ, denoted r ξ (I), or simply r if ξ and I is clear from the context. For an element Λ(d ξ ) ∈ C[d ξ ] we define the degree or order ord(Λ) to be the maximal |β| s.t. d β ξ appears in Λ(d ξ ) with non-zero coefficient.

Lemma 2 . 2 (( 1 )

 221 Graded primal-dual basis pair). Let f , ξ, Q, D, r = r ξ (f ) and δ = δ ξ (f ) be as above. Then there exists a primal-dual basis pair (B, Λ) of the local ring C[x]/Q with the following properties: The primal basis of the local ring C[x]/Q has the form (4) B := x β1 ξ , x β2 ξ , . . . , x βr ξ .

Lemma 2 . 5 .

 25 Let B = {x β1 ξ , . . . , x βr ξ } ⊂ C[x] be closed under division and ordered by degree. Let Λ = {Λ 1 , . . . , Λ r } ⊂ C[d ξ ] be a linearly independent set such that

  Fixing a basis e 1 , . . . , e s of V , the Plücker coordinates of v are the coefficients of ∆ i1,...,ir (v) of v = i1<•••<ir ∆ i1,...,ir (v) e i1 ∧ • • • ∧ e ir . When V = C[d] r-1 , a natural basis is the dual monomial basis ( d α α! ) |α|<r . The Plücker coordinates of an element v ∈ G r (C[d] r-1 ) for this basis are denoted ∆ α1,...,αr (v) where

∆Definition 3 . 1 .

 31 β1,...,βj-1,α,βj+1,...,βr ∆ β1,...,βrx βj .(if deg(x α ) ≥ r, then N (x α ) = 0). Let H r ⊂ G r (C[d]r-1 ) be the set of linear spaces D of dimension r in C[d] r-1 which are stable by the derivations ∂ di with respect to the variables d (i.e. ∂ di D ⊂ D for i = 1, . . . , n). We called H r the punctual Hilbert scheme of points of multiplicity r. If D ⊂ C[d] is stable by the derivations ∂ di , then by duality I = D ⊥ ⊂ C[x] is a vector space of C[x] stable by multiplication by x i , i.e. an ideal of C[x].

Proposition

  which shows the first implication. Conversely, let Q be a m-primary ideal such that dim C[x]/Q = r. Then by Lemma 2.1, D = Q ⊥ ⊂ C[d] t is stable by derivation and of dimension r = dim C[x]/Q. Thus D ∈ H r . This concludes the proof of the proposition.

  . The next theorem follows from Proposition 2.3 and Lemma 2.5. Theorem 3.3. For B ⊂ C[x] closed under division such that |B| = r and δ = deg(B), the following points are equivalent: (1) D ∈ H r and B t is a basis of C[x]/(D ⊥ + m t+1 ) for t = 1, . . . , δ.

  ),[START_REF] Doubilet | On the Foundations of Combinatorial Theory[END_REF] and Proposition 2.3, that∂ d k (Λ j ) = rt-1 j=1 ν i,k Λ i , k = 1, . . . , n. Thus Λ j ∈ C[d] t. This shows that D t is stable by derivation whereD t ⊂ C[d] t is the vector space spanned Λ 1 , . . . , Λ rt ∈ C[d] t . Let D = D δ . Since, by (10), B t is dual to Λ 1 , . . . , Λ rt ∈ C[d] t , we see that D∩C[d] t = D t . By Proposition 3.2, Q = D ⊥ is a m-primary ideal such that dim C[x]/Q = dim D = |B| = r. Moreover, since B t is dual to the basis {Λ 1 , . . . , Λ rt } of D t , B t is a basis C[x]/(D ⊥ + m t+1). This proves the reverse inclusion.

Lemma 4 . 1 .µ

 41 Let B = {x β1 , . . . , x βr t } closed under division, D ∈ H B and Λ = {Λ 1 , . . . , Λ r } be the unique basis of D dual to B with Λ i = |α|≤|βi| µ βi,α d α α! for i = 1, . . . , r. Then Λ 1 = 1 and for i = 2, . . . , r βi,βj +e k ∫ k Λ j . Thus, µ βi,α is a polynomial function of {µ βs,βj +e k : |β s | ≤ |β i |, |β j | < |β s |, k = 1, . . . , n} for i = 1, . . . , r, |α| < |β i |.

  . (1) ⇒ (2) Assume that D ∈ H r and that B t is a basis of C[x]/(D ⊥ + m t+1 ). Let Λ t = {Λ 1 , . . . , Λ rt } be a basis of D t dual to B t with r t = |B t |.Then, for j = r t-1 + 1, . . . , r t , Λ j ∈ D t is such that

  t-1 ) is irreducible (in fact linear) of fixed dimension for D t-1 ∈ U t-1 ⊂ H reg ht-1 , there is a unique irreducible component W t,reg of W t such that π 1 (W t,reg ) = H reg ht-1 (see eg. [29][Theorem 1.26]). We define H reg ht = π 2 (W t,reg ). It is an irreducible component of H ht , since otherwise W t,reg = π -1 2 (H reg ht ) would not be a component of W t but strictly included in one of the irreducible components of W t . Definition 3.4. Let π t : H ht → H ht-1 , D → D ∩ C[d] t-1 be the projection in degree t -1. We define by induction on t, H reg h0 = { 1 } and H reg

	ht	is the irreducible component π -1 t (H reg ht-1 ) of H

ht for t = 1, . . . , δ.

  δ}. By Theorem 3.3, H B is the set of linear spaces D ∈ H r such that D t = D ∩ C[d] t satisfy Equations (8) and (9). It is the open subset of D ∈ H h such that ∆ Bt (D t ) = 0 for t = 1, . . . , δ, where ∆ Bt := ∆ β1,...,βr t denotes the Plücker coordinate for G rt (C[d] t ) corresponding to the monomials in B t .

  to the monomials in B [t] . Let U t := {D t ∈ H reg ht : D t is regular for B t }. Then U t is either an open dense subset of the irreducible variety H reg ht or empty if ∆ B [t] (D [t] ) = 0 for all D ∈ H reg ht . We say that B is a regular basis if U t = H reg ht (or U t = ∅) for t = 1, . . . , δ. We denote by γ [t] = dim G ht (ker H t (D t-1 )) for D t-1 ∈ U t-1 and γ =

	δ t=0 γ [t] .
	If the basis B is regular and closed under division, then H reg h

  r,|α|<r parametrizes a dense subset of H reg h . Proof. Let us define, by induction on t, parameters µ t with |µ t | = t i=1 γ [i]

Table 1 .

 1 .[START_REF] Brianc ¸on | Description de Hilb n C{x, y}[END_REF]. for the instances cmbs1, cmbs2 and KSS. Size of required matrices and parameters for deflation.

	System	r/n	IM	SC #µ	OS
	cmbs1 cmbs2 mth191 decker2 Ojika2 Ojika3 KSS Capr. Cyclic-9 4/9 104 × 33 36 40 11/3 27 × 23 75 74 108 × 77 8/3 21 × 17 21 33 45 × 36 4/3 10 × 9 3 9 15 × 12 4/2 5 × 5 4 8 12 × 10 2/3 6 × 5 0 2 6 × 5 4/3 12 × 9 15 14 27 × 17 16/5 155 × 65 510 362 590 × 367 4/4 22 × 13 6 15 22 × 19 72 × 49

The convergence is seen up to machine error. If we increase the accuracy to 150 digits the rate remains quadratic for 7 iterations: . . . 3.55 • 10 -35 , 6.78 • 10 -70 , 4.15 • 10 -140 , 5.1 • 10 -281 .

Acknowledgments. This research was partly supported by the H2020-MSCA-ITN projects GRAPES (GA 860843) and POEMA (GA 813211) and the NSF grant CCF-1813340.

We shall first compute a primal basis using Algorithm 1. In the first iteration we compute the 3 × 3 matrix K 1 = K 1 (ξ 0 ). The elements in the kernel of this matrix consists of elements of the form Λ = ν 1 1

The singular values of K 1 (ξ 0 ) are (4.1421, 0.0064, 0.0012), which implies a two-dimensional kernel, since two of them are below threshold ε. The (normalized) elements in the kernel are Λ2 = d 1 -0.00117d 2 and Λ3 = d 3 -0.00235d 2 . Note that d 2 was not chosen as a leading term. This is due to pivoting used in the numeric process, in order to avoid leading terms with coefficients below the tolerance ε. The resulting primal basis B 1 = {1, x 1 , x 3 } turns out to be closed under derivation.

Similarly, in degree 2 we compute one element Λ4 =

In the next step, we have ker K 3 = {0}, since the minimum singular value is σ min = 0.21549, therefore we stop the process, since the computed dual is approximately complete (cf. Definition 2.7). We derive that the approximate multiple point has multiplicity r = 4 and one primal basis is

The parametric form of a basis of D 1 is ker

Here we incorporated [START_REF] Emiris | The DMM bound: multivariate (aggregate) separation bounds[END_REF], thus fixing some of the parameters according to primal monomials x 1 and x 3 .

The parametric form of the matrix K 2 (ξ, µ) of the integration method at degree 2 is

, where the columns correspond to the parameters in the expansion (5):

Setting Λ 4 (x 1 x 3 ) = 1 and Λ 4 (x 1 ) = Λ 4 (x 3 ) = Λ 4 (1) = 0, we obtain ν 1 1 = ν 3 1 = 0 and ν 3 2 = 1. The dual element of order 2 has the parametric form

). Overall 8 parameters are used in the representation of D 2 .

The highlighted entries of K 2 (ξ, µ) form the non-singular matrix H 2 in Definition 4.5, therefore D 2 is regular for B (cf. Definition 4.3). We obtain the polynomial parameterization µ 4,3 = µ 2,1 µ 4,2 + µ 3,1 , µ 4,4 = 1, µ 4,5 = µ 2,1 + µ 3,1 µ 4,6 with the free parameters μ = (µ 2,1 , µ 3,1 , µ 4,1 , µ 4,2 , µ 4,6 ). There is no denominator since det H 2 = 1. We now setup the numerical scheme. The overdetermined and deflated system F (x, µ) consists of 15 equations:

We now consider J F (ξ 0 , µ 0 ). This Jacobian is of full rank, and we can obtain a maximal minor by removing Λ 1 (f 2 ), Λ 1 (f 3 ), Λ 2 (f 3 ) and Λ 3 (f 3 ) from F . We obtain the square 11 × 11 system denoted by F 0 .

The initial point of the Newton iterations is ξ 0 = (0.002, 1.003, 0.004) and the approximation of the variables µ i,j provided by the numerical integration method: µ 0 = (-0.00117,-0.00235,5.9•10 -6 ,-0.00002,-0.00235,1.0, -0.00117, -0.00002) .

We now use Theorem 6.1 to certify the convergence to a singular system. We can compute for (ξ 0 , µ 0 ) the value β ≈ 0.01302. Moreover, σ min (K δ+1 (ξ 0 , µ 0 )) = 0.21549 and the minimum singular value of the highlighted submatrix of K 2 (ξ 0 , µ 0 ) is equal to one. Therefore β is at least one order of magnitude less than