Deflation and Certified Isolation of Singular Zeros of Polynomial Systems Angelos Mantzaflaris and Bernard Mourrain GALAAD, INRIA Méditerranée BP 93, 06902 Sophia Antipolis, France arXiv:1101.3140v1 [cs.SC] 17 Jan 2011 [FirstName.LastName]@inria.fr ABSTRACT We develop a new symbolic-numeric algorithm for the certification of singular isolated points, using their associated local ring structure and certified numerical computations. An improvement of an existing method to compute inverse systems is presented, which avoids redundant computation and reduces the size of the intermediate linear systems to solve. We derive a one-step deflation technique, from the description of the multiplicity structure in terms of differentials. The deflated system can be used in Newton-based iterative schemes with quadratic convergence. Starting from a polynomial system and a small-enough neighborhood, we obtain a criterion for the existence and uniqueness of a singular root of a given multiplicity structure, applying a well-chosen symbolic perturbation. Standard verification methods, based eg. on interval arithmetic and a fixed point theorem, are employed to certify that there exists a unique perturbed system with a singular root in the domain. Applications to topological degree computation and to the analysis of real branches of an implicit curve illustrate the method. Categories and Subject Descriptors G.1.5 [Mathematics of Computing]: Roots of Nonlinear Equations; I.1.2 [Computing Methodologies]: Symbolic and Algebraic Manipulation—Algebraic algorithms General Terms Algorithms, Theory Keywords root deflation, multiplicity structure, dual space, inverse system, isolated point 1. 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 [1], the intersection of parametric Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. To copy otherwise, to republish, to post on servers or to redistribute to lists, requires prior specific permission and/or a fee. Copyright 200X ACM X-XXXXX-XX-X/XX/XX ...$5.00. 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[10], [13], homotopies, and so on. At the end of the day, a numerical approximation or a box of isolation is usually computed to characterize the (real) roots of the polynomial system. But we often need to improve the numerical approximation of the roots. Numerical methods such that 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 of approximate polynomial systems seems to be 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). To tackle this difficult problem, we adopt the following strategy. We try to find a (small) perturbation of the input system such that the root we compute is an exact multiple root of this perturbed system. 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 eg. on interval arithmetic [16] or α-theorems [17]. Related work. In order to be develop Newton-type methods that converge to multiple roots, deflation techniques which consist in adding new equations in order to reduce the multiplicity have already been considered. In [14], 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 [6], a presentation 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 [7, 8], 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. In [3], this construction in related to the construction of the inverse system. The dialytic method of F.S. Macaulay [9] is revisited for this purpose. These deflation methods are applied iteratively until the root becomes regular, doubling each time the number of variables. More recent algorithms for the construction of inverse systems are described eg. in [11], reducing the size of the intermediate linear systems (and exploited in [18]), or in [12] using an integration method. In [15], 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 [20], the inverse system is constructed via Macaulay’s method; tables of multiplications are deduced and their eigenvalues are used to improve the approximated root. It is proved that the convergence is quadratic when the Jacobian has corank 1 at the multiple root. Verification of multiple roots of (approximate) polynomial equations is a difficult task. The approach proposed in [16] consists in 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 univariate case. In preparation for the multivariate case, we review some techniques used to treat singular zeros of univariate polynomials, and we present our method on a univariate instance. Let g(x) ∈ K[x] be a polynomial which attains at x = 0 a root of multiplicity µ > 1. The latter is defined as the smallest positive integer µ such that dµ g(0) 6= 0 whereas g(0) = dg(0) = · · · = dµ−1 g(0) = 0. Here we denote dk dk g(x) = dx k g(x)/k! the normalized k−th order derivative. We see that D0 = h1, d1 , . . . , dµ−1 i is the maximal space of differentials which is stable under derivation, that vanish when applied to members of Q0 , the hxi−primary component of hgi at x = 0. Consider now the symbolically perturbed equation f1 (x, ε) = g(x) + ε1 + ε2 x + · · · + εµ−2 xµ−2 (1) and apply every  basis element of D0 to arrive to the new system f (x, ε) = f1 , d1 f1 , . . . , dµ−1 f1 in µ − 1 variables. The P k−i+1 i−th equation shall be fi = di−1 f1 = di−1 g+ µ−2 εk , k=i x µ−1 i.e linear in ε, the last one being fµ = d g(x). This system deflates the root, as we see that the determinant of it’s Jacobian matrix at (0, 0) is d f dx 1 det Jf (0, 0) = .. . d f dx µ−1 d f dx µ 1 0 .. 0 . 1 0 = −µdfµ (0, 0) = −µdµ g(0) 6= 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 in ζ ∗ . Then we can form (1) 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 [16], to multi-dimensional isolated multiple roots. Our approach. It consists of the following steps: (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 (eg. via interval arithmetic). In step (a), a dual basis at the singular point is computed by means of linear algebra, based on the integration approach of [12]. We describe an improvement of this method, which yields directly a triangular dual basis with no redundant computation. This method has the advantage to reduce significantly the size of the linear systems to solve at each step, compared to Macaulay’s type methods [9, 7, 8, 3]. In the case of an approximate singular point, errors are introduced in the coefficients of the basis elements. Yet a successful computation is feasible. In particular, the support of the basis elements is revealed by this approximate process. In the deflation step (b), new equations and new variables are introduced in order to arrive to a new polynomial system where the singularity is obviated. The new variables correspond to perturbations of the initial equations along specific polynomials, which form a dual counterpart to the basis of the dual space. One of the deflated systems that we compute from the dual system is a square n × n system with a simple root. This improves the deflation techniques described in [7, 8, 3], which require additional variables and possibly several deflation steps. New variables are introduced only in the case where we want to certify the multiplicity structure. The perturbation techniques that we use extend the approach of [16] to general cases where the co-rank of the Jacobian matrix could be bigger than one. The verification step (c) is mostly a contraction condition, using eg. techniques as in [16]. This step acts on the (approximate) deflated system, since verifying a simple solution of the deflated system induces a certificate of an exact singular point of (a nearby to) the initial system. We are going to detail the different steps in the following sections, starting with notations in Sect. 2, dual basis in Sect. 3, deflation in Sect. 4, and certification in Sect. 5. In the last section, we will show examples and applications to the topology analysis of curves. 2. PRELIMINARIES AND MAIN RESULTS We denote by R = K[x] a polynomial ring over the field K of characteristic zero. Also, the dual ring R∗ is the space of linear functionals Λ : R → K. It is commonly identified as 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 ζ ∈ Kn . 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] = X α∈Nn λα d|α| g (ζ), · α1 n α1 ! · · · αn ! d1 · · · dα n (2) P 1 ∂ α ∈ K[[∂ ]]. Extending this definition for Λ(∂ ζ ) = λα α! ζ ζ 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 α ∈ Nn , α α β 1 we denote dα ζ = α! ∂ ζ . When ζ = 0, we have d0 x = 1 if α = β and 0 otherwise. More generally, (dα ζ )α∈Nn is the dual basis of ((x − ζ)α )α∈Nn . For Λ ∈ R∗ and p ∈ R, let p · Λ : q 7→ Λ(p q). We check that d (xi − ζi ) · ∂ α (∂ α (3) ζ = ζ ). d∂i,ζ This property shall be useful in the sequel. 2.1 Isolated points and differentials Let hf i = hf1 , . . . , fs i be an ideal of R, ζ ∈ K a root of f and mζ = hx1 − ζ1 , . . . , xn − ζn i the maximal ideal at ζ. Suppose that ζ is an isolated root \ of f , then a minimal priQ contains a primary mary decomposition of I = p Q prim.⊃I √ component Qζ such that Qζ = mζ and Q′ 6⊂ mζ for the other p primary components Q′ associated to I [2]. As Qζ = mζ , 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 Jf (ζ) = 0. We can now define the dual space of an isolated point. Definition 2.1. The dual space of I is the subspace of elements of K[[∂ ζ ]] that vanish on all the elements of I. It is also called the orthogonal of I and denoted by 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 lemma is an essential property that allows extraction of the local structure Dζ directly from the “global” ideal I = hf i, notably by matrix methods outlined in Sect. 3. Proposition 2.2 ([12, Th. 8]). For any isolated point ζ ∈ K of f , we have hf 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. The space Dζ is a vector space of polynomials in ∂ ζ dimension µζ , the multiplicity of ζ. As the variables (xi − ζi ) act on R∗ as derivations (see (3)), Dζ is a space of differential polynomials in ∂ ζ , which is stable by derivation. This property will be used explicitly in constructing Dζ (Sec. 3). Definition 2.3. The nilindex of Qζ is the maximal integer N ∈ N s.t. mN ζ 6⊂ Qζ . It is directly seen that the maximal degree of an element of Dζ has to be equal to N , also known as the depth of Dζ . 2.2 Quotient ring and dual structure In this section we explore the relation between the dual ring and the quotient R/Qζ 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ζ . It is convenient in terms of notation to make the assumption ζ = 0. This saves some indices, while it poses no constraint (since it implies a linear change of coordinates), and shall be adopted hereafter and in the next section. Let suppD0 be the set of exponents of monomials appearing in D0 , with a non-zero coefficient. These are of degree at most N , the nilindex of Q0 . Since ∀Λ ∈ D0 , Λ0 [p] = 0 iff p ∈ Q0 , we derive that supp D0 = {α : xα ∈ / Q0 }. In particular, we can find a basis of R/Q0 between the monomials {xα : α ∈ supp D}. This is a finite set of monomials, since their degree is bounded by the nilindex of Q0 . Given a standard basis B = (xβi )i=1,...,µ of R/Q0 and, for all monomials xγ j ∈ / Q0 , j = 1, . . . , s − µ, s = #supp D0 , with xγ j ∈ / B, the expression (normal form) of xγ j = µ X λij xβ i mod Q0 i=1 (4) in the basis B then the dual elements [12, Prop. 13] s−µ Λi = dβi + X λij dγ j , (5) j=1 for i = 1, . . . , µ form a basis of D0 . We give a proof of this fact in the following lemma. Lemma 2.4. The set of elements D = (Λi )i=1,...,µ is a basis of D0 and the normal form of any g(x) ∈ R with respect to the standard basis B = (xβi )i=1,...,µ is NF(g) = µ X Λ0i [g] xβi . (6) i=1 Proof. First note that the elements P of D are linearly independent. Now, by construction, µi=1 Λ0i [xα ] = NF(xα ) for all xα ∈ / Q0 , eg. NF(xβi ) = xβi . Also, for xα ∈ Q0 , 0 α ∀i, Λi (x ) = 0, since α ∈ / supp D. Thus the elements of D compute NF(·) on all monomials of R, and (6) follows by linearity. We deduce that D is a basis of D0 , as in Def. 2.1. Computing the normal form of the border monomials of B via (6) also yields the border basis relations and the operators of multiplication in the quotient R/Q0 (see eg. [5] for more properties). If a graded monomial ordering is fixed and B = (xβi )i=1,..,µ is the corresponding standard basis of R/Q0 , then dβi is the leading term of (5) wrt the opposite ordering [8, Th. 3.1]. Conversely, if we are given a basis D of D0 whose coefficient matrix in the dual monomials basis (dα )α∈Q / 0 is D ∈ Kµ×s , we can compute a basis of R/Q0 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 β1 Λ′1 ′ D = .. . Λ′µ      ··· 1 .. 0 βµ γ1 ··· 0 λ1,1 .. . λµ,1 ··· . 1 ··· γ s−µ  λ1,s−µ  ..  , .  λµ,s−µ (7) i.e. a basis of the form (5). Note that an arbitrary basis of D does not have the above diagonal form, nor does it directly provide a basis for R/Q0 . For t ∈ N, Dt denotes the vector space of polynomials of D of degree ≤ t. The Hilbert function h : N → N is defined by h(t) = dim(Dt ), t ≥ 0, hence h(0) = 1 and h(t) = dim D for t ≥ N . The integer h(1) − 1 = corank Jf is known as the breadth of D. 3.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 ∈ Rm and ζ ∈ Kn 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: ∀ Λ ∈ Dt , d Λ ∈ Dt−1 , d∂i i = 1, . . . , n. (8) We list existing algorithms that compute dual-space bases: • As pointed out in (3), an equivalent form of (8) is: ∀Λ ∈ Dt , Λ[xi fj ] = 0, ∀i, j = 1, . . . , n. Macaulay’s method [9] uses it to derive the algorithm that is outlined in Sect. 3.1. • In [11] they exploit (8) by forming the matrix Di of the d map : K[∂]t → K[∂]t−1 for all i = 1, . . . , n and some d∂i triangular decomposition of the differential polynomials in terms the differential variables. This approach was used in [18] to reduce the row dimension of Macaulay’s matrix, but not the column dimension. The closedness condition is also used in [21] to identify a superset of supp Dt+1 . • The integration method in [12] “integrates” elements of a basis of Dt , 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 D0 = (1) and continue by computing Di , i = 1, . . . , N, N + 1. When #DN = #DN+1 then DN 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 quotient ring basis while avoiding redundant computations. 3.1 Macaulay’s dialytic matrices This matrix construction is presented in [9, Ch. 4], a modern introduction is contained in [3], together with an implementation of the method in ApaTools1 . The idea behind the algorithm is the An eleXfollowing: λα dα under the ment of D is of the form Λ(∂) = |α|≤N condition: Λ0 evaluates to 0 at any g ∈ hf i, i.e. Λ0 (g) = 0 P Λ ( gi fi ) = 0 ⇐⇒ Λ0 (xβ fi ) = 0 for all monomials xβ . 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. Note that the only requirement is to be able to perform derivation of the input equations and evaluation at ζ = 0. Example 3.1. Let f1 = x1 −x2 +x21 , f2 = x1 −x2 +x22 . We also refer the reader to [3, Ex. 2] for a detailed demonstration of Macaulay’s method on the same instance. The matrices in order 1 and 2 are: f1 f2 1  1 0 0 d1 1 1 f1 d2 f2  −1 , x 1 f1 −1 x 1 f2 x 2 f1 x 2 f2 1 0 0  0  0 0 0  d1 1 1 0 0 0 0 d2 −1 −1 0 0 0 0 d21 1 0 1 1 0 0 http://www.neiu.edu/ zzeng/apatools.htm d1 d2 0 0 −1 −1 1 1 d22  0 1   0 .  0  −1  −1 The kernel of the left matrix gives D1 = (1, d1 + d2 ). Expanding up to order two, we get the matrix on the right, and D2 = (1, d1 + d2 , −d1 + d21 + d1 d2 + d22 ). If we expand up to depth 3 we get the same null-space, thus D = D2 . 2 3.2 Integration method This method is presented in [12]. It is an evolution of Macaulay’s method, in the sense that 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. R For Λ ∈ K[∂], we denote by k Λ the element Φ ∈ K[∂] with the property d∂d Φ(∂) = Λ(∂) and with no constant k term wrt ∂k . Theorem 3.2 ( [12, Th. 15] ). Let hΛ1 , Λ2 , . . . , Λs i be a basis of Dt−1 , that is, the subspace of D of elements of order at most t − 1. An element Λ ∈ K[∂] with no constant term lies in Dt iff it is of the form: Λ(∂) = n s X X λik i=1 k=1 R k Λi (∂1 , . . . , ∂k , 0, . . . , 0), (9) for λik ∈ K, and the following two conditions hold: (i) s X s X d d λil Λi (∂) − Λi (∂) = 0, d∂ d∂ l k i=1 i=1 for all 1 ≤ k ≤ l ≤ n . λik (ii) Λζ [fk ] = 0, for k = 1, . . . , n . Condition (i) is equivalent to d∂dk Λ ∈ Dt−1 , for 1 ≤ k ≤ n. Thus the two conditions express exactly the fact that D must be stable by derivation and his members must vanish on hf i. This gives the following algorithm to compute the dual basis: Start with D0 = h1i. Given R a basis of Dt−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 the matrix gives the possible solutions, which we add to Dt−1 to obtain Dt . If for some t there are no further new elements, then D = Dt is a basis of D. Example 3.3. Consider the instance of Ex. 3.1. We have f1 (ζ) = f2 (ζ) = 0, thus we set D0 = {1}. Equation (9) gives Λ = λ1 d1 + λ2 d2 . Condition (i) induces no constraints and (ii) yields the system  1 1 −1 −1  λ1 λ2  =0 (10) where the columns are indexed by d1 , d2 . We get λ1 = λ2 = 1 from the kernel of this matrix, thus D1 = {1, d1 + d2 }. For the second step, we compute the elements of D2 , that must be of the form Λ = λ1 d1 + λ2 d2 + λ3 d21 + λ4 (d1 d2 + d22 ). Condition (i) yields λ3 − λ4 = 0, and together with (ii) we form the system  0  1 1 0 −1 −1 1 1 0   λ1  −1   0   ..  = 0, . 1 λ4 (11) with columns indexed by d1 , d21 , d2 , d1 d2 + d22 . We get two vectors in the kernel, the first yielding again d1 + d2 and a second one for λ1 = −1, λ2 = 0, λ3 = λ4 = 1, so we deduce that −d1 + d21 + d1 d2 + d22 is a new element of D2 . In the third step we have Λ =λ1 d1 + λ2 d2 + λ3 d21 + λ4 (d1 d2 + d22 )+ λ5 (d31 − d21 ) + λ6 (d32 + d1 d22 + d21 d2 − d1 d2 ), condition (i) leads to λ3 − λ4 + (λ5 − λ6 )d1 + (λ5 − λ6 )d2 = 0, and together with condition (ii) we arrive to 0  0  1 1  0 0 −1 −1 0 1 1 0 0 −1 0 1 1 0 −1 0   −1 λ1 0  .  = 0. . 0  .  λ6 0 (12) Since the kernel of this matrix gives elements that are already in D2 , we derive that D = D2 = D3 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. 2 3.3 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 D1 ⊂ · · · ⊂ DN = DN+1 . It is more efficient to compute only new elements Λ ∈ Dt which are independent in Dt /Dt−1 at step t. Also, once dual basis is computed, one has to transform it into the form (5), 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 (7). We introduce a condition (iii) extending Th. 3.2, 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 3.4. Let Dt−1 = (Λ1 , . . . , Λk ) be a basis of Dt−1 , whose coefficient matrix is Λ1 . .. Λk  β1 ··· βk γ1 ··· γs−k 1 ∗ ∗ ∗ . .. ∗ ··· ∗ . .. ∗     0 0 .. 0 . ∗ 1 ···    ,  (13) yielding the standard basis B = (xβi )i=1,...,k . An element Λ ∈ K[∂] is not zero in Dt /Dt−1 iff in addition to (i), (ii) of Th. 3.2 we impose: (iii) Λ[xβ i ] = 0, 1 ≤ i ≤ k . Proof. Let Λ ∈ K[∂] be a non-zero functional satisfying (i), (ii) and (iii). Then Λ ∈ Dt and Λ[xβi ] = 0 for i = Pk 1, . . . , k. If Λ ∈ Dt−1 , then Λ = i=1 λi Λi . Take for i0 the minimal i such that λi 6= 0. Then Λ[xβ i0 ] = λi0 , which is in contradiction with condition (iii). Thus, the nonzero solutions of (i), (ii) and (iii) correspond to the elements which are not zero in Dt /Dt−1 . The above constraint is easy to realize; it is equivalent β to ∀i, dζ i ∈ / supp Λζ , which implies adding a row (linear constraint) for every i. In many cases this constraint is just λik = 0 for some i, k, thus we rather remove the column corresponding to λik instead of adding a row. Either way, this lemma allows to shrink the kernel of the matrix and compute only new dual elements. Let us explore our running example, to demonstrate the essence of this improvement. Example 3.5. We re-run Ex. 3.3 using Lem. 3.4. In the initialization step D0 = (1) is already in triangular form with respect to B0 = {1}. For the first step, we demand Λ[1] = 0, thus the matrix is the same as (10), yielding D1 = (1, d1 + d2 ). We extend B1 = {1, x2 }, so that D1 is triangular with respect to B1 . In the second step we remove from (11) the second column, hence  0  1 1 1 1 0   −1 λ1 0   λ3  = 0, 1 λ4 yielding a single solution −d1 + d21 + d1 d2 + d22 . We extend B1 by adding monomial x1 : B1 = {1, x2 , x1 }. For the final step, we search an element with Λ[x1 ] = Λ[x2 ] = 0 thus (12) loses two columns: 0  1  1 0  0 −1 0 1 1 0 −1 0   −1 λ3 0  .  = 0. . 0  .  λ6 0 We find an empty kernel, thus we recover the triangular basis D = D2 , which is then diagonalized to reach the form: 1 Λ1 1 Λ2  0 Λ3 0 d2 0 1 0 d1 0 0 1 d21 0 1 −1 d1 d2 0 1 −1 d22  0 1 . −1 This diagonal basis is dual to the basis B = (1, x2 , x1 ) of the quotient ring and also provides a normal form algorithm (Lem. 2.4) wrt B. In the final step we generated a 4 × 4 matrix, size smaller compared to all previous methods. 2 This technique for computing B can be applied similarly to other the matrix methods, e.g. Macaulay’s dialytic method. If h(t) − h(t − 1) > 1, ie. there are more than one elements 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 monomials of smallest degree. 3.4 Approximate dual basis In our deflation method, we assume that the multiple point is known approximately and 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. To handle safely the numerical problems which may occur, we utilize the following techniques: • At each step, the solutions of linear system (9, 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 (3.4, 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 (9, i-iii). A similar numerical technique is employed in [21], for Macaulay’s method. 4. DEFLATION OF A SINGULAR POINT We consider a system of equations f = (f1 , . . . , fs ), which has a multiple point at x = ζ. Also, let B = (b1 , . . . , bµ ) be a basis of R/Qζ and D = (Λ1 , . . . , Λµ ) its dual counterpart, with Λ1 = 1. We introduce a new set of equations starting from f , as follows: add for every fi the polynomial gk = fk + pk , pk = Pµ i=1 εi,k bi where εk = (ε1,k , . . . , εµ,k ) is a new vector of µ variables. Consider the system   Dg(x, ε) = Λ1 (∂ x )[g], . . . , Λµ (∂ x )[g] . where Λx [gk ] = Λi (∂ x )[gk ] is defined as in (2) with ζ replaced by x, ie. we differentiate gk but we do not evaluate at ζ. This is a system of µs equations, which we shall index Dg(x, ε) = (g1,1 , . . . , gµ,s ). We have x x x gik (x, ε) = Λx i [fk +pk ] = Λi [fk ]+Λi [pk ] = Λi [fk ]+pik (x, ε). Notice that pi,k (ζ, ε) = Λζi [pk ] = ε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). This way we arrive to a square system Dg(x, ε̃) (we use ε̃ for the remaining µs − n variables) of size µs × µs. We shall prove that this system vanishes on (ζ, 0) and that JDg (ζ, 0) 6= 0. By linearity of the Jacobian matrix we have JDg (x, ε) = JDf (x, ε) + JDp (x, ε) ε x (x, ε) ], (x, ε) | JDp = [ JDf (x) | 0 ] + [ JDp x ε where JDp (x, ε) (resp. JDp (x, ε)) is the Jacobian matrix of Dp with respect to x (resp. ε). ε Lemma 4.1. The Jacobian JDp (x, ε) of the linear system Dp = (p1,1 , . . . , pµ,s ) with pi,k (εk ) = Λx i [pk ](x, εk ) evaluated at (x, ε) = (ζ, 0) is the identity matrix in dimension µs. Proof of Lemma 4.1. First note that the system is block separated, i.e. every pik depends only on variables εi and not on all variables ε = (ε1 , . . . , εn ). This shows that Jpε (x, ε) is block diagonal, ε JDp (x, ε)   = J1 0 .. . 0 Jµ   . Now we claim that all these blocks are all equal to the identity matrix. To see this, consider their entry ∂εkj [pik ] for i, j = 1, . . . , µ, which is  d d 1 ,i = j x , Λx pk ] = Λi [bj ] = i−1 [pk ] = Λi [ 0 , otherwise dεkj dεjk since d dεj,k pk = d (bj εj,k ) dεj,k = bj . Lemma 4.2. The µs × n Jacobian matrix JDf (x) of the system Df (x) = (f1 , . . . , fµn ) is of full rank n at x = ζ. Proof of Lemma 4.2. Suppose that the matrix is rankdeficient. Then there is a non-trivial vector in its kernel, JDf (ζ) · v = 0. The entries of v are indexed by ∂i . This implies that a nonzero differential ∆ = v1 ∂1 + · · · + vn ∂n of order one satisfies the following relations: (∆Λi )ζ [fj ] = 0, i = 1, . . . , µ, j = 1, . . . , s. By the standard derivation rules, we have d d (∆Λi ) = vk Λi + ∆ Λi , d∂k d∂k for i = 1, . . . , µ, , k = 1, . . . , n. Since D is stable by derivation, d∂d Λi ∈ D. We deduce that the vector space spanned k by hD, ∆Di is stable by derivation and vanishes on f at ζ. By Proposition 2.2, we deduce that ∆D ⊂ D. This is a contradiction, since ∆ is of degree 1 and the the elements in D are of degree ≤ N . The columns of JDg (x, ε) are indexed by the variables (x, ε), while the rows are indexed by the polynomials gik . We construct the following systems: (a) Let Df I be a subsystem of Df such that the corresponding n rows of JDf (ζ) are linearly independent (Lem. 4.2). We denote by I = {(i1 , k1 ), . . . , (in , kn )} their indices. (b) Let Dg̃(x, ε̃) be the square system formed by removing the variables εi1 ,k1 , . . . , εin ,kn from Dg(x, ε). Therefore the Jacobian JDg̃ (x, ε̃) derives from JDg (x, ε), after purging the columns indexed by εi1 ,k1 , . . . , εin ,kn , T and it’s (ij , kj ) row becomes [∇(Λx ij g̃ij ,kj ) | 0 ]. Theorem 4.3 (Deflation Theorem 1). 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 = ζ. Proof. By construction, ζ is a solution of Df I (x) = 0. Moreover, the indices I are chosen such that det JDf I (ζ) 6= 0. This shows that ζ is a simple (thus isolated) root of the system Df I (x) = 0. Example 4.4. In our running example, we expand the rectangular Jacobian matrix of 6 polynomials in (x1 , x2 ). Choosing the rows corresponding to f1 and (d1 − d22 − d1 d2 − d21 )[f1 ], we find a non-singular minor, hence the resulting system (f1 , 2x1 ) has a regular root at ζ = (0, 0). 2 The deflated system Df I (x) = 0 is a square system in n variables. Contrarily to the deflation approach in [7, 3], we do not need to introduce new variables here and one step of deflation is sufficient. In the following theorem, we do introduce new variables to express the condition that the perturbed system has a given multiplicity structure. Theorem 4.5 (Deflation Theorem 2). Let f (x) be a n−variate polynomial system with an µ−fold isolated zero at x = ζ. The square system Dg̃(x, ε̃) = 0, as defined in (b), has a regular isolated root at (x, ε̃) = (ζ, 0). Proof. By definition of D, we have Dg̃(ζ, 0) = (Λζ1 [f ], . . . , Λζµ [f ]) = 0. Moreover, by construction of Dg̃ we get, up to a row permutation, the determinant: ±det JDg̃ (ζ, 0) = det J1 J2 0 I = det J1 6= 0, where J1 = JDf I (ζ). This shows that (ζ, 0) is regular and thus isolated point of the algebraic variety defined by Dg̃(x, ε̃) = 0. Nevertheless, this deflation does differ from the deflation strategy in [7, 3]. There, new variables are added that correspond to coefficients of differential elements, thus introducing a perturbation in the approximate dual basis, in case of approximate input. Hence the output concerns a deflated root of the given approximate system. In our method, we perturb the equations, keeping an approximate structure of a multiple point. Consequently, the certification of a root concerns a nearby system, within controlled error bounds, as it shall be described in Sect. 5. We mention that it would also be possible to use the equations (9, i-iii) to construct a deflated system on the differentials and to perturb the approximate dual structure. 5. VERIFYING APPROXIMATE SINGULAR POINTS 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, eg. [13, 10]. In the regular case, Smale’s α−theory, extending Newton’s method, can be used to answer this problem. Another option is Rump’s Theorem, also based on Newton theory. In our implementation we choose this latter approach, since it is suitable for inexact data and suits best with the perturbation which is applied. Our perturbation coincides to the numerical scheme of [16] in the univariate case. The certification test is based on the verification method of Rump [16, Th. 2.1], which we rewrite in our setting: Theorem 5.1 ([16] Rump’s Theorem). Let f ∈ Rn be a polynomial system and ζ ∗ ∈ Rn a real point. Given an interval domain Z ∈ IRn containing ζ ∗ ∈ Rn , and an interval matrix M ∈ IRn×n whose i−th column Mi satisfies ∇fi (Z) ⊆ Mi for i = 1 . . . , n, then the following holds: If the inclusion ◦ V (f , Z, ζ ∗ ) = −Jf (ζ ∗ )−1 f (ζ ∗ ) + (I − Jf (ζ ∗ )−1 M )Z ⊆Z (14) is true, then there is a unique ζ ∈ Z with f (ζ) = 0 and the Jacobian matrix Jf (ζ) ∈ M is non-singular. This theorem is applied on the perturbed system. If the test succeeds, we also get a domain for ε−variables that reflects the distance of the approximate system from a precise system with the computed local structure. Example 5.2. We start with an approximate system: f1 = 1.071x1 −1.069x2 +1.018x21 , f2 = 1.024x1 −1.016x2 +1.058x22 and the domain: Z = [−.01, .03] × [−.03, .01]. The Jacobian at x = (0, 0) evaluates to .00652, hence it is close to singular. We set the tolerance equal to .04, i.e. the size of the domain, and we consider the center of the box as our approximate point, ζ ∗ = (.01, −.01). First we compute approximate multiplicity structure at ζ ∗ , D = (1, d2 + 1.00016d22 + .99542d1 d2 + 1.03023d21 , d1 − 1.00492d22 − 1.00016d1 d2 − 1.03514d21 ) as well as (1, x2 , x1 ), a standard basis for the quotient. The latter indicates to perturb up to linear terms. Now apply the second deflation theorem 4.5 to get the 6×6 system g = (1.018 x21 +1.071 x1 +(ε12 − 1.069) x2 , ε12 − .02023, .01723 + 2.036 x1 , 1.058 x22 + (1.024 + ε23 )x1 + (ε22 − 1.016)x2 + ε21 , .04217 + 2.116 x2 + ε22 , ε23 − .03921), which has a regular root for ζ ∈ Z and parameters (ε12 , ε21 , ε22 , ε23 ). Indeed, applying Theorem 5.1 with Z ′ = Z ×[−.04, .04]4 and ◦ (ζ ∗ , 0, .., 0) we get an inclusion V (g, Z ′ , ζ ∗ ) ⊆Z ′ . 2 6. GEOMETRY AROUND A SINGULARITY As a final step in analyzing isolated singularities, we show how the local basis can be used to compute the topological degree around the singular point. If the latter is a selfintersection point of a real algebraic curve, one can deduce the number of curve branches that pass through it. Topological degree computation. Let f (x) be a square n−variate system with an µ−fold isolated zero at x = ζ. To a linear form Λ ∈ R[∂ ζ ], we associate the quadratic form QΛ : R/Q × R/Q → R , (bi , bj ) 7→ Λ(bi bj ) (15) for R/Q = hb1 , . . . , bµ i. The signature of this (symmetric and bi-linear) form is the sum of signs of the diagonal entries of any diagonal matrix representation of it. Proposition 6.1 ([4, Th. 1.2]). If QΦ , Φ ∈ D is any bi-linear symmetric form such that Φζ [det Jf (x)] > 0, then tdegζ (f ) = sgn(QΦ ). (16) This signature is independent of the bi-linear form used. We can use this result to compute the topological degree at x = ζ using the dual structure at ζ. Since a basis D is available we set Φ = ±Λi , for some basis element that is not zero on det Jf (x). Indeed, such an element can be retrieved among the basis elements, since det Jf ∈ / hf i, see [5, Ch. 0]. In practice it suffices to generate a random element of D, compute it’s matrix representation [Φ(bi bj )]ij , and then extract the signature of QΦ . Branches around a singularity. In the context of computing with real algebraic curves, the identification of selfintersection points is only the first step of determining the local topology. As a second step, one needs to calculate the number of branches attached to the singular point ζ, hereafter denoted Br(f , ζ). This information is encoded in the topological degree. An implicit curve in n−space is given by all points satisfying f (x) = 0, f = (f1 , . . . , fn−1 ). Consider p(x) = (x1 − ζ1 )2 + · · · + (xn − ζn )2 , and g(x) = det J(f ,p) (x). Then ([19] and references therein): Br(f , ζ) = 2 tdegζ (f , g). (17) This implies an algorithm for Br(f , ζ). First compute the primal-dual structure of (f , g) at ζ and then use Prop. 6.1 to get tdegζ (f , g). Example 6.2. Consider the implicit curve f (x, y) = 0, in xy−plane, with f (x, y) = x4 + 2x2 y 2 + y 4 + 3x2 y − y 3 , that looks like this . We search for the number of branches touching ζ = (0, 0). We compute g(x, y) = J(f,x2 +y 2 ) = 18xy 2 − 6x3 , and then the multiplicity structure of (f, g) at ζ, and we arrive to the standard basis B = (1, y, x, y 2 , xy, x2 , y 3 , xy 2 , x2 y). Among the 9 elements of the dual basis, we find Φ = d3y + 38 d4y + 1 2 2 d d + 38 d4x , having value Φ0 [det J(f,g) (x)] = 54 > 0 on the 8 x y Jacobian determinant. Using B and (6), we get the 9 × 9 matrix representation of QΦ (15) with ij−th entry Φ[NF(xβi xβ j )], and we compute tdegζ (f, g) = sgn QΦ = 3, thus Br(f, (0, 0)) = 6. 2 7. EXPERIMENTATION Our method is developed in Maple. It uses Mourrain’s Integration technique to compute (approximate) dual basis and derive the augmented system of Th. 4.5. Then Rump’s method is used to verify the root. Macaulay’s method is also implemented for testing purposes. Example 7.1. Consider the system [8] of 3 equations in 2 variables f1 = x31 + x1 x22 , f2 = x1 x22 + x32 , f3 = x21 x2 + x1 x22 , and the singular point (0, 0). Suppose that the point is given. Using 3.2 and 3.4 we derive the primal-dual pair D = (1, d1 , d2 , d21 , d1 d2 , d22 , d32 + d31 + d21 d2 − d1 d22 ), where d32 is underlined to show that it corresponds to x32 in the primal standard basis B. 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 x Df of 21 equations. The 21 × 2 Jacobian matrix JDf is of rank 2 and a full-rank minor consists of the rows 4 and 5. Consequently find the system (d21 [f1 ], d1 d2 [f1 ]) = (3x1 , 2x2 ) which deflates (0, 0). Note that even though both equations of the deflated system derive from f1 , the functionals used on f1 are computed using all initial equations. 2 Example 7.2. Let, as in [6, 8], f1 = 2x1 + 2x21 + 2x2 + + x23 − 1, f2 = (x1 + x2 − x3 − 1)3 − x31 , and f3 = + 2x22 + 10x3 + 5x23 + 5)3 − 1000x51 . Point (0, 0, −1) occurs with multiplicity equal to 18, in depth 7. The final matrix size with our method is 206 × 45, 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 for D1 and stop the process. We get the evaluation 1 and 2 first order functionals, which we apply to f1 . We arrive to (1[f1 ], (d2 − d1 )[f1 ], (d1 + d3 )[f1 ]) = (f1 , −4x1 + 4x2 , 2 + 4x1 + 2x3 ) and we check that the Jacobian determinant is 64, thus we have a deflated system only with a partial local structure. 2 2x22 2x31 Table 7 shows computations on the benchmark set of [3]. Multiplicity, matrix sizes at termination step are reported. Sys. cmbs1 cmbs2 mth191 decker2 Ojika2 Ojika3 KSS Caprasse Cyclic 9 DZ1 DZ2 DZ3 µ 11 8 4 4 2 4 16 4 4 131 7 5 Integration 33 × 23 21 × 17 10 × 9 5×5 6×5 24 × 9 569 × 69 34 × 13 369 × 33 1450 × 524 73 × 33 14 × 8 Macaulay 105 × 56 60 × 35 30 × 20 20 × 15 12 × 10 60 × 35 630 × 252 60 × 35 495 × 33 4004 × 1365 360 × 165 30 × 21 Table 1: Benchmark systems from [3]. Acknowledgments. This research has received funding from the EU’s 7th Framework Programme [FP7/2007-2013], Marie Curie ITN SAGA, grant no [PITN-GA-2008-214584]. 8. REFERENCES [1] L. Alberti, B. Mourrain, & J. Wintz. Topology and arrangement computation of semi-algebraic planar curves. CAGD, 25:631–651, November 2008. [2] M.F. Atiyah & I.G. MacDonald. Introduction to Commutative Algebra. Addison-Wesley, 1969. [3] B. H. Dayton & Z. Zeng. Computing the multiplicity structure in solving polynomial systems. In ISSAC ’05,, pp. 116–123, 2005. ACM. [4] D. Eisenbud & H.I. Levine. An algebraic formula for the degree of a c∞ map germ. The Annals of Mathematics, 106(1):pp. 19–44, 1977. [5] M. Elkadi & B. Mourrain. Introduction à la résolution des systèmes d’équations algébriques, vol. 59 of Mathématiques et Applications. Springer, 2007. [6] G. Lecerf. Quadratic newton iteration for systems with multiplicity. Fo. Comp. Math., 2:247–293, 2002. [7] A. Leykin, J. Verschelde, & Zhao A. Newton’s method with deflation for isolated singularities of polynomial systems. TCS, 359(1-3):111 – 122, 2006. [8] A. Leykin, J. Verschelde, & A. Zhao. Higher-order deflation for polynomial systems with isolated singular solutions. vol. 146 of The IMA Volumes in Math. and its Appl., pp. 79–97, 2008. [9] F.S. Macaulay. The algebraic theory of modular systems. Cambridge Univ. Press, 1916. [10] A. Mantzaflaris, B. Mourrain, & E. Tsigaridas. Continued fraction expansion of real roots of polynomial systems. In Proc. of SNC ’09, pp. 85–94, 2009. [11] M. G. Marinari, T. Mora, & H.M. Möller. Gröbner duality and multiplicities in polynomial system solving. In Proc. of ISSAC ’95, pp. 167–179, 1995. [12] B. Mourrain. Isolated points, duality and residues. J. of Pure & App. Alg., 117-118:469 – 493, 1997. [13] B. Mourrain & J. P. Pavone. Subdivision methods for solving polynomial equations. J. Symb. Comp., 44:292–306, March 2009. [14] T. Ojika, S. Watanabe, & T. Mitsui. Deflation algorithm for multiple roots of a system of nonlinear equations. J. of Math. An. & Appls., 96(2):463–479, 1983. [15] S.R. Pope & A. Szanto. Nearest multivariate system with given root multiplicities. JSC,44(6):606-625,2009. [16] S. Rump & S. Graillat. Verified error bounds for multiple roots of systems of nonlinear equations. Num. Algs., 54:359–377, 2010. [17] Michael Shub & Steve Smale. Complexity of Bezout’s theorem I: Geometric aspects. J. of the AMS, 6(2):459–501, 1993. [18] H. J. Stetter. Analysis of zero clusters in multivariate polynomial systems. In Proc. of ISSAC ’96, pp. 127–136, New York, NY, USA, 1996. ACM. [19] Zbigniew Szafraniec. Topological degree and quadratic forms. J. of Pure & App. Alg., 141(3):299 – 314, 1999. [20] X. Wu & L. Zhi. Computing the multiplicity structure from geometric involutive form. In Proc. of ISSAC ’08, pp. 325–332, 2008. ACM. [21] Z. Zeng. The closedness subspace method for computing the multiplicity structure of a polynomial system. vol. 496 of Contemporaty Mathematics, pp. 347–362. AMS. Providence, RI, 2009. APPENDIX We attach additional examples that did not fit page limits. Verification Example. Let f1 = (x21 x2 − x1 x22 , f2 = x1 − x22 ). The verification method of [16] applies a linear perturbation on this system, but fails to certify the root x = (0, 0). We consider an approximate point ζ ∗ = (.01, .002) and we compute the approximate multiplicity structure: D = (Λ1 , . . . , Λ4 ) = (1.0, 1.0d2 , 1.0d1 +1.0d22 , 1.0d1 d2 +1.0d32 ) The augmented system g(x) = (Λj (fi )) = (f1 , 2.0x1 x2 − 1.0x22 −1.0x1 , 2.0x1 −2.0x2 , 1.0x1 −1.0x22 , f2 , −2.0x2 , 0., 0.) has a Jacobian matrix:   .00 .016 −.99 2.0 1.0 0 0 0 Jg (ζ ∗ )T = .00 −.02 .016 −2.0 −.004 −2.0 0 0 Example 6.2 (Cont’d). Consider the implicit curve f (x, y) = 0, in xy−plane, with f (x, y) = x4 + 2x2 y 2 + y 4 + 3x2 y − y 3 , . We search for the number of branches that looks like this touching ζ = (0, 0). This point is of multiplicity 4, as the dual basis we get for f (x, y) = d d f (x, y) = f (x, y) = 0 dx dy is (1, dx , dy , d2x + d2y ), which provides no information for the number of branches. We compute g(x, y) = J(f,x2 +y 2 ) = 18xy 2 − 6x3 , f5 = x1 − x22 + ε21 + ε22 x2 + ε23 x1 + ε24 x1 x2 and then the multiplicity structure of (f, g) at ζ, and we arrive to 3 1 3 D =(1, dy , dx , d2y , dx dy , d2x , d3y + d4y + d2x d2y + d4x , 8 8 8 9 4 3 2 2 9 4 2 3 2 dx dy + 3 dx , dx dy − dy − dx dy − dx ), 8 8 8 and the standard basis Thus g(x1 , x2 , ε11 , ε12 , ε21 , ε22 , ε23 , ε24 ), computed as before, is a square system with additional equations B = (1, y, x, y 2 , xy, x2 , y 3 , xy 2 , x2 y). with a non-zero minor at the third and forth row. Using this information, we apply the following perturbation to the original system: f1 = x21 x2 − x1 x22 + ε11 + ε12 x2 Among the 9 elements of the dual basis, we find f2 = 1.0x21 − 2.0x1 x2 + 1.0ε12 f3 f4 f6 f7 f8 Φ = d3y + = 2.0x1 x2 − 1.0x22 − 1.0x1 = 2.0x1 − 2.0x2 = −2.0x2 + 1.0ε22 + 1.0x1 ε24 = 1.0ε23 + 1.0x2 ε24 = 1.0ε24 Now take the box Z1 = [−.03, .05]×[−.04, .04]×[−.01, .01]6 . We apply Th. 5.1 on g, ie. we compute V (g, Z1 , ζ ∗ ). For the variable ε21 the interval is [−.015, .15] 6⊆ (−.01, .01), therefore we don’t get an answer. We shrink a little Z1 down to Z2 = [−.03, .05]×[−.02, .02]× [−.01, .01]6 and we apply again Th. 5.1, which results in   [−.004, .004]   [−.004, .004]   [−.001, .001]    ◦  [−.007, .007]   ∗ V (g, Z2 , ζ ) =   ⊆Z 2 , [−.006, .006]     [−.009, .009]    [−.00045, .00035]  [.0, .0] thus we certify the multiple root of the original system inside Z2 .             1 y x y2 xy x2 y3 xy 2 x2 y y y2 xy y3 xy 2 x2 y 3/8y 3 − 89 x2 y 0 1/8y 3 − 38 x2 y x xy x2 xy 2 x2 y 3xy 2 0 1 3 3 2 y − x y 8 8 0 y2 y3 xy 2 3 3 y − 89 x2 y 8 0 1 3 y − 3/8x2 y 8 0 0 0 xy xy 2 x2 y 0 1 3 y − 83 x2 y 8 0 0 0 0 Multiplication table for R/hf, gi (Example 6.2). 3 4 1 2 2 3 4 d + d d + d , 8 y 8 x y 8 x having value Φ0 [det J(f,g) (x)] = 54 > 0 on the Jacobian determinant. Using B and (6), we compute the matrix NF(xβi xβj ) of multiplication in R/hf, gi, given at the end of the page. Now a representation of QΦ (15) can be computed, by applying Φ0 on the multiplication table to get:        QΦ =      0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 3/8 0 1/8 0 0 0 0 0 0 0 1/8 0 0 1 0 3/8 0 1/8 0 0 0 0 0 0 0 1/8 0 0 0 0 0 0 0 1/8 0 3/8 0 0 0 1 3/8 0 0 0 0 0 0 0 0 0 1/8 0 0 0 0 0 0 0 1/8 0 0 0 0 0 0 0      .     With a QR iteration, we find 6 positive and 3 negative eigenvalues of this representation, hence we compute tdegζ (f, g) = sgn QΦ = 6 − 3 = 3, i.e. there are 6 branches of the curve around (0, 0). x2 x2 y 3xy 2 1 3 y − 83 x2 y 8 0 3/8y 3 − 89 x2 y 0 0 0 3 3 y 8 y3 − 89 x2 y 0 0 0 0 0 0 0 xy 2 0 1/8y 3 − 3/8x2 y 0 0 0 0 0 0 1 3 y 8 x2 y − 83 x2 y 0 0 0 0 0 0 0            