Deflation and Certified Isolation of Singular Zeros of Polynomial Systems

0.23 MB·English

by  Angelos Mantzaflaris

Deflation and Certified Isolation of Singular Zeros of
Polynomial Systems
Angelos Mantzaflaris and Bernard Mourrain
GALAAD,INRIAMéditerranée
BP93,06902SophiaAntipolis,France
[FirstName.LastName]@inria.fr
1
ABSTRACT surfacesingeometricmodeling. Whenalgebraicrepresenta-
1
tions are used, this reduces to solving polynomial systems.
0 Wedevelopanewsymbolic-numericalgorithmforthecerti-
Several approaches are available: algebraic techniques such
2 ficationofsingularisolatedpoints,usingtheirassociatedlo-
calringstructureandcertifiednumericalcomputations. An asGr¨obnerbasesorborderbases,resultants,subdivisional-
n gorithms[10], [13], homotopies, and so on. At the end of
improvementofanexistingmethodtocomputeinversesys-
a the day, a numerical approximation or a box of isolation is
J temsispresented,whichavoidsredundantcomputationand
usuallycomputedtocharacterizethe(real)rootsofthepoly-
reduces thesize of theintermediate linear systems to solve.
7 nomialsystem. Butweoftenneedtoimprovethenumerical
We derive a one-step deflation technique, from the descrip-
1 approximation of the roots. Numerical methods such that
tion of the multiplicity structure in terms of differentials.

] The deflated system can be used in Newton-based iterative Newton’s iteration can be used to improve the quality of
theapproximation,providedthatwehaveasimpleroot. In
C schemes with quadratic convergence. Starting from a poly-
the presence of a multiple root, the difficulties are signifi-
S nomialsystemandasmall-enoughneighborhood,weobtain
cantly increasing. The numerical approximation can be of
. acriterionfortheexistenceanduniquenessofasingularroot
s ofagivenmultiplicitystructure,applyingawell-chosensym- verybadquality,andthemethodsusedtocomputethisap-
c proximationareconvergingslowly (ornotconverging). The
bolic perturbation. Standard verification methods, based
[
situation in practical problems, as encountered in CAGD
eg. on interval arithmetic and a fixed point theorem, are
for instance, is even worse, since the coefficients of the in-
1 employed to certify that there exists a unique perturbed
v system with a singular root in the domain. Applications to putequationsareknown,withsomeincertitude. Computing
0 topological degree computation and to the analysis of real multipleroots of approximatepolynomial systemsseems to
4 branchesof an implicit curveillustrate themethod. be an ill-posed problem, since changing slightly the coeffi-
1 cientsmaytransformamultiplerootintoaclusterofsimple
3 Categories andSubject Descriptors roots (oreven makeit disappear).
. To tackle this difficult problem, we adopt the following
1 G.1.5 [Mathematics of Computing]: RootsofNonlinear
strategy. Wetrytofinda(small) perturbationof theinput
0 Equations; I.1.2 [Computing Methodologies]: Symbolic
system such that the root we compute is an exact multi-
1 and Algebraic Manipulation—Algebraic algorithms
ple root of this perturbed system. In this way, we identify
1
General Terms themultiplicitystructureandweareabletosetupdeflation
:
v techniques which restore the quadratic convergence of the
Algorithms, Theory
i Newtonsystem. Thecertificationofthemultiplerootisalso
X
Keywords possible on the symbolically perturbed system by applying
r a fixed point theorem, based eg. on interval arithmetic [16]
a rootdeflation,multiplicitystructure,dualspace,inversesys-
or α-theorems [17].
tem, isolated point
Related work. InordertobedevelopNewton-typemeth-
1. INTRODUCTION ods that converge to multiple roots, deflation techniques
A main challenge in algebraic and geometric computing which consist in adding new equations in order to reduce
is singular point identification and treatment. Such prob- the multiplicity have already been considered. In [14], by
lems naturally occur when computing the topology of im- applyingatriangulationpreprocessingstepontheJacobian
plicit curves or surfaces [1], the intersection of parametric matrixattheapproximateroot,minorsoftheJacobianma-
trix are added to thesystem to reducethemultiplicity.
In [6], a presentation of the ideal in a triangular form in
a good position and derivations with respect to the leading
variablesareusedtoiterativelyreducethemultiplicity. This
Permission tomake digital orhardcopies ofall orpartofthis workfor processisappliedforp-adicliftingwith exactcomputation.
personalorclassroomuseisgrantedwithoutfeeprovidedthatcopiesare In [7, 8], instead of triangulating the Jacobian matrix,
notmadeordistributedforprofitorcommercialadvantageandthatcopies the number of variables is doubled and new equations are
bearthisnoticeandthefullcitationonthefirstpage.Tocopyotherwise,to
introduced,whicharelinearinthenewvariables. Theydes-
republish,topostonserversortoredistributetolists,requirespriorspecific
cribethekerneloftheJacobianmatrixatthemultipleroot.
permissionand/orafee.
Copyright200XACMX-XXXXX-XX-X/XX/XX...$5.00. In [3], this construction in related to the construction of
the inverse system. The dialytic method of F.S. Macaulay (b) Deflate the system by augmenting it with new equa-
[9]isrevisitedforthispurpose. Thesedeflationmethodsare tionsderivedfromthedualbasis,introducingadequateper-
applied iteratively until the root becomes regular, doubling turbation terms.
each time thenumberof variables. (c)Certifythesingularpointanditsmultiplicitystructure
Morerecentalgorithmsfortheconstructionofinversesys- for theperturbedsystem checkingthecontraction property
tems are described eg. in [11], reducing the size of the in- of Newton iteration (eg. viainterval arithmetic).
termediate linear systems (and exploited in [18]), or in [12]
Instep(a),adualbasisatthesingularpointiscomputed
using an integration method.
by means of linear algebra, based on the integration ap-
In [15], a minimization approach is used to reduce the
proachof[12]. Wedescribeanimprovementofthismethod,
value of the equations and their derivatives at the approxi-
which yields directly a triangular dualbasis with noredun-
materoot,assumingabasisoftheinversesystem isknown.
dant computation. This method has the advantage to re-
In [20], the inverse system is constructed via Macaulay’s
duce significantly the size of the linear systems to solve at
method;tablesofmultiplicationsarededucedandtheireigen-
each step, compared to Macaulay’s type methods [9, 7, 8,
values are used to improve the approximated root. It is
3]. In the case of an approximate singular point, errors are
provedthattheconvergenceisquadraticwhentheJacobian
introduced in the coefficients of the basis elements. Yet a
has corank 1 at themultiple root.
successful computation is feasible. In particular, the sup-
Verificationofmultiplerootsof(approximate)polynomial
port of the basis elements is revealed by this approximate
equations is a difficult task. The approach proposed in [16]
process.
consists in introducingperturbation parameters and tocer-
Inthedeflationstep(b),newequationsandnewvariables
tifying the multiple root of nearby system by using a fixed
areintroducedinordertoarrivetoanewpolynomialsystem
pointtheorem, basedonintervalarithmetic. Itappliesonly
where the singularity is obviated. The new variables corre-
to cases where theJacobian has corank equal to1.
spondtoperturbationsoftheinitialequationsalongspecific
The univariate case. In preparation for the multivariate polynomials, which form a dual counterpart to the basis of
case,wereviewsometechniquesusedtotreatsingularzeros thedualspace. Oneofthedeflatedsystemsthatwecompute
of univariate polynomials, and we present our method on a from the dual system is a square n n system with a sim-
univariate instance. ple root. This improves the deflatio×n techniques described
Let g(x) K[x] be a polynomial which attains at x = 0 in [7, 8, 3], which require additional variables and possi-
∈
a root of multiplicity µ > 1. The latter is defined as the bly several deflation steps. New variables are introduced
smallest positive integer µ such that dµg(0) = 0 whereas only in the case where we want to certify the multiplicity
g(0) = dg(0) = = dµ−1g(0) = 0. Her6e we denote structure. The perturbation techniques that we use extend
···
dkg(x)= dk g(x)/k!thenormalized k th order derivative. the approach of [16] to general cases where the co-rank of
dxk −
We see that D = 1,d1,...,dµ−1 is the maximal space the Jacobian matrix could be bigger than one. The verifi-
0
of differentials whichhisstable underiderivation,that vanish cation step (c) is mostly a contraction condition, using eg.
when applied to members of , the x primary compo- techniques as in [16]. This step acts on the (approximate)
0
nent of g at x = 0. ConsidQer now thhei−symbolically per- deflated system, since verifying a simple solution of the de-
turbed ehquiation flatedsysteminducesacertificateofanexactsingularpoint
of (a nearby to) theinitial system.
f1(x,ε)=g(x)+ε1+ε2x+ +εµ−2xµ−2 (1) Weare going to detail thedifferentsteps in thefollowing
···
andapplyeverybasiselementof toarrivetothenewsys- sections, starting with notations in Sect. 2, dual basis in
0
D Sect. 3, deflation in Sect. 4, and certification in Sect. 5. In
temf(x,ε)= f ,d f ,...,dµ−1f inµ 1variables. The
(cid:16) 1 1 1 1(cid:17) − the last section, we will show examples and applications to
ii.−etlhineeqauraitnioεn,sthhaelllbasetfoin=edbie−i1nfg1f=µd=i−d1µg−+1Pg(xµk=−).i2Txkh−isi+s1yεsk-, thetopology analysis of curves.
tem deflatestheroot, aswe seethatthedeterminantofit’s
2. PRELIMINARIESANDMAINRESULTS
Jacobian matrix at (0,0) is
We denote by R = K[x] a polynomial ring over the field
ddxf1 1 0 Kofcharacteristic zero. Also, thedual ring R∗ isthespace
detJf(0,0)=(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ddxddxf...µf−µ1 0 .0.. 1 (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ==−−µµddfµµg((00,)06=) 0. asoaerfsrelitiefhnsoeeriamnsrpadafluicffvneeacrorteifianofbtnoilaarelmlss.oaΛpTlehs:reuaRrstioew→rsseKaKv[ti[.e∂aw]I]ptdowiuisnhateclroζeemle∂mmKo=ennnlt.(ys∂T1aiod,s.esf.npo.tre,imcfi∂ienfayd)l
(cid:12) (cid:12) ∈
Nowsupposet(cid:12)hatζ∗ isanapproximatez(cid:12)ero,closetox=ζ. thatweusethepointζ,wealsodenotethesedifferentials∂ζ.
tWiveescaunpsttiollacotmhpreustheoDldζrbeylaetvivaelutaotinthgeg(exrr)oarnidntζh∗e.deTrihvean- wWehwenilladpepnlyoitnegbΛy(Λ∂ζζ[)g]∈=KΛ[[ζ∂gζ=]] tΛo(∂aζp)o[gly(xn)o]mtihaelogp(xer)a∈tioRn
we can form (1) and use verification techniques to certify
the root. Checking that the Newton operator is contract- Λζ[g]= λα d|α|g (ζ), (2)
ing shows the existence and unicity of a multiple root in a αX∈Nn α1!···αn! · dα11···dαnn
neighborhood of the input data. We are going to extend
tlahtiesdamppurlotiapclhe,rdooestcs.ribed in [16], to multi-dimensional iso- tfooraΛn(∂orζd)e=rePd sλetαα1!∂=αζ(Λ∈1K,.[[.∂.,ζΛ]].µ)ExteKn[d[∂in]]gµt,hwisedsehfianliltdioen-
D ∈
Our approach. It consists of thefollowing steps: note ζ[g]=(Λζg,...,Λζg). Insomecases,itisconvenient
(a) Compute a basis for the dual space and of the local touseDnormalize1ddifferenµtialsinsteadof∂: foranyα Nn,
quotient ring at a given (approximate) singular point. we denote dαζ = α1!∂αζ. When ζ = 0, we have dα0xβ∈= 1
if α=β and 0 otherwise. More generally, (dαζ)α∈Nn is the LetsuppD0 bethesetofexponentsofmonomialsappear-
dual basis of ((x ζ)α)α∈Nn. ing in D0, with a non-zero coefficient. These are of degree
For Λ R∗ an−d p R, let p Λ : q Λ(pq). We check atmost N,thenilindexof 0. Since Λ D0,Λ0[p]=0iff
that ∈ ∈ · 7→ p 0, we derive that supQpD0 = α∀: x∈α / 0 . In par-
∈Q { ∈Q }
d ticular,wecanfindabasisofR/ 0 betweenthemonomials
(xi−ζi)·∂αζ = d∂ (∂αζ). (3) xα : α suppD . Thisis afinQiteset of monomials, since
i,ζ { ∈ }
their degreeis boundedbythenilindex of 0.
This property shall be useful in thesequel. Givenastandardbasis =(xβi)i=1,...,µQofR/ 0and,for
2.1 Isolatedpoints and differentials all monomials xγj / 0,Bj = 1,...,s µ, s =Q#suppD0,
with xγj / , thee∈xpQression (normal fo−rm) of
Let f = f1,...,fs be an ideal of R, ζ K a root of ∈B
h i h i ∈
f and m = x ζ ,...,x ζ the maximal ideal at ζ. µ
Supposeζthathζ1is−an1isolatedn−roontioff,thenaminimalpri- xγj = λijxβi mod 0 (4)
Q
mary decomposition of = contains a primary Xi=1
I Q
Qpr\im.⊃I in thebasis then thedualelements [12, Prop. 13]
B
component suchthat =m and√ ′ m forthe
ζ ζ ζ ζ
othAesr prima=Qrymcom,Rp/onenitpssaQQfi′naistseodciimateendsitoonQIal[6⊂v2e].ctorspace. Λi =dβi +s−µλijdγj, (5)
ThempulQtiζplicityζµζ ofQζζisdefinedasthedimensionofR/ ζ. Xj=1
Apointofmultiplicityoneiscalledregularpoint,orsimQple for i=1,...,µ form a basis of D0. We give a proof of this
root,otherwisewesaythatζ isasingularisolated point,or fact in thefollowing lemma.
multiple root of f. In thelatter case we haveJ (ζ)=0.
f
Wecan now definethedual space of an isolated point. Lemma 2.4. Thesetofelements =(Λi)i=1,...,µ isaba-
sis of D0 and the normal form of anDy g(x) R with respect
eleDmeefnitnsiotfioKn[[∂2.1].] tThhaet vdaunailshspoancealolfthIeiesletmheensutsbsopface. oItf to the standard basis B=(xβi)i=1,...,µ is ∈
ζ
is also called the orthogonal of I and denoted by I⊥. I NF(g)= µ Λ0[g]xβi. (6)
i
Consider now the orthogonal of ζ, i.e. the subspace Dζ Xi=1
of elementsof R∗ that vanish on mQembers of , namely
Qζ Proof. Firstnotethattheelementsof arelinearlyin-
Q⊥ζ =Dζ ={Λ∈R∗ : Λζ[p]=0, ∀p∈Qζ}. dfoerpeanlldxenαt./Now0,,ebgy. cNonFs(txruβcit)io=n,xPβiµi.=1AΛlsD0io[,xfαo]r=xαNF(xα0),
Thefollowinglemmaisanessentialpropertythatallowsex- i, Λ0(xα)∈=Q0, since α / supp . Thus the elements∈oQf
traction of thelocal structureDζ directly from the“global” ∀compiuteNF(·)onallmo∈nomialsDofR,and(6)followsbylDi-
idealI =hfi,notablybymatrixmethodsoutlinedinSect.3. nearity. Wededucethat isabasisofD0,asinDef.2.1.
D
Proposition 2.2 ([12, Th. 8]). Foranyisolatedpoint Computing the normal form of the border monomials of
B
ζ K of f, we have f ⊥ K[∂ ]=D . via (6) also yields the border basis relations and the opera-
ζ ζ
∈ h i ∩ tors of multiplication in the quotient R/ 0 (see eg. [5] for
In other words, we can identify Dζ =Q⊥ζ with thespace of more properties). Q
polynomial differential operators that vanish at ζ on every Ifagradedmonomialorderingisfixedand =(xβi)i=1,..,µ
eleTmheentspoafcIe.Dζ is a vector space of polynomials in ∂ζ di- ilseatdhiencgotrerremspoonfd(i5n)gwstratntdhaerodpbpaossiisteofoRrd/eQri0n,Bgth[8e,nTdhβ.i3is.1t].he
amcetnosnionR∗µζa,stdheerimvautlitoipnlsic(isteyeo(f3ζ)).,ADsζtihseavsapraiacbeleosf(dxiffie−reζni)- effiCcoiennvtermsealtyr,ixifinwethaeredugiavlenmoanobmasiiaslsDbaosfisD(0dαw)hose cois-
tial polynomials in ∂ζ, which is stable by derivation. This D Kµ×s, we can compute a basis of R/ 0 by choαo∈/sQin0g µ
propertywillbeusedexplicitlyinconstructingDζ (Sec. 3). ind∈ependent columns of D, say those indQexed by dβi, i =
1,...,µ . If G Kµ×µ is the (invertible) matrix formed by
tegDereNfinitNiosn.t.2.m3.NThe ni.lindex of Qζ is the maximal in- these columns,∈thenD′ :=G−1D, is
∈ ζ 6⊂Qζ
β β γ γ
It is directly seen that themaximal degree of an element of 1 ··· µ 1 ··· s−µ
Dζ has tobe equalto N, also known as thedepth of Dζ. Λ′1  1 0 λ1,1 ··· λ1,s−µ
2.2 Quotient ring anddual structure D′ = ...  ... ... ... , (7)
rinIgnatnhdistsheectqiounotwieentexRp/lorewthheerreelatioins tbheetwpereimnatrhyecdouma-l Λ′µ  0 1 λµ,1 ··· λµ,s−µ
ζ ζ
ponent of the isolated poQint ζ. WeQshow how to extract a i.e. abasisoftheform(5). NotethatanarbitrarybasisofD
basis of this quotient ring from the support of the elements does not have theabove diagonal form, nor does it directly
of Dζ and how Dζ can be used to reduce any polynomial providea basis for R/ 0.
modulo . For t N, D denotQes the vector space of polynomials of
ζ t
ItiscoQnvenientintermsofnotationtomaketheassump- D of deg∈ree t. The Hilbert function h:N N is defined
tion ζ =0. This saves some indices, while it poses no con- by h(t)=dim≤(D ),t 0, hence h(0)=1 and→h(t)=dimD
t
≥
straint (sinceit impliesalinearchangeof coordinates),and fort N. Theintegerh(1) 1=corankJ isknownasthe
f
shall beadopted hereafter and in thenextsection. bread≥th of D. −
3.COMPUTINGLOCALRINGSTRUCTURE The kernel of the left matrix gives = (1,d +d ). Ex-
1 1 2
D
The computation of a local basis, given a system and a panding up to order two, we get the matrix on the right,
panodintc,onissedqouneentelysseitntciaalnlybbeycamrraiterdixo-kuetrnneulmceormicpaulltya,tieovnesn, aunpdtoDd2e=pth(13,dw1e+gde2t,t−hde1sa+mde21n+uldl-1sdp2a+ce,dt22h).usIf w=e ex2p.an2d
D D
when the point or even the system is inexact. Throughout
3.2 Integration method
thesection we supposef Rm and ζ Kn with f(ζ)=0.
∈ ∈
Severalmatrixconstructionshavebeenproposed,thatuse This method is presented in [12]. It is an evolution of
differentconditionstoidentifythedualspaceasanull-space. Macaulay’s method, in the sense that the matrices are not
They are based on thestability property of thedual basis: indexed by all differentials, but just by elements based on
knowledge of the previous step. This performs a computa-
d
∀Λ∈Dt, d∂iΛ∈Dt−1 , i=1,...,n. (8) ttiroicnesa.dapted to the given input and results in smaller ma-
We list existing algorithms that compute dual-spacebases: For Λ K[∂], we denote by Λ the element Φ K[∂]
Aspointedoutin(3),anequivalentformof (8)is: Λ with the∈property d Φ(∂) = ΛRk(∂) and with no co∈nstant
D•,Λ[x f ]=0, i,j =1,...,n. Macaulay’smethod[9]∀use∈s d∂k
ittto deirivjethea∀lgorithm that is outlined in Sect. 3.1. term wrt ∂k.
In[11]theyexploit (8)byformingthematrixD ofthe Theorem 3.2 ( [12, Th. 15] ). Let Λ ,Λ ,...,Λ be
i 1 2 s
ma•p d : K[∂] K[∂] for all i = 1,...,n and some abasisofDt−1, thatis, thesubspace ofDhofelements ofior-
triangdu∂lair decomtp→osition to−f1the differential polynomials in tdeerrmatliemsoisntDt−iff1.itAins oefletmheenftorΛm:∈ K[∂] with no constant
t
terms the differential variables. This approach was used in
s n
[18] toreducetherow dimension of Macaulay’s matrix, but
Λ(∂)= λ Λ (∂ ,...,∂ ,0,...,0), (9)
not thecolumn dimension. Theclosedness condition is also ik k i 1 k
used in [21] toidentify asuperset of suppD . Xi=1Xk=1 R
The integration method in [12]“integratt+es1”elements of for λik K, and the following two conditions hold:
ab•asisofDt,andobtainsa priori knowledgeoftheform of s∈ d s d
elements in degree t+1 (Sect.3.2). (i) λ Λ (∂) λ Λ (∂)=0,
All methods are incremental, in thesense that they start Xi=1 ikd∂l i −Xi=1 ild∂k i
by setting = (1) and continue by computing , i = for all 1 k l n .
D0 Di ≤ ≤ ≤
o1f,.D..,,aNn,dNN+is1.thWehneilnin#deDxNof=#. DN+1 then DN is a basis (ii) Λζ[fk]=0, for k=1,...,n .
baWsisefsohralDl,reavniedwthtwenodoefscthriebsQee aanppirmoapcrhovesemtoencto,mthpautteala- CTohnudsittihoent(wi)oiscoenqduiitviaolnesntetxoprde∂dsksΛex∈acDtlty−1t,hfeorfa1ct≤tkha≤t nD.
lowssimultaneouscomputationofaquotientringbasiswhile must be stable by derivation and his members must vanish
avoiding redundantcomputations. on f .
h i
This gives the following algorithm to compute the dual
3.1 Macaulay’s dialyticmatrices
basis: Startwith = 1 . GivenabasisofD wegener-
0 t−1
This matrixconstruction is presentedin [9, Ch.4], amo- ate the ns candidDate elhemients Λ (∂ ,...,∂ ,0,...,0).
k i−1 1 k
dern introduction is contained in [3], together with an im- Conditions (i) and (ii) give a linRear system with unknowns
plementation of themethod in ApaTools1. λ . The columns of the corresponding matrix are indexed
ik
The idea behind the algorithm is the following: An ele- by the candidate elements. Then the kernel of the matrix
ment of D is of the form Λ(∂) = λαdα under the givesthepossiblesolutions,whichweaddto t−1 toobtain
D
|αX|≤N t. If for some t there are no further new elements, then
condition: Λ0 evaluates to 0 at any g f , i.e. Λ0(g) = D= t is a basis of D.
Λ0( g f ) = 0 Λ0(xβf ) = 0 fo∈r ahllimonomials xβ. D D
i i ⇐⇒ i Example 3.3. ConsidertheinstanceofEx.3.1. Wehave
If wPe apply this condition recursively for α N we get f (ζ)=f (ζ)=0,thusweset = 1 . Equation(9)gives
| | ≤ 1 2 0
a vector of coefficients (λα)|α|≤N in the (right) kernel of Λ=λ1d1+λ2d2. Condition (iD) indu{ce}s no constraints and
the matrix with rows indexed by constraints Λ0[xβf ]=0, (ii) yields thesystem
i
β N 1.
| N|≤ote t−hat the only requirement is to be able to perform (cid:20) 11 −−11 (cid:21)(cid:20) λλ12 (cid:21)=0 (10)
derivation of theinputequations and evaluation at ζ =0.
wherethecolumnsareindexedbyd ,d . Wegetλ =λ =
Example 3.1. Letf =x x +x2,f =x x +x2. We 1 2 1 2
alsoreferthereaderto[13,Ex.12−]fo2rad1eta2iledd1e−mo2nstr2ation 1 fFroomrththeesekceornndelsotfepth,iwsemcaotmripxu,ttehuthseDe1le=me{n1t,sd1of+d2,}t.hat
of Macaulay’s method on the same instance. The matrices mustbeoftheformΛ=λ d +λ d +λ d2+λ (d Dd2+d2).
in order1 and 2 are: 1 1 2 2 3 1 4 1 2 2
Condition (i) yields λ λ = 0, and together with (ii) we
3 4
1 d1 d2 d21 d1d2 d22 form the system −
ff12 (cid:20)100 d111 −−d211(cid:21) , xxxff11221fff121 00000 01100 −−00011 11001 −−10011 −10001.  011 −−011 110 −011  λλ...41 =0, (11)
x2f2 0 0 0 0 1 −1 with columns indexed by d1,d21,d2,d1d2+d22. We get two
1http://www.neiu.edu/ zzeng/apatools.htm vectors in the kernel, the first yielding again d +d and a
1 2
second one for λ = 1,λ =0,λ =λ =1, so we deduce corresponding to λ instead of adding a row. Either way,
1 2 3 4 ik
that d +d2+d d−+d2 is a new element of . this lemma allows to shrink the kernel of the matrix and
− 1 1 1 2 2 D2
In thethird step we have computeonly new dual elements.
Let us explore our running example, to demonstrate the
Λ=λ d +λ d +λ d2+λ (d d +d2)+
1 1 2 2 3 1 4 1 2 2 essence of thisimprovement.
λ (d3 d2)+λ (d3+d d2+d2d d d ),
5 1− 1 6 2 1 2 1 2− 1 2 Example 3.5. Were-runEx. 3.3 using Lem. 3.4.
aconnddtiotigoenth(ie)rlweaidthstcoonλd3i−tioλn4+(ii()λw5−e λar6r)idv1e+to(λ5−λ6)d2=0, forImntwhiethinrietsiapleicztattioon st=epD10.=Fo(r1t)hiesfiarlrsetasdteypi,nwteridaenmgualnadr
0
B { }
Λ[1] = 0, thus the matrix is the same as (10), yielding
 100 −001 101 −001 −101 −001  λ...1 =0. (12) Dtr1ian=gu(l1a,rd1w+ithd2re).spWecet teoxte1n.d B1 = {1,x2}, so that D1 is
 1 −1 0 1 0 0  λ6  mnIn,htehnecesecond step we reBmove from (11) the second colu-
Since the kernel of this matrix gives elements that are al-
readyin ,wederivethat = = andthealgorithm 0 1 −1 λ1
D2 D D2 D3  1 1 0  λ3 =0,
terminates.
1 0 1 λ4
Notethat for thisexampleMacaulay’s method endswith   
amatrixofsize12 10, insteadof4 6inthisapproach. 2 yieldingasingle solution d +d2+d d +d2. Weextend
× × − 1 1 1 2 2
byaddingmonomial x : = 1,x ,x .
1 1 1 2 1
3.3 Computing a primal-dual pair BFor the final step, we seaBrch a{n elemen}t with Λ[x ] =
1
Λ[x ]=0 thus(12) loses two columns:
2
In this section we provide a process that allows simulta-
neous computation of a basis pair( , ) of D and R/ . 0 0 1 −1 λ3
Computing a basis of D degreeDbyBdegree involvesQdu-  11 −01 −01 00  ... =0.
plicatedcomputations. Thesuccessivespacescomputedare
D D =D . Itismoreefficienttocomputeonly  0 1 0 0  λ6 
1 N N+1
new⊂el·e·m·⊂ents Λ Dt which are independent in Dt/Dt−1 at We find an empty kernel, thus we recover the triangular
∈
step t. basis = 2,whichisthendiagonalized toreachtheform:
D D
Also,oncedualbasisiscomputed,onehastotransformit
intotheform(5),inordertoidentifyabasisofR/ aswell. 1 d2 d1 d21 d1d2 d22
This transformation can be done a posteriori, byQfinding a Λ1 1 0 0 0 0 0
sub-matrix of full rank and then performing Gauss-Jordan ΛΛ32 00 01 10 −11 −11 −11.
elimination overthis sub-matrix,to reach matrix form (7).  
Weintroduceacondition(iii)extendingTh.3.2,thatad- This diagonal basis is dual to the basis = (1,x2,x1) of
B
dressesthesetwoissues: Itallowsthecomputationofatotal thequotientringandalsoprovidesanormalformalgorithm
of µ independent elements throughout execution, and re- (Lem. 2.4) wrt . In the final step we generated a 4 4
turnsa“triangular basis”, e.g. a basis of R/ is identified. matrix, size smaBller compared toall previousmethods. ×2
Q
Lemma 3.4. Let t−1 =(Λ1,...,Λk) be a basis of Dt−1, This techniquefor computing can be applied similarly to
whose coefficient maDtrix is otherthematrixmethods,e.g.BMacaulay’sdialyticmethod.
β1 ··· βk γ1 ··· γs−k meInfths(itn)s−tehp(tt,−th1e)n>the1,chieo.ictehoefremaorneommioarlesttohaanddontoe eleis-
Λ1 1 ∗ ∗ ∗ ··· ∗ obtained by extracting a non-zero maximal minor fromBthe
...  0 ... ∗ ... ... , (13) cthoeeffimcoiennotmmiaalstroixf simna(dlleαs)t.dIengrperea.ctice, we will look first at
Λk  0 0 1 ∗ ··· ∗  3.4 Approximatedual basis
ΛyieldKin[∂g]tihsensotatnzdearordinbaDsis/DB = (ixffβiin)i=a1d,d...i,tkio.nAtnoe(lie)m, (eini)t In our deflation method, we assume that the multiple
∈ t t−1 pointisknownapproximatelyandweuseimplicitlyTaylor’s
of Th. 3.2 we impose:
expansion of the polynomials at this approximate point to
(iii) Λ[xβi]=0, 1 i k . deduce the dual basis, applying the algorithm of the previ-
≤ ≤ oussection. Tohandlesafely thenumericalproblemswhich
Proof. Let Λ K[∂]beanon-zerofunctionalsatisfying may occur, we utilize the following techniques:
(i), (ii) and (iii).∈Then Λ Dt and Λ[xβi] = 0 for i = At each step, the solutions of linear system (9, i-iii)
1,...,k. If Λ ∈ Dt−1, the∈n Λ = ki=1 λiΛi. Take for are•computed via Singular Value Decomposition. Using a
i0 the minimal i such that λi 6= 0.PThen Λ[xβi0] = λi0, given threshold, we determine the numerical rank and an
whichisincontradictionwithcondition(iii). Thus,thenon- orthogonalbasisofthesolutionsfromthelastsingularvalues
zerosolutionsof(i),(ii)and(iii)correspondtotheelements and thelast columnsof theright factor of theSVD.
which are not zero in D /D . Forthecomputationofthemonomialswhichdefinethe
t t−1 •
equations (3.4, iii) at the next step, we apply QR decom-
The above constraint is easy to realize; it is equivalent position on thetranspose of the basis to extract a non-zero
to i, dβi / suppΛζ, which implies adding a row (linear maximal minor. The monomials indexing this minor are
∀ ζ ∈
constraint) foreveryi. Inmanycases thisconstraint isjust used to determine constraints (9, i-iii). A similar numerical
λ = 0 for some i,k, thus we rather remove the column techniqueis employed in [21], for Macaulay’s method.
ik
4. DEFLATIONOF A SINGULARPOINT Theentriesofv areindexedby∂ . Thisimpliesthatanon-
i
haWs aemcounlstiidpeleraposiynstteamt xof=eqζu.atAiolsnos,fle=t (f=1,(.b.1.,,.f.s.),,bwµh)ibche tzhereofdoilffloewreinngtiarlel∆at=ionvs1:∂1(∆+Λ··i·)ζ+[fvjn]∂=n o0f,oird=er1o,n..e.s,aµt,isjfie=s
abasisofR/ and =(Λ ,...,Λ )itsBdualcounterpart, 1,...,s. By thestandard derivation rules, we have
ζ 1 µ
Q D
with Λ1 =1. d d
(∆Λ )=v Λ +∆ Λ ,
We introduce a new set of equations starting from f, as d∂ i k i d∂ i
k k
follows: add for everyf thepolynomial g =f +p ,p =
µ ε b where ε =i(ε ,...,ε ) is ak newkvectkor okf µ for i = 1,...,µ,,k = 1,...,n. Since D is stable by deriva-
i=1 i,k i k 1,k µ,k tion, d Λ D. We deducethat the vectorspace spanned
vParCiaobnlseisd.er thesystem by Dd,∂∆kDi ∈is stable by derivation and vanishes on f at ζ.
ByhProposiition 2.2, we deduce that ∆D D. This is a
⊂
g(x,ε)= Λ1(∂x)[g],...,Λµ(∂x)[g] . contradiction,since∆isofdegree1andthetheelementsin
D (cid:16) (cid:17) D are of degree N.
where Λx[gk] = Λi(∂x)[gk] is defined as in (2) with ζ re- ≤
placed by x, ie. we differentiate g but we do not evaluate ThecolumnsofJDg(x,ε)areindexedbythevariables(x,ε),
k
whiletherowsareindexedbythepolynomialsg . Wecon-
atζ. Thisisasystemofµsequations,whichweshallindex ik
struct thefollowing systems:
g(x,ε)=(g ,...,g ). Wehave
1,1 µ,s
D
gik(x,ε)=Λxi[fk+pk]=Λxi[fk]+Λxi[pk]=Λxi[fk]+pik(x,ε). (a) sLpeotnDdifnIgbneraowsusbosfysJtem(ζo)f aDrfe lsiunceharltyhaitndtehpeecnodrernet-
Df
iNsodtuicaeltthoatp.i,k(ζ,ε)=Λζi[pk]=εi,kbecauseD=(Λ1,..,Λµ) (tLheeimr.in4d.i2c)e.s.We denote by I = {(i1,k1),...,(in,kn)}
B
Asthefirstbasis element of is 1(theevaluationat the
D (b) Let g˜(x,ε˜)bethesquaresystemformedbyremoving
root), thefirst s equations are g(x,ε)=0. D
the variables ε ,...,ε from g(x,ε). There-
Notethatthissystemisunder-determined,sincethenum- i1,k1 in,kn D
ber of variables is µs+n and the number of equations is foretheJacobianJDg˜(x,ε˜)derivesfromJDg(x,ε),af-
ter purging the columns indexed by ε ,...,ε ,
µs. Weshallprovideasystematicwaytochoosenvariables and it’s (i ,k ) row becomes [ (Λxg˜ i1,k1)T 0 ].in,kn
and purge them (orbetter, set them equalto zero). j j ∇ ij ij,kj |
Thiswaywearrivetoasquaresystem g(x,ε˜)(weuseε˜ Theorem 4.3 (Deflation Theorem 1). Let f(x) be
D
for the remaining µs n variables) of size µs µs. We a n variate polynomial system with anµ foldisolated zero
shall prove that this −system vanishes on (ζ,0)×and that at x−= ζ. Then the n n system fI(x−) = 0, defined in
JDg(ζ,0)=0. (a), has a simple root a×t x=ζ. D
6
By linearity of theJacobian matrix we have Proof. By construction, ζ is a solution of fI(x) = 0.
D
JDg(x,ε)=JDf(x,ε)+JDp(x,ε) Moreover, theindices I are chosen such that detJDfI(ζ)6=
=[J (x) 0 ]+[Jx (x,ε) Jε (x,ε)], 0. This shows that ζ is a simple (thus isolated) root of the
Df | Dp | Dp system fI(x)=0.
whereJx (x,ε)(resp. Jε (x,ε))istheJacobian matrixof D
Dp Dp Example 4.4. In our running example, we expand the
p with respect to x (resp. ε).
D rectangular Jacobian matrix of 6 polynomials in (x1,x2).
aDtpL(x=e,mε(mp)1a=,1,4(.ζ..1,..0,)Tpµihs,es)tJhwaeictiohdbeipnai,ntkit(JyεDεkmp)(a=xtr,Λiεxxi)i[nopfkd]t(ihxme,elεninks)eioaenrvasµlyussa.tetemd Csdy21hs)ot[feo1ms]i,n(gwfte,h2fiexnrdo)whasacnsooarnr-ersesipgnuoglnuadlrainrrogmottionafot1r,ζanh=den((0cde,10−t)h.de22r−esdu1ldt2in−2g
1 1
Proof of Lemma 4.1. Firstnotethatthesystemisblock
The deflated system fI(x)=0 is a square system in n
separated,i.e. everypikdependsonlyonvariablesεiandnot D
on all variables ε = (ε ,...,ε ). This shows that Jε(x,ε) variables. Contrarily to the deflation approach in [7, 3], we
1 n p
is block diagonal, do not need to introduce new variables here and one step
of deflation is sufficient. In the following theorem, we do
J1 0 introduce new variables to express the condition that the
JDεp(x,ε)= ... . perturbedsystem has agiven multiplicity structure.
 0 Jµ  Theorem 4.5 (Deflation Theorem 2). Let f(x) be
Now we claim that all these blocks are all equal to the a n variate polynomial system with anµ foldisolated zero
− −
identity matrix. To see this, consider their entry ∂ [p ] at x = ζ. The square system g˜(x,ε˜) = 0, as defined in
εkj ik D
for i,j =1,...,µ, which is (b), has a regular isolated root at (x,ε˜)=(ζ,0).
Proof. By definition of , we have
d Λx [p ]=Λx[ d p ]=Λ [b ]= 1 ,i=j , D
dε i−1 k i dε k i j (cid:26) 0 ,otherwise g˜(ζ,0)=(Λζ[f],...,Λζ[f])=0.
kj jk D 1 µ
since d p = d (b ε )=b . Moreover, by construction of g˜ we get, up to a row per-
dεj,k k dεj,k j j,k j mutation, thedeterminant: D
Lemma 4.2. The µs n Jacobian matrix J (x) of the
system Df(x)=(f1,...×,fµn) is of full rank nDaft x=ζ. ±det JDg˜(ζ,0)=det(cid:12)(cid:12) JJ12 0I (cid:12)(cid:12)=det J1 6=0,
Proof of Lemma 4.2. Supposethatthematrixisrank- (cid:12) (cid:12)
deficient. Then thereis anon-trivial vectorin its kernel, where J1 = JDfI(ζ). This(cid:12) shows (cid:12)that (ζ,0) is regular
and thus isolated point of the algebraic variety defined by
J (ζ) v=0. g˜(x,ε˜)=0.
Df
· D
Nevertheless, this deflation does differ from the deflation 1.016)x +ε , .04217+2.116x +ε , ε .03921), which
2 21 2 22 23
−
strategy in [7, 3]. There, new variables are added that cor- hasaregularrootforζ Zandparameters(ε ,ε ,ε ,ε ).
12 21 22 23
respondtocoefficientsofdifferentialelements,thusintrodu- Indeed,applyingTheor∈em5.1withZ′ =Z [ .04,.04]4and
cing a perturbation in the approximate dual basis, in case × −◦
(ζ∗,0,..,0) we get an inclusion V(g,Z′,ζ∗) Z′. 2
ofapproximateinput. Hencetheoutputconcernsadeflated ⊆
root of the given approximate system. In our method, we
perturb theequations, keepingan approximatestructureof 6. GEOMETRYAROUNDASINGULARITY
a multiple point. Consequently, the certification of a root
Asafinalstepinanalyzingisolatedsingularities,weshow
concerns a nearby system, within controlled error bounds,
how the local basis can be used to compute the topological
as it shall be described in Sect.5.
degree around the singular point. If the latter is a self-
Wementionthatitwouldalsobepossibletousetheequa-
intersection point of a real algebraic curve, one can deduce
tions (9, i-iii) to construct a deflated system on the differ-
thenumberof curvebranchesthat pass through it.
entials and toperturb theapproximatedual structure.
Topological degree computation. Letf(x)beasquare
5. VERIFYINGAPPROXIMATESINGULAR n variatesystemwithanµ foldisolatedzeroatx=ζ. To
a−linear form Λ R[∂ ], we−associate thequadraticform
POINTS ∈ ζ
Q : R/ R/ R , (b ,b ) Λ(b b ) (15)
In real-life applications it is common to work with ap- Λ Q× Q→ i j 7→ i j
proximate inputs. Also, there is the need to (numerically) for R/ = b ,...,b . The signature of this (symmetric
1 µ
decide if an (approximate) system possesses a single (real) andbi-Qlinearh)formisthiesumofsignsofthediagonalentries
rootinagivendomain,notablyforuseinsubdivision-based of any diagonal matrix representation of it.
algorithms, eg. [13, 10].
In the regular case, Smale’s α theory, extending New- Proposition 6.1 ([4, Th. 1.2]). If Q , Φ is any
ton’smethod,canbeusedtoansw−erthisproblem. Another bi-linear symmetric form such that Φζ[detJΦ (x)]∈>D0, then
f
option is Rump’s Theorem, also based on Newton theory.
Inourimplementationwechoosethislatterapproach,since tdegζ(f)=sgn(QΦ). (16)
itissuitableforinexactdataandsuitsbestwiththepertur-
bation which is applied. Our perturbation coincides to the This signature is independentof thebi-linear form used.
numerical scheme of [16] in theunivariate case. We can use this result to compute the topological degree
The certification test is based on the verification method at x = ζ using the dual structure at ζ. Since a basis is
D
of Rump[16, Th. 2.1], which we rewrite in our setting: availablewesetΦ= Λi,forsomebasiselementthatisnot
±
zeroondetJ (x). Indeed,suchanelementcanberetrieved
f
Theorem 5.1 ([16] Rump’s Theorem). Let f Rn amongthebasiselements,sincedetJ / f ,see[5,Ch.0].
be a polynomial system and ζ∗ Rn a real point. G∈iven In practice it suffices to generatefa∈rahndiom element of
an interval domain Z IRn con∈taining ζ∗ Rn, and an D, compute it’s matrix representation [Φ(b b )] , and then
interval matrix M ∈IR∈n×n whose i−th colum∈n Mi satisfies extract thesignature of QΦ. i j ij
f (Z) M for i=1...,n, then the following holds:
∇ i ⊆ i Branches around a singularity. In the context of com-
If the inclusion
puting with real algebraic curves, the identification of self-
V(f,Z,ζ∗)=−Jf(ζ∗)−1f(ζ∗)+(I−Jf(ζ∗)−1M)Z ⊆Z◦ ilnoctearlsteocptioolnogpyo.inAtssaisseocnolnydthsteepfi,rsotnsetenpeeodfsdteotecramlciunliantge tthhee
(14)
number of branches attached to the singular point ζ, here-
is true, then there is a unique ζ Z with f(ζ)=0 and the after denoted Br(f,ζ). This information is encoded in the
∈
Jacobian matrix J (ζ) M is non-singular. topological degree.
f
∈
An implicit curve in n space is given by all points sati-
This theorem is applied on the perturbed system. If the sfyingf(x)=0,f =(f ,−...,f ). Considerp(x)=(x
1 n−1 1
test succeeds,wealso get adomain forε variablesthatre- ζ )2+ +(x ζ )2,andg(x)=detJ (x). Then([19−]
flectsthedistanceoftheapproximatesys−temfromaprecise a1nd ref·e·r·encesnt−hernein): (f,p)
system with thecomputedlocal structure.
Br(f,ζ)=2tdeg (f,g). (17)
ζ
Example 5.2. Westartwithanapproximatesystem: f =
1
1.071x 1.069x +1.018x2, f =1.024x 1.016x +1.058x2 This implies an algorithm for Br(f,ζ). First compute the
1− 2 1 2 1− 2 2 primal-dual structure of (f,g) at ζ and then use Prop. 6.1
andthedomain: Z =[ .01,.03] [ .03,.01]. TheJacobian
− × − to get tdeg (f,g).
atx=(0,0)evaluatesto.00652,henceitisclosetosingular. ζ
We set the tolerance equal to .04, i.e. the size of the
Example 6.2. Considertheimplicitcurvef(x,y)=0,in
domain,andweconsiderthecenteroftheboxasourappro-
ximate point,ζ∗ =(.01, .01). xy−plane,with f(x,y)=x4+2x2y2+y4+3x2y−y3,that
First we compute app−roximate multiplicity structure at looks like this . We search for the number of branches
ζ∗, = (1, d +1.00016d2 +.99542d d +1.03023d2, d touchingζ =(0,0).
a1.0s0tDa4n92dda22rd−b1a2.s0i0s01fo6rd1tdh2e−q2u1o.0ti3e5n1t4.d21T)1hae2slwatetlelrasin(d1i,cx1a2te,sx11t−)o, theWmeuclotmippliucitteygs(txr,uyc)tu=reJo(ff,(xf2+,gy)2)a=tζ1,8axnyd2−we6xar3r,iavnedtoththene
perturb upto linear terms. standard basis = (1,y,x,y2,xy,x2,y3,xy2,x2y). Among
Now apply the second deflation theorem 4.5 to get the the 9 elementsBof the dual basis, we find Φ = d3y + 38d4y +
6 6systemg=(1.018x2+1.071x +(ε 1.069)x , ε 1d2d2+3d4,havingvalueΦ0[detJ (x)]=54>0onthe
.0×2023, .01723+2.036x ,11.058x2+1 (1.01224−+ε )x +2(ε12− J8acxobyian8dexterminant. (f,g)
1 2 23 1 22−
Using and(6),wegetthe9 9matrixrepresentationof 8. REFERENCES
QΦ (15)Bwith ij th entryΦ[NF×(xβixβj)],andwecompute
tdeg (f,g)=sg−nQ =3, thusBr(f,(0,0))=6. 2
ζ Φ [1] L. Alberti,B. Mourrain, & J. Wintz. Topology and
arrangement computation of semi-algebraic planar
7. EXPERIMENTATION
curves.CAGD, 25:631–651, November2008.
Our method is developed in Maple. It uses Mourrain’s [2] M.F. Atiyah & I.G. MacDonald. Introduction to
Integration technique to compute (approximate) dual basis Commutative Algebra. Addison-Wesley,1969.
and derivetheaugmented system of Th.4.5. Then Rump’s [3] B. H.Dayton & Z. Zeng. Computingthemultiplicity
methodisusedtoverifytheroot. Macaulay’smethodisalso structurein solving polynomial systems. In
implemented for testingpurposes. ISSAC ’05,,pp.116–123, 2005. ACM.
Example 7.1. Considerthesystem[8]of3equationsin2 [4] D.Eisenbud & H.I. Levine.An algebraic formula for
variablesf =x3+x x2, f =x x2+x3, f =x2x +x x2, thedegree of a c∞ map germ. The Annals of
1 1 1 2 2 1 2 2 3 1 2 1 2
and thesingular point (0,0). Mathematics, 106(1):pp. 19–44, 1977.
Suppose that the point is given. Using 3.2 and 3.4 we [5] M. Elkadi& B. Mourrain. Introduction `a la r´esolution
derive the primal-dual pair = (1,d ,d ,d2,d d ,d2,d3+ des syst`emes d’´equations alg´ebriques, vol.59 of
D 1 2 1 1 2 2 2
d3+d2d d d2),whered3isunderlinedtoshowthatitcor- Math´ematiques et Applications. Springer,2007.
re1spon1ds2−to x132in theprim2al standard basis . The biggest [6] G. Lecerf. Quadratic newton iteration for systems
matrix used,i2ndepth4, wasofsize9 8, whBileMacaulay’s with multiplicity. Fo. Comp. Math.,2:247–293, 2002.
×
method terminates with a matrix of size 30 15. [7] A.Leykin,J. Verschelde,& ZhaoA. Newton’s method
×
To deflate the root, we construct the augmented system with deflation for isolated singularities of polynomial
f of 21 equations. The 21 2 Jacobian matrix Jx is of systems. TCS,359(1-3):111 – 122, 2006.
D × Df
rank 2 and a full-rank minor consists of the rows 4 and 5. [8] A.Leykin,J. Verschelde,& A.Zhao. Higher-order
Consequently find thesystem (d21[f1],d1d2[f1])=(3x1,2x2) deflation for polynomial systems with isolated singular
whichdeflates(0,0). Notethateventhoughbothequations solutions. vol. 146 of The IMA Volumes in Math. and
of the deflated system derive from f1, the functionals used its Appl., pp.79–97, 2008.
on f1 are computed using all initial equations. 2 [9] F.S. Macaulay. The algebraic theory of modular
systems. Cambridge Univ.Press, 1916.
Example 7.2. Let, as in [6, 8], f = 2x +2x2+2x +
1 1 1 2 [10] A.Mantzaflaris, B. Mourrain, & E. Tsigaridas.
2x2 + x2 1, f = (x +x x 1)3 x3, and f =
2x32+2x32−+10x2 +5x21+5)32 −10030−x5. − 1 3 Continuedfraction expansion of real roots of polyno-
1 2 3 3 − 1 mial systems. In Proc. of SNC ’09, pp.85–94, 2009.
Point (0,0, 1) occurs with multiplicity equal to 18, in
− [11] M. G. Marinari, T. Mora, & H.M. M¨oller. Gr¨obner
depth7. Thefinalmatrix sizewithourmethodis206 45,
× duality and multiplicities in polynomial system
while Macaulay’s method ends with a 360 165 matrix.
× solving. In Proc. of ISSAC ’95,pp.167–179, 1995.
Iftheobjectiveistodeflateasefficientlyaspossible,then
one can go step by step: First compute a basis for D and [12] B. Mourrain. Isolated points, duality and residues. J.
1
stop the process. We get the evaluation 1 and 2 first order of Pure & App. Alg.,117-118:469 – 493, 1997.
functionals,whichweapplytof . Wearriveto(1[f ],(d [13] B. Mourrain & J. P.Pavone. Subdivisionmethods for
1 1 2
d )[f ],(d +d )[f ])=(f , 4x +4x ,2+4x +2x )an−d solving polynomial equations. J. Symb. Comp.,
1 1 1 3 1 1 1 2 1 3
wecheckthattheJacobian d−eterminantis64,thuswehave 44:292–306, March 2009.
a deflated system only with a partial local structure. 2 [14] T. Ojika, S.Watanabe, & T. Mitsui. Deflation algori-
thm for multipleroots of a system of nonlinear equa-
Table 7 shows computations on the benchmark set of [3]. tions. J. of Math. An. & Appls.,96(2):463–479, 1983.
Multiplicity,matrixsizes attermination steparereported.
[15] S.R.Pope & A.Szanto. Nearest multivariatesystem
with given root multiplicities. JSC,44(6):606-625,2009.
Sys. µ Integration Macaulay
[16] S.Rump& S.Graillat. Verified error boundsfor
cmbs1 11 33 23 105 56
× × multiplerootsof systemsof nonlinearequations.Num.
cmbs2 8 21 17 60 35
× × Algs.,54:359–377, 2010.
mth191 4 10 9 30 20
× × [17] Michael Shub& SteveSmale. Complexity of Bezout’s
decker2 4 5 5 20 15
× × theorem I: Geometric aspects. J. of the AMS,
Ojika2 2 6 5 12 10
6(2):459–501, 1993.
× ×
Ojika3 4 24 9 60 35
× × [18] H.J. Stetter. Analysisof zero clusters in multivariate
KSS 16 569 69 630 252
× × polynomial systems. In Proc. of ISSAC ’96,pp.
Caprasse 4 34 13 60 35 127–136, New York,NY,USA,1996. ACM.
× ×
Cyclic 9 4 369 33 495 33
× × [19] Zbigniew Szafraniec. Topological degree and quadratic
DZ1 131 1450 524 4004 1365 forms. J. of Pure & App. Alg.,141(3):299 – 314, 1999.
× ×
DZ2 7 73 33 360 165
× × [20] X.Wu & L. Zhi. Computingthemultiplicity structu-
DZ3 5 14 8 30 21 refrom geometric involutiveform. In Proc. of ISSAC
× ×
’08,pp.325–332, 2008. ACM.
Table 1: Benchmark systems from [3].
[21] Z.Zeng. The closedness subspacemethod for
Acknowledgments. This research has received funding computingthemultiplicity structureof a polynomial
fromtheEU’s7th FrameworkProgramme[FP7/2007-2013], system.vol. 496 of Contemporaty Mathematics, pp.
Marie CurieITN SAGA,grantno [PITN-GA-2008-214584]. 347–362. AMS.Providence, RI,2009.
APPENDIX Example 6.2 (Cont’d).
Considertheimplicitcurvef(x,y)=0,inxy plane,with
We attach additional examples that did not fit page limits. −
f(x,y)=x4+2x2y2+y4+3x2y y3,
Verification Example. −
Let f = (x2x x x2,f = x x2). The verification
1 1 2 − 1 2 2 1− 2 thatlookslikethis . Wesearchforthenumberofbranches
method of[16]appliesalinearperturbationon thissystem,
touchingζ =(0,0).
but fails to certify theroot x=(0,0).
Thispointisofmultiplicity4,asthedualbasiswegetfor
We consider an approximate point ζ∗ = (.01, .002) and
we computethe approximatemultiplicity structure: d d
f(x,y)= f(x,y)= f(x,y)=0
dx dy
=(Λ ,...,Λ )=(1.0,1.0d ,1.0d +1.0d2,1.0d d +1.0d3)
D 1 4 2 1 2 1 2 2 is (1,d ,d ,d2+d2),which providesnoinformation for the
The augmented system g(x) = (Λ (f )) = (f , 2.0x x x y x y
1.0x2 1.0x ,2.0x 2.0x , 1.0x 1j.0xi2, f , 12.0x ,10.,20−.) numberof branches.
2− 1 1− 2 1− 2 2 − 2 Wecompute
has a Jacobian matrix:
Jg(ζ∗)T =(cid:20) ..0000 .0.1062 −.0.1969 22.0.0 1.0.004 02.0 00 00 (cid:21) and then thge(xm,uyl)ti=pliJc(ift,yx2s+tyru2)ct=ur1e8xoyf2(f−,g6)x3a,t ζ, and we
− − − −
arrive to
with a non-zero minor at the third and forth row. Using
thisinformation,weapplythefollowingperturbationtothe =(1,d ,d ,d2,d d ,d2,d3+ 3d4+ 1d2d2+ 3d4,
original system: D y x y x y x y 8 y 8 x y 8 x
9 3 9
f =x2x x x2+ε +ε x d d2+3d3,d2d d4 d2d2 d4),
1 1 2− 1 2 11 12 2 x y x x y− 8 y− 8 x y− 8 x
f =x x2+ε +ε x +ε x +ε x x
5 1− 2 21 22 2 23 1 24 1 2 and thestandard basis
Thusg(x1,x2,ε11,ε12,ε21,ε22,ε23,ε24),computedasbefore, =(1,y,x,y2,xy,x2,y3,xy2,x2y).
is a square system with additional equations B
Amongthe9 elements of thedualbasis, we find
f =1.0x2 2.0x x +1.0ε
2 1− 1 2 12 3 1 3
f3 =2.0x1x2−1.0x22−1.0x1 Φ=d3y+ 8d4y+ 8d2xd2y+ 8d4x,
f =2.0x 2.0x
4 1− 2 having value Φ0[detJ (x)] = 54 > 0 on the Jacobian
(f,g)
f = 2.0x +1.0ε +1.0x ε
6 − 2 22 1 24 determinant.
f7 =1.0ε23+1.0x2ε24 Using and (6), we compute the matrix NF(xβixβj) of
B
f =1.0ε multiplicationinR/ f,g ,givenattheendofthepage. Now
8 24 h i
a representation of Q (15) can be computed, by applying
Φ
NowtaketheboxZ1 =[ .03,.05] [ .04,.04] [ .01,.01]6. Φ0 on themultiplication table to get:
WeapplyTh.5.1ong,ie.−wecompu×te−V(g,Z ,×ζ∗−). Forthe
1
variable ε the interval is [ .015,.15] ( .01,.01), there- 0 0 0 0 0 0 1 0 0
21
− 6⊆ − 0 0 0 1 0 0 3/8 0 1/8
fore we don’t get an answer.  
0 0 0 0 0 0 0 1/8 0
WeshrinkalittleZ downtoZ =[ .03,.05] [ .02,.02]
[−.01,.01]6 and we ap1ply again T2h. 5−.1, which×re−sults in × QΦ= 00 01 00 30/8 1/08 10/8 00 00 00 .
 0 0 0 1/8 0 3/8 0 0 0 
[ .004,.004]  
−  1 3/8 0 0 0 0 0 0 0 
 [−.004,.004]   0 0 1/8 0 0 0 0 0 0 
[ .001,.001]  0 1/8 0 0 0 0 0 0 0 
V(g,Z2,ζ∗)= [[[−−−...000000967,,,...000000967]]] ⊆Z◦2, eigWenivtahluaesQoRf tihteisrarteiporne,sewnetafitniodn,6hpenocsietiwvee caonmdp3utneegative
 [−.−000[.405,,.0.0]0035]  tdegζ(f,g)=sgnQΦ=6−3=3,
 
i.e. thereare 6 branchesof thecurvearound (0,0).
thuswecertifythemultiplerootoftheoriginalsysteminside
Z .
2
1 y x y2 xy x2 y3 xy2 x2y
y y2 xy y3 xy2 x2y 3y3− 9x2y 0 1y3− 3x2y
 8 8 8 8 
x xy x2 xy2 x2y 3xy2 0 1/8y3−3/8x2y 0
 y2 y3 xy2 3y3− 9x2y 0 1y3− 3x2y 0 0 0 
 xy xy2 x2y 8 08 1y3− 3x2y 8 08 0 0 0 
 x2 x2y 3xy2 1y3−3/8x2y 8 08 3/8y3− 9x2y 0 0 0 
 xyy32 3/8y3−0 89x2y 81y3−038x2y 8 00 00 00 8 00 00 00 
 x2y 1/8y3− 3x2y 0 0 0 0 0 0 0 
 8 
Multiplication table for R/ f,g (Example 6.2).
h i