Angelos Mantzaflaris
email: angelos.mantzaflaris@inria.fr

An overview of Geometry plus Simulation Modules

Keywords: C++, B-splines, NURBS, isogeometric analysis, geometric design

We give an overview of the open-source library "G+Smo". G+Smo is a C++ library that brings together mathematical tools for geometric design and numerical simulation. It implements the relatively new paradigm of isogeometric analysis, which suggests the use of a unified framework in the design and analysis pipeline. G+Smo is an objectoriented, cross-platform, fully templated library and follows the generic programming principle, with a focus on both efficiency and ease of use. The library aims at providing access to high quality, open-source software to the community of numerical simulation and beyond.

Introducing isogeometric analysis

Isogeometric analysis suggests the use of common spline basis for geometric modeling and finite element analysis of a given model. Engineering process begins with designers encapsulating their concepts on the computer using Computer Aided Design (CAD) software. With the finite element analysis (FEM) technology, the CAD model must be first meshed before any numerical simulation can be performed on it. The mesh generation results in a piecewise linear approximation of the original model, that both increases the data volume and introduces unwanted approximation errors in the geometric description. Moreover, the mesh generation phase involves tedious manual work for producing finite element meshes that are suitable for analysis. Last but not least, the results of the analysis should be projected back to the original CAD model, for allowing the designer to update it accordingly. The data transfer between models suitable for design (CAD) and analysis (FEM) poses a severe problem in industry today.

To address these issues Hughes and co-workers introduced in [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF] an analysis framework which is based on NURBS (Non-Uniform Rational B-Splines), which is standard technology employed in CAD systems. They propose to keep the exact CAD geometry by NURBS surfaces, and apply numerical simulation using the same NURBS as basis functions for the analysis [START_REF] Cottrell | Isogeometric Analysis: Toward Integration of CAD and FEA[END_REF].

General concept of the library

The design and software development philosophy of G+Smo was triggered by the observation that simply augmenting existing libraries for geometric design or numerical analysis will not succeed, since the implementation of isogeometric methods requires a fundamentally new approach to unlock the full set of benefits of the isogeometric paradigm.

The goal of G+Smo is to realize the seamless integration of Finite Element Analysis (FEA) and Computer-aided design (CAD) with open-source code from and to the isogeometric analysis community.

Three general guidelines have been set for the development process. Firstly, we promote both efficiency and ease of use; secondly, we ensure code quality and cross-platform compatibility and, thirdly, we always explore new strategies better suited for isogeometric analysis before adopting FEA practices.

G+Smo is an object-oriented, templated C++ library, that implements a generic concept for IgA, based on abstract classes for geometry map, discretization basis, assemblers, solvers etc. It makes use of object polymorphism and inheritance techniques in order to support a variety of different discretization bases, namely B-spline, Bernstein, NURBS bases, hierarchical and truncated hierarchical B-spline bases of arbitrary polynomial order. The library was first announced in the extended abstract [START_REF] Juettler | Geometry + simulation modules: Implementing isogeometric analysis[END_REF], a more detailed description is contained in [START_REF] Langer | Multipatch discontinuous Galerkin isogeometric analysis[END_REF].

Open-source license

The library is licensed under the Mozilla Public License v2.0 (MPL). The MPL is a simple copy-left license. The MPL's "file-level" copy-left is designed to encourage contributors to share their modifications with the library, while still allowing them to combine the library with code under other licenses (open or proprietary) with minimal restrictions. The latest revision can be downloaded from GitHub: https://github.com/gismo.

Modules

The library is partitioned into smaller entities, called modules. Examples of available modules include the dimension-independent NURBS module, the data fitting and solid segmentation module, the PDE discretization module and the adaptive spline module, based on hierarchical splines of arbitrary dimension and polynomial degree.

At the present point, the B-spline, Bernstein, NURBS bases, and also the truncated and non-truncated hierarchical B-spline bases (cf. [START_REF] Giannelli | THBsplines: An effective mathematical technology for adaptive refinement in geometric design and isogeometric analysis[END_REF]) are in a stable and functional state, for general dimension. From the simulation side, isogeometric simulation algorithms with the Galerkin approach, notably over nonconforming multi-patch physical domains are implemented and in a stable state. A module treating a large class of linear and non-linear elasticity problems was also released recently [START_REF] Shamanskiy | Isogeometric simulation of thermal expansion for twin screw compressors[END_REF]. A module dedicated to isogeometric multigrid solvers on multipatch domains is also available [START_REF] Takacs | Fast multigrid solvers for conforming and non-conforming multi-patch isogeometric analysis[END_REF].

Another interesting module is the assembler module based on expression templates, which are a powerful C++ technique for the development of user-friendly and versatile code. They have the advantage of bringing the implementation of a numerical method much closer to the mathematical notation. This allows the easy formation of isogeometric Galerkin matrices using a high-level language, without sacrificing efficiency.

The library takes advantage of modern, efficient move semantics, which are available in recent C++ versions, for more efficiency and cleaner code.

Code management and documentation

In order to deliver a successful library, quality of the code needs to be assured. To this end, the CMAKE cross-platform compilation system is used in the library, which allows for seamless installation of the library on different platforms. Moreover Jenkins and a CDASH testing servers are in place for executing regular nightly builds, memory checks and unit-tests on different platforms (Windows, Linux and MacOSX). Documentation is done mostly in the form of in-source text, using Doxygen as well as in the Wiki pages http://gs.jku.at/gismo. The documentation is currently available online in the form of HTML pages to all the users, and is updated regularly, as we actively continue adding new material, see https://gismo.github.io.

Third-party dependencies

Apart from the C++ Standard Library, we use open-source third party software for common tasks, for example the Eigen (http://eigen.tuxfamily.org) linear algebra library for (sparse) matrices and linear algebra computations, as well as tools for argument parsing, input and output of XML formatted files, the openNurbs library for reading Rhino's 3DM commercial CAD format (https: //www.rhino3d.com/opennurbs), Trilinos (https://trilinos.github.io) for high performance solvers etc.

Plugins and extensions

Concerning plugins, visualization and user interfaces, we developed a plug-in of G+Smo for the Axl geometric modeling platform ( [START_REF] Christoforou | Axl, a geometric modeler for semi-algebraic shapes[END_REF], http://axl.inria.fr. This provides an interactive GUI (graphical user interface) that allows the user to manipulate multi-patch geometries by their control points and work with complicated CAD models. Nevertheless, our primary approach to visualization is based on the well-established Paraview software (https://www.paraview.org/)and has been enhanced, e.g., with trimmed surface visualization, gradient fields and absolute error plots, etc. Another connection has been setup with Siemens' Parasolid geometric CAD kernel (https://www.plm.automation.siemens.com), which is widely used in industry. Currently, this allows us to exchange data and to use the advanced modeling utilities available in Parasolid from within G+Smo.

Acknowledgement G+Smo is jointly developed by several contributors at the Johannes Kepler University, at the RICAM Institute of the Austrian Academy of Sciences (in the frame of the Austrian Science Fund NFN project S117) at INRIA and at the Department of Applied Mathematics of TU Delft (The Netherlands). More contributions have been made by developers from other institutions. The full list of contributors is available at https://github.com/gismo/gismo/wiki/ About--G-Smo.


Angelos Mantzaflaris

Bernard Mourrain

Singular Zeros of Polynomial Systems

come

Introduction

A main challenge in algebraic and geometric computing is singular point identification and treatment. Such problems naturally occur when computing the topology of implicit curves or surfaces [START_REF] Alberti | Topology and arrangement computation of semi-algebraic planar curves[END_REF], the intersection of parametric surfaces in geometric modeling. When algebraic representations are used, this reduces to solving polynomial systems. Several approaches are available: algebraic techniques such as Gröbner bases or border bases, resultants, subdivision algorithms [START_REF] Mantzaflaris | Continued fraction expansion of real roots of polynomial systems[END_REF], [START_REF] Mourrain | Subdivision methods for solving polynomial equations[END_REF], homotopies, and so on. At the end of the day, a numerical approximation or a box of isolation is usually computed to identify every real root of the polynomial system. But we often need to improve the numerical approximation of the roots. Numerical methods such as Newton's iteration can be used to improve the quality of the approximation, provided that we have a simple root. In the presence of a multiple root, the difficulties are significantly increasing. The numerical approximation can be of very bad quality, and the methods used to compute this approximation are converging slowly (or not converging). The situation in practical problems, as encountered in CAGD for instance, is even worse, since the coefficients of the input equations are known with some incertitude. Computing multiple roots and root multiplicities of approximate polynomial systems is an ill-posed problem, since changing slightly the coefficients may transform a multiple root into a cluster of simple roots (or even make it disappear).

For instance Newton's method converges only linearly to such a point, if it converges at all [START_REF] Gilbert | Newton's method for multiple roots[END_REF]. Also, certification tests for the existence of roots on a domain do not directly treat these cases. On the other hand, computing the local multiplicity structure around a singularity breaks down to stable linear algebra methods, which can be run approximately. One can use this local structure to deflate the root, and thus restore super-linear convergence of Newton iteration, or use standard verification techniques to certify a singular root of the original system. In case of inexact coefficients, known up to a certain tolerance, an exact singular root no longer exists. Nevertheless, a well chosen symbolic perturbation, combined with deflation, allows the certification of a nearby system, within a controlled neighborhood of the original one, which attains a single singular point.

The numerical treatment of singular zeroes is a difficult task, mainly because of the ill-posedness of the problem. The following strategy can however be adopted. Find a perturbation of the input system such that the root is a deformation of an exact multiple root. Certainly, there is not a single multiple system, if the input data is approximate. But using the knowledge of the dual structure and interval arithmetic, our method aims at providing a controlled deformation that is compatible with the input.

In this way, we identify the multiplicity structure and we are able to setup deflation techniques which restore the quadratic convergence of the Newton system. The certification of the multiple root is also possible on the symbolically perturbed system by applying a fixed point theorem, based e.g. on interval arithmetic [START_REF] Rump | Verified error bounds for multiple roots of systems of nonlinear equations[END_REF] or α-theorems ( [START_REF] Giusti | On location and approximation of clusters of zeros: Case of embedding dimension one[END_REF] and references therein).

This approach has already been explored in the past. The first algebraic work on the analysis of singular points may be due to F.S. Macaulay [START_REF] Macaulay | The algebraic theory of modular systems[END_REF], who introduced the terminology of "inverse system". His so-called dialytic method has been exploited in [START_REF] Leykin | Newton's method with deflation for isolated singularities of polynomial systems[END_REF][START_REF] Leykin | Higher-order deflation for polynomial systems with isolated singular solutions[END_REF][START_REF] Dayton | Computing the multiplicity structure in solving polynomial systems[END_REF] to construct the inverse system of a multiple point.

Another construction of inverse systems is described e.g. in [START_REF] Marinari | Gröbner duality and multiplicities in polynomial system solving[END_REF], reducing the size of the intermediate linear systems (and exploited in [START_REF] Stetter | Analysis of zero clusters in multivariate polynomial systems[END_REF]).

In [START_REF] Mourrain | Isolated points, duality and residues[END_REF], another approach to construct the dual basis at the singular point which is based on an integration strategy, has been proposed.

Regarding deflation techniques, in [START_REF] Ojika | Deflation algorithm for the multiple roots of a system of nonlinear equations[END_REF], by applying a triangulation preprocessing step on the Jacobian matrix at the approximate root, minors of the Jacobian matrix are added to the system to reduce the multiplicity.

In [START_REF] Lecerf | Quadratic newton iteration for systems with multiplicity[END_REF], a representation of the ideal in a triangular form in a good position and derivations with respect to the leading variables are used to iteratively reduce the multiplicity. This process is applied for p-adic lifting with exact computation.

In [START_REF] Leykin | Newton's method with deflation for isolated singularities of polynomial systems[END_REF][START_REF] Leykin | Higher-order deflation for polynomial systems with isolated singular solutions[END_REF], instead of triangulating the Jacobian matrix, the number of variables is doubled and new equations are introduced, which are linear in the new variables. They describe the kernel of the Jacobian matrix at the multiple root. The process is iterative, yet for some practical applications, the root may already be deflated with a few iterations.

In [START_REF] Dayton | Computing the multiplicity structure in solving polynomial systems[END_REF], the deflation method is applied iteratively until the root becomes regular, doubling each time the number of variables.

In [START_REF] Pope | Nearest multivariate system with given root multiplicities[END_REF], a minimization approach is used to reduce the value of the equations and their derivatives at the approximate root, assuming a basis of the inverse system is known.

In [START_REF] Wu | Computing the multiplicity structure from geometric involutive form[END_REF], the inverse system is constructed via Macaulay's method; tables of multiplications are deduced and their eigenvalues are used to improve the approximated root. They show that the convergence is quadratic at the multiple root.

Verification of multiple roots of (approximate) polynomial equations is a difficult task. The approach proposed in [START_REF] Rump | Verified error bounds for multiple roots of systems of nonlinear equations[END_REF] consists of introducing perturbation parameters and to certifying the multiple root of nearby system by using a fixed point theorem, based on interval arithmetic. It applies only to cases where the Jacobian has corank equal to 1.

The goal of this paper is to review different techniques that can be used to handle efficiently the following tasks:

(a) Compute a basis for the dual space and of the local quotient ring at a given (approximate) singular point.

(b) Deflate the system by augmenting it with new equations derived from the dual basis, introducing adequate perturbation terms.

(c) Certify the singular point and its multiplicity structure for the perturbed system checking the contraction property of Newton iteration (e.g. via interval arithmetic).

These tools can be applied to improve the quality of approximation of a multiple isolated solution of a system of (polynomial) equations, but they can also be used to solve geometrical problems, such as for instance computing the number of real branches at a singular point of an algebraic curve. For more details on these applications, we refer to [START_REF] Mantzaflaris | Deflation and Certified Isolation of Singular Zeros of Polynomial Systems[END_REF] and references therein.

Preliminary considerations

In this section we present some definitions together with the main tools that we shall need in the sequel.

We denote by R = K[x], x = (x 1 , . . . , x n ), a polynomial ring over the field K of characteristic zero. The dual ring R * is the space of linear functionals Λ : R → K. It is commonly identified to the space of formal series K[[∂]] where ∂ = (∂ 1 , . . . , ∂ n ) are formal variables. Thus we view dual elements as formal series in differential operators at a point ζ ∈ K n . To specify that we use the point ζ, we also denote these differentials

∂ ζ . When applying Λ(∂ ζ ) ∈ K[[∂ ζ ]] to a polynomial g(x) ∈ R we will denote by Λ ζ [g] = Λ ζ g = Λ(∂ ζ )[g(x)] the operation Λ ζ [g] = α∈N n λ α α 1 ! • • • α n ! • d |α| g dx α1 1 • • • dx αn n (ζ), (1)
for

Λ(∂ ζ ) = λ α 1 α! ∂ α ζ ∈ K[[∂ ζ ]]. Extending this definition to an ordered set D = (Λ 1 , . . . , Λ µ ) ∈ K[[∂]] µ , we shall denote D ζ [g] = (Λ ζ 1 g, . . . , Λ ζ µ g).
In some cases, it is convenient to use normalized differentials instead of ∂: for any

α ∈ N n , we denote d α ζ = 1 α! ∂ α ζ .
In particular, with the use of this notation we recover the nice property that, if ζ = 0, we have d α 0 x β = 1 if α = β and 0 otherwise.

More generally, (d α ζ ) α∈N n is the dual basis of ((x -ζ) α ) α∈N n , i.e., a nonzero root implies a linear transformation of the variables, so that the root is translated to (0, 0). Example 1. Consider the integral of a polynomial function g ∈ R over the unit hypercube. Since this is a linear map, it may be expressed in terms of differentials, i.e.:

g -→ [0,1] n g(x) dx 1 • • • dx n = α∈sup(g) c α d α [g] , where d α [g] = 1 α! ∂ |α| g ∂x α (
0) and sup(g) stands for the support of g. Indeed, it can be verified using simple calculations that the (unique) coefficients are given by c α = n i=1

1 α i + 1 .
For Λ ∈ R * and p ∈ R, let us define the operation p • Λ : q → Λ(p q). We check that

(x i -ζ i ) • ∂ α ζ = d d∂ i,ζ (∂ α ζ ), (2)
and R * obtains the structure of an R-module. This property shall be useful in the sequel.

Isolated points and differentials

Let I = f 1 , . . . , f s be an ideal of R and ζ ∈ K n a root of the polynomial system f = (f 1 , . . . , f s ). We call ζ an isolated zero of V (I) if, in a primary decomposition of I, the radical of one of the primary components is the maximal ideal

m ζ = x 1 -ζ 1 , . . . , x n -ζ n defining ζ and no other primary component is contained in m ζ .
Suppose that ζ is an isolated root of f , then a minimal primary decomposition of

I = Q prim.⊃I Q contains a primary component Q ζ such that Q ζ = m ζ and √ Q ⊂ m ζ for the other primary components Q associated to I [2].
As

Q ζ = m ζ , it follows that R/Q ζ is a finite dimensional vector space. The multiplicity µ ζ of ζ is defined as the dimension of R/Q ζ .
A point of multiplicity one is called regular point, or simple root, otherwise we say that ζ is a singular isolated point, or multiple root of f . In the latter case we have

J f (ζ) = 0. Example 2. Consider the ideal I = x 1 -x 2 + x 2 1 , x 1 -x 2 + x 2 2
, and the root ζ = (0, 0). Then a minimal primary decomposition of I is

I = x 3 2 , x 1 -x 2 + x 2 2 ∩ -2 + x 2 , 2 + x 1 .
Among the two factors we find the maximal ideal of ζ given by the radical ideal

x 3 2 , x 1 -x 2 + x 2 2 = x 1 , x 2 .
We can now define the dual space of an ideal. The dual space is known to be isomorphic to the quotient R/I. Consider now the orthogonal of Q ζ , i.e. the subspace D ζ of elements of R * that vanish on members of Q ζ , namely

Q ⊥ ζ = D ζ = {Λ ∈ R * : Λ ζ [p] = 0, ∀p ∈ Q ζ }.
The following is an essential property that allows extraction of the local structure D ζ directly from the "global" ideal I = f , notably by matrix methods that will be outlined in Section 3.

Proposition 1 ([18, Th. 8]). For any isolated point ζ ∈ K of f , we have

I ⊥ ∩ K[∂ ζ ] = D ζ .
In other words, we can identify D ζ = Q ⊥ ζ with the space of polynomial differential operators that vanish at ζ on every element of I. Also note that

D ⊥ ζ = Q ζ .
The space D ζ has dimension µ ζ , the multiplicity at ζ. As the variables (x iζ i ) act on R * as derivations (see (2)), D ζ is a space of differential polynomials in ∂ ζ , which is stable under derivation. This property will be used explicitly

in constructing D ζ (Sect. 3). Definition 2 The nilindex of Q ζ is the maximal integer N ∈ N such that m N ζ ⊂ Q ζ .
It is directly seen that the maximal order of elements in D ζ is equal to N , also known as the depth of the space.

Quotient ring and dual structure

In this section we explore the relation between the dual ring and the quotient It is convenient in terms of notation to make the assumption ζ = 0. This poses no constraint, since it implies only a linear change of coordinates.

R/Q ζ ,
Let suppD 0 be the set of exponents of monomials appearing in D 0 , with a non-zero coefficient. These are of degree at most N , the nilindex of Q 0 . Since

(∀Λ ∈ D 0 , Λ 0 [p] = 0) iff p ∈ D ⊥ 0 = Q 0 ,
we derive that supp D 0 = {γ : x γ / ∈ Q 0 }. In particular, we can find a basis of R/Q 0 between the monomials {x γ : γ ∈ supp D 0 }. This is a finite set of monomials, since their degree is bounded by the nilindex of Q 0 . Now let Let x γj , j = 1, . . . , s be an enumeration of these monomials. It is clear that these are finitely many, since Q 0 is zero-dimensional. Given a monomial basis B = (x βi ) i=1,...,µ of R/Q 0 and, for all monomials x γj / ∈ Q 0 , the expression (normal form)

x γj = µ i=1 λ ij x βi (mod Q 0 ) (3)
of x γj in the basis B, then the dual elements [START_REF] Mourrain | Isolated points, duality and residues[END_REF]Prop. 13]

Λ i (d) = d βi + s-µ j=1 λ ij d γj , (4)
for i = 1, . . . , µ form a basis of D ζ . We give a proof of this fact in the following lemma.

Lemma 1. The set of elements D = (Λ i ) i=1,...,µ defined in ( 4) is a basis of D ζ and the normal form of any g(x) ∈ R with respect to the monomial basis B = (x βi ) i=1,...,µ is

NF(g) = µ i=1 Λ ζ i [g] x βi . (5)
Proof. First note that the elements of D are linearly independent, since d βi appears only in Λ i (d). Now, by construction,

µ i=1 Λ ζ i [x α ]x βi = NF(x α ) , for all x α / ∈ Q ζ , e.g. NF(x βi ) = x βi . Also, for x α ∈ Q ζ , ∀i, Λ ζ i (x α ) = 0, since α / ∈ supp D.
Thus the elements of D compute NF(•) on all monomials of R, and (5) follows by linearity. We deduce that D generates the dual, as in Definition 1.

It becomes clear that with the knowledge of the dual basis at ζ, we are able to compute any g ∈ R modulo Q ζ by applying the basis elements to the monomials of g (formal derivation plus evaluation at ζ). This lemma also shows an isomorphism between the dual D ζ and the quotient ring R/Q ζ , since it implies a one-to-one mapping between the primal and dual basis.

Example 3. Consider f (x, y) = x 4 + 2x 2 y 2 + y 4 + 3x 2 y -y 3 and g(x, y) = 18xy 2 -6x 3 . The common zero ζ = (0, 0) yields the local dual space D = (1, d x , d y , d x 2 , d x d y , d y 2 , d x 3 + 1 3 d x d y 2 , d x 2 d y + 3 d y 3 , d x 4 + 1 3 d x 2 d y 2 + d y 4 + 8 3 d 3 y ),
therefore ζ is a singular zero with multiplicity m = 9. The primal counterpart is B = (1, x, y, x 2 , xy, y 2 , x 3 , x 2 y, x 4 ). The relation between D and B is revealed in the following construction:

              1 x y x 2 xy y 2 x 3 x 2 y x 4 xy 2 y 3 x 2 y 2 y 4 1 1 0 0 0 0 d x 1 0 0 0 0 d y 1 0 0 0 0 d 2 x 1 0 0 0 0 d x d y 1 0 0 0 0 d 2 y 1 0 0 0 0 d 3 x 1 1/3 0 0 0 d 2 x d y 1 0 3 0 0 d 4 x 1 0 8/3 1/3 1               . (6)
The dual monomial of every row couples with a primal monomial in the corresponding column. From the rows we read the coefficients basis elements in D.

The leftmost 9 × 9 block is the identity matrix, implying the duality between D and B. Then in the last four columns there are some extra monomials; these do not belong to the basis B, yet appear with a non-zero coefficient in D. These monomials are not in Q ζ , but they can be reduced modulo B: the last four columns yield the normal form of these monomials with respect to B. For example, using column 11 we find

y 3 = 3x 2 y + 8 3 x 4 (mod Q ζ ).
Using the normal form formula (5), we can derive the table of multiplication by x and y in the quotient algebra represented by B. To do this, it suffices to compute NF(y x αi y βi ) and NF(x x αi y βi ), for all monomials x αi y βi ∈ B. This computation can be done by looking up the normal form of each monomial from the rows of matrix [START_REF] Gilbert | Newton's method for multiple roots[END_REF]. The coefficients of these normal forms fill the i-th rows of the matrices

M x =                      
0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1/3 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

                      and M y =                      
0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1/3 0 0 0 0 0 0 0 0 0 3 8/3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1/3 0 0 0 0 0 0 0 0 0

                      .
Computing the normal form of the border monomials of B via (5) also yields the border basis relations and the operators of multiplication in the quotient R/Q ζ (see e.g. [START_REF] Elkadi | Introduction à la résolution des systèmes d'équations algébriques[END_REF] for more properties).

If a graded monomial ordering is fixed and B = (x βi ) i=1,..,µ is the corresponding monomial basis of R/Q 0 , then d βi is the leading term of (4) with respect to the reversed ordering (that is, we reverse the outcome of the comparison of two monomials, keeping equality unchanged) [START_REF] Leykin | Higher-order deflation for polynomial systems with isolated singular solutions[END_REF]Th. 3.1].

Conversely, if we are given a basis D of D ζ whose coefficient matrix in the dual monomials basis

(d α ) α / ∈Q ζ is D ∈ K µ×s , we can compute a basis of R/Q ζ by choosing µ independent columns of D, say those indexed by d βi , i = 1, . . . , µ . If G ∈ K µ×µ is the (invertible) matrix formed by these columns, then D := G -1 D, is D =    β 1 • • • β µ γ 1 • • • γ s-µ Λ 1 1 0 λ 1,1 • • • λ 1,s-µ . . . . . . . . . . . . Λ µ 0 1 λ µ,1 • • • λ µ,s-µ   , (7)
i.e. a basis of the form (4). Note that an arbitrary basis of D does not have the above diagonal form, nor does it directly provide a basis for R/Q ζ . However, a basis of this form has the desired property

Λ i [x βj ] = 1 , if i = j 0 , if i = j ,
for all i = 1, . . . , µ.

For t ∈ N, D t denotes the vector space of polynomials of D of degree ≤ t. The Hilbert function h :

N → N is defined by h(t) = dim(D t ), t ≥ 0, hence h(0) = 1 and h(t) = dim D for t ≥ N . The integer h(1) -1 = corank J f is known as the breadth of D.

Computing local ring structure

The computation of a local basis, given a system and a point, is done essentially by matrix-kernel computations, and consequently it can be carried out numerically, even when the point or even the system is inexact. Throughout the section we suppose f ∈ R m and ζ ∈ K n with f (ζ) = 0.

Several matrix constructions have been proposed that use different conditions to identify the dual space as a null-space. They are based on the stability property of the dual basis:

∀ Λ ∈ D t , d d∂ i Λ ∈ D t-1 , i = 1, . . . , n. (8)
We list existing algorithms that compute dual-space bases:

• As pointed out in (2), an equivalent form of ( 8) is

∀Λ ∈ D t , Λ[g i f i ] = 0 , ∀g i ∈ R ⇐⇒ Λ[x β • f i ] = 0 , ∀β ∈ N n (9)
Macaulay's method [START_REF] Macaulay | The algebraic theory of modular systems[END_REF] uses this equivalent characterization to derive the algorithm that is outlined in Sect. 3.1.

• In [START_REF] Marinari | Gröbner duality and multiplicities in polynomial system solving[END_REF] they exploit (8) by forming the matrix D i of the map

d d∂ i : K[∂] t → K[∂] t-1
for all i = 1, . . . , n and some triangular decomposition of the differential polynomials in terms of differential variables. This approach was used in [START_REF] Stetter | Analysis of zero clusters in multivariate polynomial systems[END_REF] to reduce the row dimension of Macaulay's matrix, but not the column dimension.

• The closedness subspace method of Zeng [START_REF] Zeng | The closedness subspace method for computing the multiplicity structure of a polynomial system[END_REF], uses the same condition to identify a superset of supp D t+1 when a basis of D t is computed, and thus reduces the column dimension of the matrix.

• The integration method in [START_REF] Mourrain | Isolated points, duality and residues[END_REF] "integrates" elements of a basis of D t , and obtains a priori knowledge of the form of elements in degree t + 1 (Sect. 3.2).

All methods are incremental, in the sense that they start by setting D 0 = (1) and continue by computing

D i , i = 1, . . . , N, N + 1. When #D N = #D N +1
then D N is a basis of D, and N is the nilindex of Q.

We shall review two of these approaches to compute a basis for D, and then describe an improvement, that allows simultaneous computation of a monomial basis of the quotient ring, while avoiding redundant computations.

Macaulay's dialytic matrices

This matrix construction is presented in [START_REF] Macaulay | The algebraic theory of modular systems[END_REF]Ch. 4], a modern introduction is contained in [START_REF] Dayton | Computing the multiplicity structure in solving polynomial systems[END_REF], together with an implementation of the method in ApaTools1 .

The idea behind the algorithm is the following: An element of D is of the form

Λ(d) = |α|≤N λ α d α
under the condition: Λ 0 evaluates to 0 at any g ∈ f , that is,

Λ 0 (g) = Λ 0 g i f i = 0 ⇐⇒ Λ 0 (x β f i ) = 0 , for all monomials x β , β ∈ N.
If we apply this condition recursively for |α| ≤ N , we get a vector of coefficients (λ α ) |α|≤N in the (right) kernel of the matrix with rows indexed by constraints Λ 0 [x β f i ] = 0, |β| ≤ N -1. A basis of D N is given by the kernel of this matrix in depth N . The method consists in computing the kernel of these matrices for N = 1, 2, . . .; when N reaches the nilindex of I, For some value of N , this kernel stabilizes and the generating vectors form a basis of D.

Note that the only requirement is to be able to perform derivation of the input equations and evaluation at ζ = 0.

Example 4. Let f 1 = x 1 -x 2 + x 2 1 , f 2 = x 1 -x 2 + x 2 2 .
We also refer the reader to [4, Ex. 2] for a detailed demonstration of this instance. The matrices in order 1 and 2 are:

1 d 1 d 2 f 1 0 1 -1 f 2 0 1 -1 ,         1 d 1 d 2 d 2 1 d 1 d 2 d 2 2 f 1 0 1 -1 1 0 0 f 2 0 1 -1 0 0 1 x 1 f 1 0 0 0 1 -1 0 x 1 f 2 0 0 0 1 -1 0 x 2 f 1 0 0 0 0 1 -1 x 2 f 2 0 0 0 0 1 -1        
.

The kernel of the left matrix gives

D 1 = (1, d 1 + d 2 )
. Expanding up to order two, we get the matrix on the right, and

D 2 = (1, d 1 +d 2 , -d 1 +d 2 1 +d 1 d 2 +d 2 2
). If we expand up to depth 3 we get the same null-space, thus D = D 2 .

Integration method

This method is presented in [START_REF] Mourrain | Isolated points, duality and residues[END_REF]. It is an evolution of Macaulay's method, since the matrices are not indexed by all differentials, but just by elements based on knowledge of the previous step. This performs a computation adapted to the given input and results in smaller matrices.

For

Λ ∈ K[∂],
Λ(∂) = s i=1 n k=1 λ ik k Λ i (∂ 1 , . . . , ∂ k , 0, . . . , 0), ( 10
)
for λ ik ∈ K, and the following two conditions hold:

(i) s i=1 λ ik d d∂ l Λ i (∂) - s i=1 λ il d d∂ k Λ i (∂) = 0, for all 1 ≤ k < l ≤ n . (ii) Λ ζ [f k ] = 0, for k = 1, . . . , m . Condition (i) is equivalent to d d∂ k Λ ∈ D t-1
, for all k. Thus the two conditions express exactly the fact that D must be stable under derivation and its members must vanish on f . This gives the following algorithm to compute the dual basis: Start with D 0 = 1 . Given a basis of D t-1 we generate the ns candidate elements k Λ i-1 (∂ 1 , . . . , ∂ k , 0, . . . , 0). Conditions (i) and (ii) give a linear system with unknowns λ ik . The columns of the corresponding matrix are indexed by the candidate elements. Then, the kernel of this matrix gives a basis of D t , which we use to generate new candidate elements. If for some t we compute a kernel of the same dimension as D t-1 , then we have a basis of D.

Example 5. Consider the instance of Ex. 4,

f 1 = x 1 -x 2 +x 2 1 , f 2 = x 1 -x 2 +x 2 2 . We have f 1 (ζ) = f 2 (ζ) = 0, thus we set D 0 = {1}. Equation (10) gives Λ = λ 1 d 1 + λ 2 d 2 .
Condition (i) induces no constraints and (ii) yields the system

1 -1 1 -1 λ 1 λ 2 = 0 ( 11
)
where the columns are indexed by d 1 , d 2 . We get λ 1 = λ 2 = 1 from the kernel of this matrix, thus

D 1 = {1, d 1 + d 2 }.
For the second step, we compute the elements of D 2 , that must be of the form

Λ = λ 1 d 1 + λ 2 d 2 + λ 3 d 2 1 + λ 4 (d 1 d 2 + d 2 2
). Condition (i) yields λ 3 -λ 4 = 0, and together with (ii) we form the system

  0 0 1 -1 1 -1 1 0 1 -1 0 1      λ 1 . . . λ 4    = 0, (12)
with columns indexed by

d 1 , d 2 , d 2 1 , d 1 d 2 +d 2 2 .
We get two vectors in the kernel, the first yielding again d 1 + d 2 and a second one for

λ 1 = -1, λ 2 = 0, λ 3 = λ 4 = 1, so we deduce that -d 1 + d 2 1 + d 1 d 2 + d 2 2 is a new element of D 2 .
In the third step we have

Λ = λ 1 d 1 + λ 2 d 2 + λ 3 d 2 1 + λ 4 (d 1 d 2 + d 2 2 ) + ( 13
)
λ 5 (d 3 1 -d 2 1 ) + λ 6 (d 3 2 + d 1 d 2 2 + d 2 1 d 2 -d 1 d 2 ), condition (i) leads to λ 3 -λ 4 + (λ 5 -λ 6 )(d 1 + d 2 )
= 0, and together with condition (ii) we arrive at

    0 0 0 0 1 -1 0 0 1 -1 0 0 1 -1 1 0 -1 0 1 -1 0 1 0 0        λ 1 . . . λ 6    = 0, (14)
of size 4 × 6, having two kernel elements that are already in D 2 . We derive that D = D 2 = D 3 and the algorithm terminates. Note that for this example Macaulay's method ends with a matrix of size 12 × 10, instead of 4 × 6 in this approach.

Computing a primal-dual pair

In this section we provide a process that allows simultaneous computation of a basis pair (D, B) of D and R/Q.

Computing a basis of D degree by degree involves duplicated computations. The successive spaces computed are

D 1 ⊂ • • • ⊂ D N = D N +1 . It is more efficient to produce only new elements Λ ∈ D t , independent in D t /D t-1 , at step t.
Also, once a dual basis is computed, one has to transform it to the form (4), in order to identify a basis of R/Q as well. This transformation can be done a posteriori, by finding a sub-matrix of full rank and then performing Gauss-Jordan elimination over this sub-matrix, to reach matrix form [START_REF] Giusti | On location and approximation of clusters of zeros: Case of embedding dimension one[END_REF].

We introduce a condition (iii) extending Th. 3, that addresses these two issues: It allows the computation of a total of µ independent elements throughout execution, and returns a "triangular" basis, e.g. a basis of R/Q is identified.

Lemma 2. Let D t-1 = (Λ 1 , . . . , Λ k ) be a basis of D t-1 , whose coefficient matrix is

   β 1 • • • β k γ 1 • • • γ s-k Λ 1 1 * * * • • • * . . . 0 . . . * . . . . . . . . . Λ k 0 0 1 * • • • *   , (15)
yielding the monomial basis B t-1 = (x βi ) i=1,...,k . Also, let Λ ∈ K[∂] be of the form [START_REF] Krawczyk | Newton-algorithmen zur bestimmung von nullstellen mit fehlerschranken[END_REF], satisfying (i-ii) of Th. 3.

If we impose the additional condition:

(iii) Λ ζ [x βi ] = 0, 1 ≤ i ≤ k,
then the kernel of the matrix implied by (i-iii) is isomorphic to D t /D t-1 . Consequently, it extends D t-1 to a basis of D t .

Proof. Let S be the kernel of the matrix implied by (i-iii), and let Λ ∈ K[∂] be a non-zero functional in S. We have

Λ ∈ D t and Λ ζ [x βi ] = 0 for i = 1, . . . , k. First we show that Λ / ∈ D t-1 . If Λ ∈ D t-1 , then Λ = k i=1 λ i Λ i .
Take for i 0 the minimal i such that λ i = 0. Then Λ ζ [x βi 0 ] = λ i0 , which contradicts condition (iii). Therefore, S ∩ D t-1 = {0}, and S can be naturally embedded in

D t /D t-1 , i.e. dim S ≤ dim D t -dim D t-1 .
It remains to show that dim S is exactly dim D t -dim D t-1 . This is true, since with condition (iii) we added k = dim D t-1 equations, thus we excluded from the initial kernel (equal to D t ) of (i-ii) a subspace of dimension at most k = dim D t-1 , so that dim S ≥ dim D t -dim D t-1 .

We deduce that S ∼ = D t /D t-1 , thus a basis of S extends D t-1 to a basis of D t .

The above condition is easy to realize; it is equivalent to ∀i, d βi / ∈ supp Λ, which implies adding a row (linear constraint) for every i. If we choose the elements of B with a "reversed" total degree ordering (if a monomial compares total-degree "less than" another one, then it compares "bigger than" the same monomial in the reversed order), then in many cases this constraint becomes λ ik = 0 for some i, k. In this case we rather remove the column corresponding to λ ik instead of adding a row. Hence this lemma allows to shrink the kernel (but also the dimension) of the matrix and compute only new dual elements, which are reduced modulo the previous basis. For a detailed size comparison, see 1.

Let us explore our running example, to demonstrate the essence of this improvement.

Example 6. We re-run Ex. 5 using Lem. 2. In the initialization step D 0 = (1) is already in triangular form with respect to B 0 = {1}. For the first step, we demand Λ[1] = 0, thus the matrix is the same as [START_REF] Lecerf | Quadratic newton iteration for systems with multiplicity[END_REF], yielding D 1 = (1, d 1 + d 2 ). We extend B 1 = {1, x 2 }, so that D 1 is triangular with respect to B 1 .

In the second step we remove from 12 the second column, hence we are left with the 3 × 3 system

  0 1 -1 1 1 0 1 0 1     λ 1 λ 3 λ 4   = 0, yielding a single solution -d 1 + d 2 1 + d 1 d 2 + d 2 2 .
We extend B 1 by adding the monomial

x 1 : B 2 = {1, x 2 , x 1 }.
For the final step, we search an element of the form (13) with Λ[x 1 ] = Λ[x 2 ] = 0, and together with (i-ii) we get:

    0 0 1 -1 1 -1 0 0 1 0 -1 0 0 1 0 0        λ 3 . . . λ 6    = 0.
We find an empty kernel, thus we recover the triangular basis D = D 2 , which can be diagonalized to reach the form:

  1 d 2 d 1 d 2 1 d 1 d 2 d 2 2 Λ 1 1 0 0 0 0 0 Λ 2 0 1 0 1 1 1 Λ 3 0 0 1 -1 -1 -1   .
This diagonal basis is dual to the basis B = (1, x 2 , x 1 ) of the quotient ring and also provides a normal form algorithm (Lem. 1) with respect to B. In the final step we generated a 4 × 4 matrix, of smaller size compared to all previous methods.

Another example is treated in figure 1, with the aid of pictures.

This technique for computing B can be applied similarly to other matrix methods, e.g. Macaulay's dialytic method.

If h(t) -h(t -1) > 1, i.e. there is more than one element in step t, then the choice of monomials to add to B is obtained by extracting a non-zero maximal minor from the coefficient matrix in (d α ). In practice, we will look first at the minimum monomials with respect to a fixed term ordering.

Deflation of a singular point

Deflation techniques allow to transform a system of equations defining a singular solution into a new system where the solution corresponds to a simple point. Usually this is done by adding new variables and new equations so that a simple isolated solution of the extended system projects onto the singular solution of the initial system. We will illustrate different types of deflation. x 2

x 1 (e) In degree 4 we get 2d 4 1 -5d1d 2 2 , 2d 4 2 -5 d 2 1 d2, and we choose x 4 1 .

x 2

x 1 (f) In degree 5, we get 2d therefore x 5 1 is a primal monomial.

x 2

x 1 (g) We find no element of degree 6 and the algorithm stops: ζ has multiplicity 11.

x 2

x 1 (h) Different choices at each step lead to different bases; here is another possible result.

The univariate case

In preparation for the multivariate case, we review an approach for treating singularities of univariate polynomials. Let g(x) ∈ K[x] be a polynomial which attains at x = 0 a root of multiplicity µ > 1. The latter is defined as the positive integer µ such that d µ g(0) = 0 whereas g(0

) = dg(0) = • • • = d µ-1 g(0) = 0. Here we denote by d k g(x) = 1 k! d k dx k g(x)
the normalized derivative of order k with respect to x.

We see that D 0 = 1, d, . . . , d µ-1 is the maximal space of differentials which is stable under derivation, that vanish when applied to members of Q 0 , the x -primary component of g at x = 0.

Example 7. Let g(x) = (x -1) 4 , and ζ = 1. First we check that the space of differentials that vanish on the solution include all linear combinations of

D ζ = 1, d ζ , d 2 ζ , d 3 ζ . For instance, we compute d 2 [g] = 6(x -1) 2 | ζ = 0. Now d 4 is not a member of D ζ since d 4 [g] = 1 = 0 does not vanish.
Similarly, for all i ≥ 4, there exists a member of the ideal generated by g which does not evaluate to zero when we apply the differential d i , namely

d i [x i-4 g] = 1.
We conclude that the local dual space is exactly D ζ , and verify that ζ 4-fold zero of f . Consider now the symbolically perturbed equation

f 1 (x, ε) = g(x) + ε 1 + ε 2 x + • • • + ε µ-1 x µ-2 (16)
and apply every basis element of D 0 to arrive at the new system

f (x, ε) = f 1 , df 1 , . . . , d µ-1 f 1
in µ -1 variables. The i-th equation has the form

f i = d i-1 f 1 = d i-1 g + µ-1 k=i k -1 i -1 x k-i ε k ,
i.e linear in ε, the last one being f µ = d µ-1 g(x). This system deflates the root, since the determinant of its Jacobian matrix at (0, 0) is

det J f (0, 0) = d dx f 1 . . . d dx f µ-1 1 0 . . . 0 1 d dx f µ 0 = -df µ (0) = -µd µ g(0) = 0.
Now suppose that ζ * is an approximate zero, close to x = ζ. We can still compute D ζ by evaluating g(x) and the derivatives up to a threshold relative to the error present in ζ * . Then we can form [START_REF] Mantzaflaris | Continued fraction expansion of real roots of polynomial systems[END_REF] and use verification techniques to certify the root. Checking that the Newton operator is contracting shows the existence and unicity of a multiple root in a neighborhood of the input data. We are going to extend this approach, described in [START_REF] Rump | Verified error bounds for multiple roots of systems of nonlinear equations[END_REF], to multi-dimensional isolated multiple roots.

Deflation using the dialytic approach

Let us consider a system of equations f = (f 1 , . . . , f s ),

f k ∈ R[x], s ≥ n,
which has an isolated zero ζ.

If the Jacobian matrix

J f (x) =    ∂ x1 f 1 (x) • • • ∂ xn f 1 (x) . . . . . . ∂ x1 f s (x) • • • ∂ xn f s (x)   
at the point ζ is of (maximal) rank n, then the root ζ is simple. Moreover the iteration

x (n+1) = x (n) + J f (x (n) ) + f (x (n) )
converges to ζ as soon as the initial point x (0) is close enough to ζ [START_REF] Kantorovich | Functional analysis and applied mathematics[END_REF]. When the root is not simple, then the rank of J f (ζ) is r 1 < n and there are n -r 1 linearly independent differentials of order 1 of the form

u 1 ∂ x1 + • • • + u n ∂ xn which satisfy u 1 ∂ x1 f (ζ) + • • • + u n ∂ xn f (ζ) = 0, (17)
or equivalently

J f (ζ)    u 1 . . . u n    = 0.
To fix a solution of this system, we can choose n -r 1 random vectors r j = (r 0,j , r 1,j , . . . , r j,n ) for j = 1, . . . , n -r 1 and consider the equations

u 1 r 1,j + • • • + u n r j,n + r j,0 = 0, j = 1, . . . , n -r 1 . (18)
Extending the initial system of equations f (x) with the new equations ( 17), [START_REF] Mourrain | Isolated points, duality and residues[END_REF], we obtain a new system of equations denoted f [START_REF] Alberti | Topology and arrangement computation of semi-algebraic planar curves[END_REF] (x [START_REF] Alberti | Topology and arrangement computation of semi-algebraic planar curves[END_REF] ) in the variables

x [1] = (x 1 , . . . , x n , u 1 , . . . , u n ).
This system is called a system deflated from f . By construction, if ζ is an isolated root of f (x) = 0 and rank J f (ζ) = r 1 , there is a unique u * satisfying equations ( 17), [START_REF] Mourrain | Isolated points, duality and residues[END_REF]. Thus x *

[1] = (ζ, u * ) is an isolated solution of the system f [1] (x [1] ) = 0.
If the root x * [1] of this system is simple, then Newton iteration applied on f [START_REF] Alberti | Topology and arrangement computation of semi-algebraic planar curves[END_REF] will converge quadratically to x *

[1] for an initial point x (0) [START_REF] Alberti | Topology and arrangement computation of semi-algebraic planar curves[END_REF] in its neighborhood.

If the root is not simple, the deflation can be applied to the system f [START_REF] Alberti | Topology and arrangement computation of semi-algebraic planar curves[END_REF] (x [START_REF] Alberti | Topology and arrangement computation of semi-algebraic planar curves[END_REF] ) = 0 and we obtain a new system f [START_REF] Atiyah | Introduction to Commutative Algebra[END_REF] (x [START_REF] Atiyah | Introduction to Commutative Algebra[END_REF] ) = 0 in 4 n variables.

As shown in the next result, the process can be applied inductively until the root becomes simple: Theorem 4 ( [START_REF] Leykin | Newton's method with deflation for isolated singularities of polynomial systems[END_REF], [START_REF] Dayton | Multiple zeros of nonlinear systems[END_REF]). If ζ is an isolated root of the system f (x) = 0, there exists a number k ∈ N such that f [k] (x [k] ) = 0 has a simple root x *

[k] whose projection on the first n coordinates is ζ.

It is proved in [START_REF] Leykin | Newton's method with deflation for isolated singularities of polynomial systems[END_REF] (or in [START_REF] Dayton | Multiple zeros of nonlinear systems[END_REF]) that the number k of iterations is at most the depth of the multiplicity of f at ζ, that is the maximum degree of a differential polynomial of the inverse system of f at ζ.

Notice that the number of variables of the system

f [k] is n × 2 k . Example 8. Consider the system f 1 (x 1 , x 2 ) = 2x 1 x 2 2 + 5x 4 1 , f 2 (x 1 , x 2 ) = 2x 2 1 x 2 + 5x 4
2 and the singular point ζ = (0, 0). Since J f (ζ) = 0, we apply a first deflation step, i.e. we compute the equations J f (x) u 1 u 2 , and two random linear equations:

g 1 (x [1] ) = f 1 = 2 x 1 x 2 2 + 5 x 1 4 , g 2 (x [1] ) = f 2 = 2 x 1 2 x 2 + 5 x 2 4 g 3 (x [1] ) = 2 x 2 2 + 20 x 1 3 u 1 + 4 x 1 x 2 u 2 g 4 (x [1] ) = 4 x 1 x 2 u 1 + 2 x 1 2 + 20 x 2 3 u 2 g 5 (x [1] ) = 16 u 1 + u 2 -1, g 6 (x [1] ) = 70 u 1 + 77 u 2
The new Jacobian matrix

J g (ζ [1] ) is rank-defect, with ζ [1] = (x, u) = 0, 0, 11 166 , - 5
83 is zero, and the multiplicity has dropped from 11 to 4. Therefore we repeat the procedure for this new system using random equations from the kernel of J (g) (ζ [START_REF] Alberti | Topology and arrangement computation of semi-algebraic planar curves[END_REF] ) (18):

h1(x [2] ) = 2 x1x2 2 + 5 x1 4 , h2(x [2] ) = 2 x1 2 x2 + 5 x2 4 h3(x [2] ) = 2 x2 2 + 20 x1 3 u1 + 4 x1x2u2 h4(x [2] ) = 4 x1x2u1 + 2 x1 2 + 20 x2 3 u2 h5(x [2] ) = 16 u1 + u2 -1, h6(x [2] ) = 70 u1 + 77 u2 h7(x [2] ) = 2 x2 2 + 20 x1 3 v1 + 4 x1x2v2, h8(x [2] ) = 4 x1x2v1 + 2 x1 2 + 20 x2 3 v2 h9(x [2] ) = 60 x1 2 u1 + 4 x2u2 v1 + (4 x2u1 + 4 x1u2) v2 + 2 x2 2 + 20 x1 3 v3 + 4 x1x2v4 h10(x [2] ) = (4 x2u1 + 4 x1u2) v1 + 4 x1u1 + 60 x2 2 u2 v2 + 4 x1x2v3 + 2 x1 2 + 20 x2 3 v4 h11(x [2] ) = 16 v3 + v4, h12(x [2] ) = 70 v3 + 77 v4 h13(x [2] ) = 53 v1 + 12 v2 + 19 v3 + 63 v4 -1, h14(x [2] ) = 40 v1 + 90 v2 + 3 v3 + 49 v4
We obtain a new system which has a regular root at

ζ [2] = (x, u, v) = 0, 0, 11 166 , -5 83 , 3 143 , - 4 429
, 0, 0 , so deflation is achieved in two steps. We see that the number of (equations and) variables increased exponentially, from 2 to 8, and the system is no longer square.

Deflation using the inverse system

We consider again a system of equations f = (f 1 , . . . , f s ),

f k ∈ R[x],
which has an isolated root ζ of multiplicity µ.

In this section, we will also extend the initial system by introducing new variables so that the extended system has a simple isolated root, which projects onto the multiple point ζ. Contrarily to the deflation technique described in Section 4.2, the number of new variables will be directly related with the multiplicity µ of the point.

Let b = ((x -ζ) β1 , . . . , (x -ζ) βµ ) be a basis of R/Q ζ and D = (Λ 1 , . . . , Λ µ ) its dual counterpart, with β 1 = (0, . . . , 0), Λ 1 = 1.
We introduce a new set of equations starting from f , as follows: add to every f k the polynomial

g k = f k + p k , p k = µ i=1 ε i,k (x -ζ) βi where ε k = (ε k,1 , . . . , ε k,µ ) is a new vector of µ variables.
Consider the system

Dg(x, ε) = Λ 1 (∂ x )[g], . . . , Λ µ (∂ x )[g] .
where 1) with ζ replaced by x, i.e. we differentiate g k but we do not evaluate at ζ. This is a system of µs equations, which we shall index Dg(x, ε) = (g 1,1 , . . . , g µ,s ). We have

Λ x [g k ] = Λ i (d x )[g k ] is defined as in (
g ik (x, ε) = Λ x i [f k + p k ] = Λ x i [f k ] + Λ x i [p k ] = Λ x i [f k ] + p i,k (x, ε). Notice that p i,k (ζ, ε) = Λ ζ i [p k ] = ε i,k because D = (Λ 1 , .., Λ µ ) is dual to b.
As the first basis element of D is 1 (the evaluation at the root), the first s equations are g(x, ε) = 0. Note that this system is under-determined, since the number of variables is µ s + n and the number of equations is µs. We shall provide a systematic way to choose n variables and purge them (or better, set them equal to zero).

By linearity of the Jacobian matrix we have

J Dg (x, ε) = J Df (x, ε) + J Dp (x, ε) = [ J Df (x) | 0 ] + [ J x Dp (x, ε) | J ε Dp (x, ε) ], (19)
where J x Dp (x, ε) (resp. J ε Dp (x, ε)) is the Jacobian matrix of Dp with respect to x (resp. ε). By construction the Jacobian matrix J ε Dp (x, ε) of the system p = (Λ x i (p j )) 1≤i,j≤µ is, up to a reordering of the rows and columns, a block diagonal matrix with s blocks of the form

(Λ x i [b j ]) 1≤i,j≤µ .
As D is dual to the basis b, (Λ x i [b j ](ζ, 0)) 1≤i,j≤µ is the identity matrix, the Jacobian J ε Dp (x, ε) evaluated at (ζ, 0) is, up to a reordering of the crows and columns, the identity matrix of dimension µs.

Using decomposition [START_REF] Mourrain | Subdivision methods for solving polynomial equations[END_REF], we easily deduce the following property:

Lemma 3. The µ s×µ s Jacobian matrix J ε Dg (x, ε) is of full rank µ s at (ζ, 0). Another interesting property is the following [START_REF] Mantzaflaris | Deflation and Certified Isolation of Singular Zeros of Polynomial Systems[END_REF]: Lemma 4. The µs × n Jacobian matrices J x Dg (x, ε) and J x Df (x, ε) are of full rank n at (ζ, 0).

We are going to use these properties to construct sub-systems of Dg with a simple root "above" ζ.

The columns of J Dg (x, ε) are indexed by the variables (x, ε), while the rows are indexed by the polynomials g ik . We construct the following systems:

(a) Let Df I be a subsystem of Df s.t. the corresponding n rows of J Df (ζ) are linearly independent (Lemma 4 implies that such a subset exists). We denote by I = {(i 1 , k 1 ), . . . , (i n , k n )} their indices.

(b) Let D g(x, ε) be the square system formed by removing the variables ε k1,i1 , . . . , ε kn,in from Dg(x, ε). Therefore the Jacobian J D g (x, ε) derives from J Dg (x, ε), after purging the columns indexed by ε k1,i1 , . . . , ε kn,in and it's (i j , k j )-th row becomes [∇(Λ x ij gij,kj ) T | 0 ]. A first consequence is the following result, giving a n × n system deduce from the initial system f , with a simple root at ζ: Theorem 5 (Deflation Theorem 1 [START_REF] Mantzaflaris | Deflation and Certified Isolation of Singular Zeros of Polynomial Systems[END_REF]). Let f (x) be a n-variate polynomial system with an µ-fold isolated zero at x = ζ. Then the n × n system Df I (x) = 0, defined in (a), has a simple root at x = ζ. Example 9. In our running example, we expand the rectangular Jacobian matrix of 6 polynomials in (x 1 , x 2 ). Choosing the rows corresponding to f 1 and

(d 1 -d 2 2 -d 1 d 2 -d 2 1 )[f 1 ]
, we find a non-singular minor, hence the resulting system (f 1 , 2x 1 ) has a regular root at ζ = (0, 0).

The deflated system Df I (x) = 0 is a square system in n variables. Contrarily to the deflation approach in [START_REF] Leykin | Newton's method with deflation for isolated singularities of polynomial systems[END_REF][START_REF] Dayton | Computing the multiplicity structure in solving polynomial systems[END_REF], we do not introduce new variables and one step of deflation is provably sufficient. The trade-off is that here we assume that exact dual elements are pointed at by indices I, so as to be able to compute the original multiple root with high accuracy.

On the other hand, when the coefficients are machine numbers, an exact multiple root is unlikely to exist. In the following theorem, we introduce new variables that will allow us later to derive an approximate deflation method. The need to introduce new variables comes from the fact that in practice the exact root is not available, or even worse, the input coefficients contain small error. Therefore, our method shall seek for a slightly perturbed system with an exact multiple zero within a controlled neighborhood of the input, that fits as close as possible to the approximate multiplicity structure of the input system and point.

Theorem 6 (Deflation Theorem 2 [START_REF] Mantzaflaris | Deflation and Certified Isolation of Singular Zeros of Polynomial Systems[END_REF]). Let f (x) be a n-variate polynomial system with a µ-fold isolated root at x = ζ. The square system D g(x, ε) = 0, as defined in (b), has a regular isolated root at (x, ε) = (ζ, 0).

Nevertheless, this deflation does differ from the deflation strategy in [START_REF] Leykin | Newton's method with deflation for isolated singularities of polynomial systems[END_REF][START_REF] Dayton | Computing the multiplicity structure in solving polynomial systems[END_REF]. There, new variables are added that correspond to coefficients of differential elements, thus introducing a perturbation in the dual basis. This is suitable for exact equations, but, in case of perturbed data, the equations do not actually define a true singular point.

Example 10. Consider the system [START_REF] Leykin | Higher-order deflation for polynomial systems with isolated singular solutions[END_REF] of 3 equations in 2 variables

f 1 = x 3 1 + x 1 x 2 2 , f 2 = x 1 x 2 2 + x 3 2 , f 3 = x 2 1 x 2 + x 1 x 2 2
, and the singular point (0, 0) of multiplicity equal to 7.

Suppose that the point is given. Using 3 and 2 we derive the primal-dual pair

D = (1, d 1 , d 2 , d 2 1 , d 1 d 2 , d 2 2 , d 3 2 + d 3 1 + d 2 1 d 2 -d 1 d 2 2 ) ,
where

d 3 2 is underlined to show that it corresponds to x 3 2 in the primal mono- mial basis B = (1, x 1 , x 2 , x 2 1 , x 1 x 2 , x 2 2 , x 3
2 ) The biggest matrix used, in depth 4, was of size 9 × 8, while Macaulay's method terminates with a matrix of size 30 × 15.

To deflate the root, we construct the augmented system Df of 21 equations. The 21 × 2 Jacobian matrix J Df (x) is of rank 2 and a full-rank minor consists of the rows 4 and 5. Therefore, we find the system (d

2 1 [f 1 ], d 1 d 2 [f 1 ]) = (3x 1 , 2x
2 ) which deflates (0, 0). Note that even though both equations of the deflated system derive from f 1 , the functionals used on f 1 are computed using all initial equations.

The perturbed equations are then

g 1 = f 1 + ε 1,1 + ε 1,2 x 1 + ε 1,3 x 2 + +ε 1,4 x 2 2 + ε 1,5 x 3 2 g 2 = f 2 + ε 2,1 + ε 2,2 x 1 + ε 2,3 x 2 + +ε 2,4 x 2 1 + ε 2,5 x 1 x 2 + ε 2,6 x 2 2 + ε 2,7 x 3 2 g 3 = f 3 + ε 3,1 + ε 3,2 x 1 + ε 3,3 x 2 + +ε 3,4 x 2 1 + ε 3,5 x 1 x 2 + ε 3,6 x 2 2 + ε 3,7 x 3 2
and the resulting system Dg has a simple root at (ζ, 0).

Approximate multiple point

In real-life applications it is common to work with approximate inputs. Also, there is the need to (numerically) decide if an (approximate) system possesses a single (real) root in a given domain, notably for use in subdivision-based algorithms, e.g. [START_REF] Mourrain | Subdivision methods for solving polynomial equations[END_REF][START_REF] Mantzaflaris | Continued fraction expansion of real roots of polynomial systems[END_REF].

In the regular case, Smale's α-theory, extending Newton's method, can be used to answer this problem, also partially extended to singular cases in [START_REF] Giusti | On location and approximation of clusters of zeros: Case of embedding dimension one[END_REF], using zero clustering. Another option is to used the following certification test, based on the verification method of Rump [22, Th. 2.1]:

Theorem 7 ([10, 22] Krawczyk-Rump Theorem). Let f ∈ R n , R = K[x]
, be a polynomial system and ζ * ∈ R n a real approximate regular isolated point. Given an interval domain Z ∈ IR n containing ζ * ∈ R n , and an interval matrix M ∈ IR n×n whose i-th column M i satisfies ∇f i (Z) ⊆ M i for i = 1 . . . , n then the following holds: If the interval domain

V f (Z, ζ * ) := -J f (ζ * ) -1 f (ζ * ) + (I -J f (ζ * ) -1 M )Z ( 20
)
is contained in • Z, the interior of Z, then there is a unique ζ ∈ Z with f (ζ) = 0 and the Jacobian matrix J f (ζ) ∈ M is non-singular.

In our implementation we use this latter approach, since it is suitable for inexact data and suits best with the perturbation which is applied. In particular, it coincides with the numerical scheme of [START_REF] Rump | Verified error bounds for multiple roots of systems of nonlinear equations[END_REF] in the univariate case.

In the case of an isolated multiple point of a polynomial system, we applied a deflation to transform it into a regular root of an extended system. The theorem is applied to the system of 6, using an (approximate) structure D. The resulting range of the ε-parameters encloses a system that attains a single multiple root of that structure. Hence the domain for ε-variables reflects the distance of the input system from a precise system with local structure D. Therefore, we obtain a perturbed system in a neighborhood of the input together with a numerically controlled bound on the perturbation coefficients, with a unique multiple root having a prescribed multiplicity.

If the multiple point is known approximately, we use implicitly Taylor's expansion of the polynomials at this approximate point to deduce the dual basis, applying the algorithm of the previous section. The following computation can be applied:

• At each step, the solutions of linear system (10, i-iii) are computed via Singular Value Decomposition. Using a given threshold, we determine the numerical rank and an orthogonal basis of the solutions from the last singular values and the last columns of the right factor of the SVD.

• For the computation of the monomials which define the equations (2, iii) at the next step, we apply QR decomposition on the transpose of the basis to extract a non-zero maximal minor. The monomials indexing this minor are used to determine constraints (10, i-iii). A similar numerical technique is employed in [START_REF] Zeng | The closedness subspace method for computing the multiplicity structure of a polynomial system[END_REF], for Macaulay's method.

Example 11. Let f 1 = x 2 1 x 2 -x 1 x 2 2 , f 2 = x 1 -x 2 2 . The verification method of [START_REF] Rump | Verified error bounds for multiple roots of systems of nonlinear equations[END_REF] applies a linear perturbation to this system, but fails to certify the root ζ = (0, 0).

We consider an approximate point ζ * = (.01, .002) and we compute the approximate multiplicity structure D = (Λ 1 , . . . , Λ 4 ) = (1.0, 1.0d 2 , 1.0d 1 + 1.0d 2 2 , 1.0d 1 d 2 + 1.0d 3 2 ).

The augmented system g(x) = (Λ j [f i ]) = (f 1 , 2.0x 1 x 2 -1.0x 2 2 -1.0x 1 , 2.0x 1 -2.0x 2 , 1.0x 1 -1.0x 2 2 , f 2 , -2.0x 2 , 0., 0.) has a Jacobian matrix:

J g (ζ * ) T = .00 .016 -.99 2.0 1.0 0 0 0 .00 -.02 .016 -2.0 -.004 -2.0 0 0 with a non-zero minor at the third and forth row. Using this information, we apply the following perturbation to the original system: 6 . We apply Th. 7 on g, i.e. we compute V g (Z 1 , ζ * ). For the variable ε 21 the interval is [-.015, .15] ⊆ (-.01, .01), therefore we don't get an answer.

g 1 = x 2 1 x 2 -
We shrink a little Z 1 down to Z 2 = [-.03, .05] × [-.02, .02] × [-.01, .01] 6 and we apply again Th. 7, which results in

V g (Z 2 , (ζ * , 0) ) =                     ⊆ • Z2,
thus we certify that the input equations admit a perturbation of magnitude of .01, so that the perturbed system has a unique exact root within the interval [-.03, .05] × [-.02, .02]. and their derivation with respect to D.

We consider a box Z with center = ζ * and length = .004 at each side. Also, we allow a range E = [-.004, .004] 8 for the variables ε. Applying 7 we get a verified inclusion V g (Z × E, (ζ * , 0)) inside Z × E and we deduce that a unique specialization ε ∈ E "fits" the approximate system f to the multiplicity structure D.

Indeed, one iteration of Newton's method on g(x, ε) gives the approximate point ζ = (1.505535473, .365266196) and corresponding values for ε 0 ∈ E, such that ζ is a 9-digit approximation of the multiple root of the perturbed system g(x, ε 0 ).

In Table 1 we run dual basis computation on the benchmark set of [START_REF] Dayton | Computing the multiplicity structure in solving polynomial systems[END_REF]. Multiplicity, matrix sizes at termination step and computation time is reported. One sees that there is at least an order of gain in the running time using the primal-dual approach.

Definition 1

 1 The dual space of I is the subspace of elements of K[[∂ ζ ]] (formal series of the variables ∂ ζ ), ζ ∈ V (I), that vanish on all the elements of I. It is also called the orthogonal of I and is denoted as I ⊥ .

x 2 x 1 (

 1 a) ζ = (0, 0) nullifies the system, so monomial 1 is inserted in B.

x 2 x 1 (

 1 b) In degree one we have d1 and d2, so x1 and x2 are added to B.

x 2 x 1 (

 1 c) In degree two, we add the dual monomials of d 2 1 , d1d2, and d 2 2 to B.

x 2 x 1 (

 1 d) In degree three, two dual basis elements appear d 3 1 , d 3 2 .

Fig. 1 .

 1 Fig. 1. Discovering a primal-dual basis pair for the root ζ = (0, 0) of the bivariate system {2x1x 2 2 + 5x 4 1 , 2x 2 1 x2 + 5x 4 2 }.

g 8 = 1 .0ε 24

 8124 Now take the box Z 1 = [-.03, .05]×[-.04, .04]×[-.01, .01]

g 1 , 1 3 , g 2 , 1 = f2 + 5 i=1ε 2

 1132152 = f1 + ε 1,1 + ε 1,2 (x 1 -ζ * 1 ) + ε 1,4 (x 1 -ζ * 1 ) ,i (x 1 -ζ * 1 ) i+1

  where Q ζ is the primary component of the isolated point ζ. We show how to extract a basis of this quotient ring from the support of the elements of D ζ and how D ζ can be used to reduce any polynomial modulo Q ζ .

  we denote by k Λ the element Φ ∈ K[∂] with the property d d∂ k Φ(∂) = Λ(∂) and with no constant term with respect to ∂ k . Theorem 3 ( [18, Th. 15]). Let {Λ 1 , Λ 2 , . . . , Λ s } be a basis of D t-1 , that is, the subspace of D of elements of order at most t -1. An element Λ ∈ K[∂] with no constant term lies in D t iff it is of the form:

  x 1 x 2 2 + ε 11 + ε 12 x 2 g 5 = x 1 -x 2 2 + ε 21 + ε 22 x 2 + ε 23 x 1 + ε 24 x 1 x 2 Thus g(x 1 , x 2 , ε 11 , ε 12 , ε 21 , ε 22 , ε 23 , ε 24), computed as before, is a square system with additional equations:g 2 = 1.0x 2 1 -2.0x 1 x 2 + 1.0ε 12 g 3 = 2.0x 1 x 2 -1.0x 2 2 -1.0x 1 g 4 = 2.0x 1 -2.0x 2 g 6 = -2.0x 2 + 1.0ε 22 + 1.0x 1 ε 24 g 7 = 1.0ε 23 + 1.0x 2 ε 24

Table 1 .

 1 Benchmark systems from[START_REF] Dayton | Multiple zeros of nonlinear systems[END_REF], reporting matrix size at the last step of computing a dual local basis, and overall time for primal-dual computation. The computations are done using Maple. Observe that Macaulay's method results in a

	System µ/n	MM'11	Mourrain'97	Macaulay
	cmbs1 11/3 27 × 23 .18s 27 × 33 .95s 105 × 56 1.55s
	cmbs2	8/3 21 × 17 .08s 21 × 24 .39s	60 × 35 .48s
	mth191 4/3 10 × 9 .03s 10 × 12 .07s	30 × 20 .14s
	decker2 4/2	5 × 5 .02s	5 × 8	.05s	20 × 15 .10s
	Ojika2	2/3	6 × 5 .02s	6 × 6	.03s	12 × 10 .04s
	Ojika3	4/3 12 × 9 .07s 12 × 12 .27s	60 × 35 .59s
	KSS	16/5 155 × 65 8.59s 155 × 80 40.41s 630 × 252 70.03s
	Capr.	4/4 22 × 13 .28s 22 × 16 .47s	60 × 35 2.34s
	Cyclic-9 4/9 104 × 33 1.04s 104 × 36 5.47s 495 × 220 31.40s
	DZ1	131/4 700 × 394 14m 700 × 524 26m 4004 × 1365 220m
	DZ2	16/3 43 × 33 .68s 43 × 48 4.38s 360 × 165 25.72s
	DZ3	5/2	6 × 6 .04s 6 × 10 .23s	30 × 21 .79s
	matrix of size n	p -1 + n p -1	×	p + n p	, in contrast to a matrix of size (	n(n -1) 2	µ+

n) × µ(n -1) + 1 for the primal-dual approach.

http://www.neiu.edu/∼zzeng/apatools.htm

Acknowledgements. This research has received funding from the EU's 7 th Framework Programme [FP7/2007-2013], Marie Curie Initial Training Network SAGA, grant n o [PITN-GA-2008-214584].

Experimentation

We have implemented the presented algorithms in Maple. It can compute (approximate) dual bases by means of Macaulay's method as well as the integration method, and it can derive the augmented system defined in Th. 6. Then Krawczyk-Rump's interval method is used to verify the root.

Example 12. Let, as in [START_REF] Lecerf | Quadratic newton iteration for systems with multiplicity[END_REF][START_REF] Leykin | Higher-order deflation for polynomial systems with isolated singular solutions[END_REF],

The point (0, 0, -1) occurs with multiplicity equal to 18, in depth 7. The final matrix size with our method is 54 × 37, while Macaulay's method ends with a 360 × 165 matrix.

If the objective is to deflate as efficiently as possible, then one can go step by step: First compute a basis of D 1 and stop the process. We get the evaluation 1 and 2 first order functionals, which we apply to f 1 . We arrive at

and we check that the Jacobian determinant is 64, thus we have a deflated system only with a partial local structure. The condition number of the Jacobian matrix is also very satisfactory, with a value of around 5.55.

The recent paper [START_REF] Hao | Algorithm 931: An algorithm and software for computing multiplicity structures at zeros of nonlinear systems[END_REF], implementing the dialytic deflation method produces a deflated system of size 75 × 48 for this instance, with a condition number of order 10 6 .

Example 13. Consider the equations (taken from [START_REF] Dayton | Computing the multiplicity structure in solving polynomial systems[END_REF]DZ3]):

and take an approximate system f with those coefficients rounded to 6 digits. A 5-fold zero of f rounded to 6 digits is ζ * = (1.50551, .365278).

Starting with the approximate system and with a tolerance of .001, we compute the basis

having 4 correct digits, with respect to the initial exact system, and the primal counterpart B = (1, x 1 , x 2 1 , x 3 1 , x 4 1 ). We form the deflated system (b), with I = {(3, 1), (5, 1)}, i.e. the 3rd and 5th dual element on f 1 have non-null Jacobian. By adding 8 new variables, the system is perturbed as:


-Size-selection (slow vs. fast life history) -Competitor (presence vs. absence) -Light intensity (high vs. low) -Size-selection (slow vs. fast life history) -Competitor (presence vs. absence) -Light intensity (high vs. low) Methods: Fish traits Results: Community abundances 1) Selection for small size not only changes life history -but behaviour and stoichiometric traits 2) Growth, behavioural and stoichiometric traits covary -Size-selection for small size <-> fast pace of life 3) Individual traits variation in feeding, C:N and excretion led to changes in 3 invertebrate species abundance Take-home message Thank you for listening

  28.85%) PC 1 (28.85%) PC 2 (21.97%) v-test = -2.39, P = 0.016 Results: Fish traits Competitor presence -> Slow pace of life • -Carbon: +C:N & C:P •smaller effect v-test = -2.06, P = 0.038 Competitor presence -> Slow pace of life • -Carbon: +C:N & C:P • -Slow growth due to lower access to resources *** But smaller effect v-test = -2.06, P = 0.038 Competitor presence -> Slow pace of life • -Carbon: +C:N & C:P • -