A Companion Curve Tracing Method for Rank-deficient Polynomial Systems
AA COMPANION CURVE TRACING METHOD FORRANK-DEFICIENT POLYNOMIAL SYSTEMS
WENYUAN WU ∗ AND
CHANGBO CHEN ∗† Abstract.
We propose a method for tracing implicit real algebraic curves defined by polynomialswith rank-deficient Jacobians. For a given curve f − (0), it first utilizes a regularization techniqueto compute at least one witness point per connected component of the curve. We improve thisstep by establishing a sufficient condition for testing the emptiness of f − (0). We also analyzethe convergence rate and carry out an error analysis for refining the witness points. The witnesspoints are obtained by computing the minimum distance of a random point to a smooth manifoldembedding the curve while at the same time penalizing the residual of f at the local minima. Totrace the curve starting from these witness points, we prove that if one drags the random point alonga trajectory inside a tubular neighborhood of the embedded manifold of the curve, the projection ofthe trajectory on the manifold is unique and can be computed by numerical continuation. We thenshow how to choose such a trajectory to approximate the curve by computing eigenvectors of certainmatrices. Effectiveness of the method is illustrated by examples. Key words. rank-deficiency, real algebraic curve, curve tracing, numerical continuation, penaltyfunction, Tikhonov regularization
AMS subject classifications.
1. Introduction.
Given an implicit real algebraic curve in R n defined by a finiteset of polynomials f ⊂ R [ x , . . . , x n ], producing a polygonal chain approximation ofit is a classical problem. Existing numerical continuation methods [1] for solving thisproblem often require that the Jacobian of f , denoted by J f , is of nullity one at all(or almost all) points of the real zero set of f , that is f − (0). It remains a challengeto trace the curve defined by f when J f is rank-deficient, that is of nullity greaterthan one, which is the dimension of the curve. One typical example is when f is asum of squares of polynomials or more generally when f is nonnegtive.For a smooth curve in R n defined by f = { f , . . . , f n − } , where J f is of fullrank at any point of the curve, the main technical challenge would be to identifyall branches of f − (0) and make sure that there is no jumping during curve tracing.Identifying all branches of f − (0) is the problem of computing witness points for everyconnected component of f − (0) [22, 12, 25]. Techniques for presenting or detectingcurve jumping also exist [3, 2, 19, 27, 26].For an almost smooth curve in R n defined by f = { f , . . . , f n − } , where J f is offull rank at nearly all points of the curve, one further difficulty is to trace across thesingular points and get the correct topology around the singular points. Many workexist for handling such problems [13, 8, 6, 11, 5, 14].In this paper, we are interested in the case that J f is rank-deficient at every pointof f − (0), such as when f consists of a polynomial in sum of squares. When J f isrank-deficient, one of the main obstacles is that it is hard to find a tracing directionsince the dimension of the nullspace of J f is at least two. Our starting point forsolving this problem is a penalty function based method for computing witness pointsfor every connected component of f − (0) [24]. We then extend it to a companioncurve tracing method. The main idea of the penalty function is to embed the curve ina high-dimensional smooth manifold defined by a system g with full rank Jacobians ∗ Chongqing Institute of Green and Intelligent Technology, Chinese Academy of Sciences,Chongqing, China ([email protected], [email protected], ). † Corresponding author. 1 a r X i v : . [ m a t h . NA ] J a n W. WU AND C. CHEN and compute the minimum distance of a random point to the manifold while at thesame time penalizing the residual of f at the local minima. To render the curve,one natural idea is to move the random point (as a guiding point) along a trajectoryand hopefully the corresponding minima will be dragged continuously. To implementthis idea, there are two main challenges to overcome. One is the potential occurenceof discontinuity, which indeed may happen as illustrated in Section 3 and provedin Section 5. Another is to make sure the minima do move along the curve as onedrags the guiding point along the trajectory. We show that both challenges can beovercome and the idea of moving the guiding point is indeed feasible if the point ismoved along the directions defined by eigenvectors of certain matrices and inside atubular neighborhood of the embedded manifold of the curve. The trajectory formedby moving the random point is called a companion curve of f − (0).Interestingly, the penalty function method, although initially proposed in a com-pletely different context, turns out to be closely related to Tikhonov regularization forsolving rank-deficient nonlinear least squares problems [23, 9, 10], which is explainedin the preliminary section.The paper is structured as follows. In Section 2, we recall how to compute atleast one witness point for every connected component of an arbitrary real algebraicvariety, whose defining system allows to be rank-deficient, via the so-called penaltyfunction method [24]. In Section 3, we illustrate by a simple example the challengefor tracing real algebraic curves defined by rank-deficient systems, with the initialwitness points provided by the penalty function method, as well as the main ideaof our companion curve method for handling this problem. One weakness of thepenalty function method for computing witness points is that it always return a non-empty set of points even the given real variety is empty. In Section 4, we provide asufficient criterion for testing emptiness of a variety by the penalty function method.To successfully tracing a real algebraic curve, it is important to make the approximatewitness points close enough to the curve. In Section 5, we propose a homotopy methodfor improving the precision of witness points. With all these preparations and anothertool from differential geometry, namely Tubular Neighborhood Theorem, we proposea companion curve method for curve tracing in Section 6. The effectiveness of themethod is illustrated by several examples in Section 7. Finally, in Section 8, we drawthe conclusion and propose several ways to improve the current method.
2. Preliminaries.
In this section, we recall some preliminary results that wereintroduced in [24] for computing witness points of rank-deficient polynomial systems.Let x = ( x , . . . , x n ) and let f = { f , ..., f k } ⊂ R [ x ]. Let V R ( f ) be the zero set of f in R n . Let a = ( a , ..., a n ) / ∈ V R ( f ) be a point in x -space and consider the minimaldistance from V R ( f ) to this point:min n (cid:88) i =1 ( x i − a i ) (2.1) s.t. f ( x ) = 0 . Clearly every semi-algebraically connected component of V R ( f ) has at least one pointattaining the local minimum. These points make the following matrix lose full rank: A := ∂f /∂x · · · ∂f /∂x n ... . . . ... ∂f k /∂x · · · ∂f k /∂x n x − a · · · x n − a n . OMPANION CURVE k rows of A are exactly the the Jacobian matrix of f w.r.t. x ,denoted by J f . If J f is rank-deficient at every point of V R ( f ), the matrix A auto-matically loses full rank and thus does not provide any extra helpful information oncomputing the semi-algebraically connected components of V R ( f ).To overcome this algebraic rank deficiency problem, the paper [24] introduces apenalty function based approach. Instead of solving the optimization problem (2.1),one considers the following unconstrained optimization problem:min µ = ( β · ( f + · · · + f k ) + n (cid:88) i =1 ( x i − a i ) ) / . (2.2)Note that as β approaches infinity, f i , i = 1 , . . . , k are forced to be zero. Intuitively,this provides an approximate solution to the problem (2.1) for large enough β . Indeed,the paper [24] proves the following result justifying such intuitions. Proposition p be a local minimum of (2.1). Thereexists a local minimum p (cid:48) of (2.2) for sufficiently large β , such that (cid:107) p − p (cid:48) (cid:107) can bearbitrarily small.Moreover, one can get a rough estimation of the distance between p and p (cid:48) via thenotion of degree index . Definition v (cid:54) = 0 ∈ R n , let f v := f ( x = vt ).Denote by deg min ( f v ) the trailing degree of f v . We call deg ind ( f ) = max v deg min ( f v )the degree index of f . Given a point p ∈ V R ( f ), the degree index of f at p is definedas deg ind ( f ( x + p )). Theorem
For a random point a ∈ R n and a sufficientlylarge β , suppose that p ∈ V R ( f ) attains the local minimal distance to a . Then thereis a solution p (cid:48) of Equation (2.3) such that (cid:107) p (cid:48) − p (cid:107) = O ( I − (cid:112) /β ) , where I =max { deg ind ( f i ( x + p )) , i = 1 , ..., k } . The local minima of (2.2) are exactly the points that vanish the gradient of µ ,that is satisfying the following equation:(2.3) x ... x n + β · J t · f ... f k = a ... a n . where the n × k matrix J t is the transpose of the Jacobian of f .The left hand side of Equation (2.3) defines a smooth mapping M : R n → R n . Lemma
For almost all points a = ( a , ..., a n ) / ∈ V R ( f ) , M − ( a ) is a nonempty finite set and every point of M − ( a ) is a regular point of M . Proposition 2.1 and Lemma 2.4 together show that for almost all points a of R n , M − ( a ) contains points meeting every semi-algebraically connected componentof V R ( f ). We warn that M − ( a ) may contain extra points not belonging to V R ( f ). Inparticular, even V R ( f ) is empty, M − ( a ) is always nonempty. Numerically testing ifa real variety is empty in general is a difficult problem. We provide a partial answerto this problem in Section 4.Sometimes, it is useful to use the following two equivalent formulations to (2.2).Let z = ( z , .., z k ) be k slack variables and g = { f + z , f + z , ..., f k + z k } . Note W. WU AND C. CHEN that we have g ⊂ R [ x, z ] and V R ( g ) ⊆ R n + k .min µ = ( β · ( z + · · · + z k ) + n (cid:88) i =1 ( x i − a i ) ) / s.t. g ( x, z ) = 0 . Let z i = w i / √ β , i = 1 , . . . , k , and substitute them into (2.4). Let h = { f + w / √ β, . . . , f k + w k / √ β } .min ( w + · · · + w k + n (cid:88) i =1 ( x i − a i ) ) / s.t. h = 0 . One nice thing about (2.4) and (2.5) is that the real varieties defined by g and h are smooth submanifolds of R n + k . Remark β = 1 /t , we obtain another equivalent formulation:min µ = (( f + · · · + f k ) + t n (cid:88) i =1 ( x i − a i ) ) / . (2.6)Such a formulation is exactly the Tikhonov regularization for nonlinear least squaresproblems [23, 9, 10], where t (cid:80) ni =1 ( x i − a i ) ) / (cid:107) f (cid:107) . Here, one must be cautious about the choice of a ,since Lemma 2.4 does not exclude the possibility that the solution of problem (2.6) isnot regular for certain a , which indeed poses a challenge for curve tracing as explainedin next section.
3. An introductory example.
In this section, we illustrate by an example thechallenge of a pure numerical method for tracing algebraic curves defined by rank-deficient polynomials as well as the main idea of our companion curve method.Let f := x − x x + x = ( x − x ) . Recall from Section 2 that, one can ob-tain approximate witness points of V R ( f ) by solving Equation (2.3). Choosing β = 10 and a = (0 , − { x − x x + 180000 x x − x x + x , − x + 60000 x x − x x + 20000 x + x + 1 } . Bythe homotopy continuation method, say the one implemented in Hom4PS-2.0 [18], oneobtains three approximate witness points: ( − . , − . − . , − . , − . x , x )-space of the following three local minima of the optimization problem (2.5):( − . , − . , − . , ( − . , − . , − . , (0 , − . , − . . These three points attain the local minimum distance from the point (0 , − ,
0) to themanifold f + w/
100 = 0, as illustrated by Figure 1.
Next we would like to generate more points of V R ( f ) bycurve tracing with the three initial points. One natural idea is to consider the followinglinear homotopy:(3.1) H τ ( x, τ ) = x + β J t · f − ( τ p + (1 − τ ) p ) ≡ , OMPANION CURVE Fig. 1 . (Color online) Left: the random point a = (0 , − , (in red • ), the three local minima(in blue + ), the smooth surface defined by f + w/
100 = 0 and its intersection with w = 0 (goldcurve). Right: the random point a = (0 , − (in red • ), the three approximate witness points (inblue + ), and the curve f = 0 (gold curve). where β = 10000 and one moves a from p = (0 , −
1) to p = ( − . , − .
5) in aline. As long as the number of solutions of H ( x, a ) = x + β J t · f − a in x remainsunchanged and the graphs of the solutions of H ( x, a ) as functions of a remain disjointand smooth, one should encounter no much difficulty during curve tracing. However,as illustrated by the left subfigure of Fig. 2, the path of a crosses the discriminantlocus [16] of H ( x, a ), which can be obtained by a Gr¨obner basis computation [16]by Lemma 5.1 in Section 5. As illustrated by the right subfigure of Fig. 2, when a approaches the discriminant locus, the solution curve starting with (0 , − . − . , − . • ) get close to eachother gradually and finally are no longer solutions of H ( x, a ). Interestingly, if we applyNewton iteration now starting with these non-solution points, they finally converge tothe third solution curve traced starting with the local minimum ( − . , − . (cid:63) ). Therefore, curve jumping happened when the parameter a crossed thediscriminant locus. Note that in Fig. 1, there are threepoints on the smooth surface attaining the local minimum distance to the given ran-dom point (0 , − , a = p is the initial random point and ( x, w ) = ( q , r ) is anapproximate local minimum of the optimization problem (2.5). Suppose that (cid:107) r (cid:107) issmall enough such that x = q can be seen as an approximate witness point of V R ( f ).We first move the point a to another point p on the segment p q such that ( p ,
0) isclose to the surface and hope that there is only one point ( q , r ) on the surface withlocal minimum distance to ( p ,
0) nearby. We then pick a well chosen direction, say (cid:126)v ,by some eigenvectors computation and move a from p to p (cid:48) . Accordingly, x is moved W. WU AND C. CHEN
Fig. 2 . (Color online) Left: the discriminant locus of H ( x, a ) (two “V” curves in darkgreen) and the segment p p (in red). Right: the traced curve starting with the initial point ( − . , − . (gold (cid:63) ), the traced curve starting with the initial point (0 , − . (blue + ),and the traced curve starting with the initial point ( − . , − . (purple • ). Fig. 3 . (Color online) Left: illustrating the idea of companion curve tracing method. Right:an approximation of V R ( f ) (polygonal chain in gold • ) and its companion curve (polygonal chain inblue + ). The initial random point (0 , − (in red • ) at the bottom was moved to a point (in blue (cid:3) )close to V R ( f ) . in the same direction from q to q (cid:48) (may be further refined to q (cid:48)(cid:48) by Newton iterationif necessary). It is possible that ( p (cid:48) ,
0) is now outside the tubular neighborhood ofthe surface and one may repeat the previous step to drag p (cid:48) to p and produce q .And then move q to q (cid:48) , p to p (cid:48) , so on so forth. The polygonal chain formed by thesequence of points p , p , . . . is called the companion curve of V R ( f ), as illustrated inFig. 3.
4. Emptiness of real variety.
In this section, we propose a criterion for testingthe emptiness of a given real variety.Lemma 2.4 implies that all the solutions of Equation (2.3) can be obtained by
OMPANION CURVE (cid:107) z (cid:107) (cid:28)
1. It is possible that such points do not exist,which then provides strong evidence that V R ( f ) is empty. Intuitively, this is because if V R ( f ) is not empty, increasing the penalty factor β will force (cid:107) z (cid:107) close to zero. Thus,the minimal value of µ will be slightly larger than the distance from a to V R ( f ).To study the relationship between µ min and the emptiness of V R ( f ), we ho-mogenize the system f by adding a variable x satisfying a new equation ¯ f k +1 = (cid:80) ni =0 x i − homogenized system ¯ f except for the new inhomogeneousequation. The corresponding unconstrained optimization problem ismin ¯ µ = ( β · ( ¯ f + · · · + ¯ f k + ¯ f k +1 ) + n (cid:88) i =0 ( x i − a i ) ) / . (4.1)where a is chosen randomly in the unit ball B (0; 1) of R n +1 . Proposition µ min be the global minimal value of ¯ µ in the optimizationproblem (4.1). If ¯ µ min >
2, then V R ( f ) = ∅ . Proof.
We prove it by contradiction. Suppose V R ( f ) (cid:54) = ∅ . Then V R ( ¯ f ) (cid:54) = ∅ andlet p ∈ V R ( ¯ f ). We know that p, a ∈ B (0; 1) which implies (cid:107) p − a (cid:107) ≤
2. Thus,¯ µ min ≤ β · (cid:107) p − a (cid:107) / ≤
2. It contradicts the assumption ¯ µ min > (cid:3) Example f = x + y − xy + 3 /
2. We can verify that f = 159784 x + 1564 y + 350 + (cid:18) − xy (cid:19) + (cid:18) x − y (cid:19) > . Homogenizing f yields ¯ f = { x + y − h xy + 3 / h , h + x + y − } . Choose a = (0 . , . , .
3) and β = 10000 and solve the corresponding system of Equation(2.3) by Hom4Ps2 . It gives 23 real roots among 111 complex ones. The minimal valueof ¯ µ is 28 . x = 0 . , y = 0 . , h = 0 . V R ( f ) = ∅ .To apply this lemma we may choose sufficiently large β to make ¯ µ min > β is. Example f = ( xy − + y and it is positive but arbitrarily closeto zero. Fixing a = (0 , , µ = β (cid:0) h − h xy + h y + x y (cid:1) + β (cid:0) h + x + y − (cid:1) + h + x + y . Numerical computation shows that ¯ µ min = 0 . x = 2 . × − , y =0 . , h = 8 . × − ) when β = 10 . Actually, ¯ µ min must be no greater than 1 since¯ µ (0 , ,
0) = 1 for any β . It means that we cannot tell if the real variety of f is emptyor not by Proposition 4.1.Proposition 4.1 only gives a sufficient condition for V R ( f ) = ∅ . If ¯ µ min < V R ( f ) (cid:54) = ∅ .
5. Refinement.
By Proposition 2.1, theoretically we can use Equation (2.3) toupdate the approximate root x (cid:48) by increasing β . Suppose we have all the real solutionsof Equation (2.3) for β = β denoted by R β . When we increase β from β to β , W. WU AND C. CHEN there are two ways to obtain R β . One way is to solve Equation (2.3) in the complexfield and then keep the real solutions. Alternatively, we may trace the real curves ofthe following homotopy starting from all points of R β by moving β from β to β continuously.(5.1) H ( x, β ) = ( x − a ) + β · J t · f ≡ . But we have to be aware of singular Jacobian for a successful tracing.
Lemma
For any polynomial system f ⊂ R [ x ] , let F = { x − a + β · J t · f } and G = { F, det( ∂F∂x ) } ⊂ R [ x, a , β ] . Then there is a nonzero polynomial φ ( a , β ) ∈ (cid:104) G (cid:105) . Proof.
When (cid:104) G (cid:105) = (cid:104) (cid:105) , it is trivial.Otherwise, let Gb be a Gr¨obner basis of G with respect to lex order x (cid:31) a (cid:31) β . Then Gb (cid:48) = Gb ∩ R [ a , β ] is a Gr¨obner basis of the elimination ideal. If Gb (cid:48) = 0, then fora generic ( a , β ) ∈ C n +1 , it can be extended to a solution z of G = 0. But by Sard’slemma for varieties (Chap. 3 in [21]), this point is a regular value of the smoothmapping with nonsingular Jacobian. Thus, z does not satisfies det( ∂F∂x ) = 0 whichcontradicts the assumption z is a solution of G = 0.Therefore, Gb (cid:48) must contain nonzero polynomials. (cid:3) This lemma shows that even a is chosen randomly, it is still possible to encountersingular Jacobian when β increases continuously from β to infinity. Example f := ( x − x ) in Section 3. If we choose a = (0 , − (cid:26) β x − β x y + 18 β x y − β x y − a1 + x = 0 − β x + 6 β x y − β x y + 2 β y − a2 + y } = 0 . Applying Lemma 5.1, we obtain φ ( a , β ) = 2082930190011 β − β − β − β − β − , which hasonly one positive real root β ≈ .
28. So if we choose β ≤
23, singular Jacobian willbe encountered when increasing β from β to infinity.Next we will study such probability. After we fix the value of a , Gb (cid:48) is a set ofunivariate polynomials and it is generated by a single polynomial denoted by g ∈ R [ β ].By rescaling, the all coefficients of g are in [ − , Lemma
Let g = (cid:80) ni =0 c i β i ∈ R [ β ] with degree n . Suppose the coefficients { c , ..., c n } are i.i.d. random variables with uniform distribution in [ − , . Let k bethe largest absolute value of all the roots of g . Then the probability Pr( k < N + 1) > − /N . Proof.
By the root bound in Chap. 8 of [20], we have k < | c n | + max( | c n − | , ..., | c | ) | c n | . Then clearly, it implies Pr( k < N + 1) > Pr(max( | c n − | , ..., | c | ) ≤ | c n | N ).Let c = max( | c n − | , ..., | c | ) which is a random variable in [0 , c ≤ t ) = t n if t <
1, and Pr( c ≤ t ) = 1 if t ≥
1. HencePr( c ≤ | c n | N ) = Pr( | c n | ≥ /N ) + Pr( | c n | < /N ∧ c ≤ | c n | N )= 1 − /N + (cid:90) /N ( | c n | N ) n d | c n | = 1 − /N + 1 / ( n + 1) /N > − /N. OMPANION CURVE k < N + 1) > − /N . (cid:3) This lemma indicates that the probability of singularity occurring during the ho-motopy path tracking (5.1) is less than 1 / ( β − β e.g. β = 10 .Proposition 2.1 shows the existence of an approximate point p (cid:48) for any criticalpoint p by solving Equation (2.3). With the distribution assumption in Lemma 5.3,we can also show the uniqueness of p (cid:48) with a quite high probability. Corollary a ∈ R n and a sufficiently large β , suppose p ∈ V R ( f ) attains the local minimal distance to a . Then there is a unique solution p (cid:48) ofEquation (2.3) such that (cid:107) p (cid:48) − p (cid:107) = O ( I − (cid:112) /β ), where I = max { deg ind ( f i ( x + p )) , i =1 , ..., k } with probability at least 1 − / ( β − Proof.
Compared with Theorem 2.3, we only need to show the uniqueness of p (cid:48) .Suppose there is another point p (cid:48)(cid:48) also close to p such that (cid:107) p (cid:48)(cid:48) − p (cid:107) = O ( I − (cid:112) /β ).By Lemma 5.3, when β → ∞ , with probability at least 1 − / ( β − p (cid:48) and p (cid:48)(cid:48) approach p as β moves to ∞ continuously. And it means that for a generic a , φ ( a , β ) → β → ∞ , where φ is a nonzero polynomial given in Lemma 5.1.This indicates that g ( a ) = 0 where g is the leading coefficient of φ with respect to β . But a is generic, g ( a ) (cid:54) = 0. It leads to a contradiction. Therefore, p (cid:48) is unique forsufficiently large β . (cid:3) Replacing β with 1 /t in (5.1), we obtain a homotopy(5.2) H ( x, t ) := t ( x − a ) + J t · f = 0 . Let ( t , x ) be an initial point satisfying H ( x , t ) = 0. As t approaches zero, withthe assumption in Lemma 5.3, the homotopy path x ( t ) will approach V R ( f ) in highprobability (about (1 − t )). As a direct consequence of Theorem 2.3, we have thefollowing error estimation for the refinement process. Corollary
Let d be the degree of f . If we reduce the value of t by a halfat each step and assume that the step sizes t s , s = 1 , , . . . are small enough to avoidcurve jumping, then after s steps of path tracking, the error of root is reduced to O ( τ s δ ) , where τ = 2 − / (2 d − , and δ is the initial error (cid:107) p (cid:48) − p (cid:107) .Remark g := { f + z , . . . , f k + z k } as a perturbed system f (cid:48) of the input f when fixingthe values of z i . More precisely, a critical point p (cid:48) which is a solution of Equation (2.3)is an exact solution of f (cid:48) = { f + z , ..., f k + z k } where z i = − f i ( p (cid:48) ). By Theorem 2.3,we have the following result. Corollary p (cid:48) be a local minimum of µ in (2.4) with a sufficiently large β .Then the residual error satisfies(5.3) (cid:88) z i = O (1 /β I +12 I − ) . Proof.
Since µ ( p (cid:48) ) is a local minimal value, then the perturbations z i are very smallbecause µ ( p (cid:48) ) ≤ µ ( p ) ⇒ β ( (cid:88) z i ) + (cid:107) p (cid:48) − a (cid:107) ≤ (cid:107) p − a (cid:107) ⇒ β ( (cid:88) z i ) ≤ (cid:107) p (cid:48) − p (cid:107) W. WU AND C. CHEN where p is a solution of Problem (2.1) and the forward error (cid:107) p (cid:48) − p (cid:107) = O ( I − (cid:112) /β )when β is sufficiently large by Theorem 2.3. It implies that the backward error (cid:80) z i = O (1 /β I +12 I − ). (cid:3)
6. Tracing algebraic curves defined by rank-deficient systems.
Previ-ously we have presented a method to approximate and refine the real witness pointsof a general system combining critical point and penalty function techniques. Furtherapplications may require to produce more points on the real variety, especially whenthe variety is an algebraic curve. It is quite challenging for rank-deficient systems,since the Jacobian is always singular along the whole curve. A straightforward methodfor generating more points by moving the random point a to generate more criticalpoints may fail as illustrated in Section 3.Thanks to one of the fundamental theorems in differential geometry: A smoothmanifold is contained in an open tubular neighborhood in which every point can beuniquely projected onto the manifold following a normal line, we can track the criticalpoint along the curve if a is always contained in this tubular neighborhood. First we give a brief review of somerelated concepts from differential geometry. For each x ∈ R n , the tangent space T x R n is canonically identified with R n and the tangent bundle T R n is canonicallydiffeomorphic to R n × R n . Let M ⊆ R n be an embedded k -dimensional submanifold.For each x ∈ M , the normal space to M at x is defined to be the ( n − k )-dimensionalsubspace N x M ⊆ T x R n consisting of all vectors orthogonal to T x M with respect tothe Euclidean inner product. The normal bundle N M := { ( x, v ) | x ∈ M, v ∈ N x M } of M is an embedded n -dimensional submanifold of T R n by Theorem 6.23 in [17].We define a smooth map E : N M → R n by E ( x, v ) = x + v where x ∈ M and v ∈ N x M . A tubular neighborhood of M is a neighborhood U of M in R n which is the diffeomorphic image under E of an open subset V ⊆ N M of the form V = { ( x, v ) ∈ N M : (cid:107) v (cid:107) < δ ( x ) } for some positive continuous function δ : M → R . Theorem
Every embedded subman-ifold of R n has a tubular neighborhood. Tubular neighborhood theorem is also true for complex analytic manifold (See forinstance Theorem 6.2 in [28]).
Lemma
Let h = { h , ..., h k } ∈ R [ x ] with V R ( h ) (cid:54) = ∅ , where k < n . Supposethat the Jacobian matrix J of h attains full rank at any point of V R ( h ) . Then thereexists a tubular neighborhood U of V R ( h ) such that every point b ∈ U has a uniqueprojection x b ∈ V R ( h ) of minimum distance to b . Moreover, x b is the projection of aunique isolated simple zero of the Lagrangian system of min ( (cid:107) x − b (cid:107) ) / s.t. h ( x ) = 0 . Proof.
The Lagrangian system of (6.1) is as below: x − b ... x n − b n = ∂h /∂x · · · ∂h k /∂x ... . . . ... ∂h /∂x n · · · ∂h k /∂x n n × k · λ ... λ k (6.2) h ( x ) = 0 . Since the Jacobian matrix of h is of full rank at any point of V R ( h ), the set V R ( h )is a submanifold of R n . Theorem 6.1 guarantees that there is a tubular neighborhood OMPANION CURVE V R ( h ) such that b ∈ U . Next we show that x b is unique. Suppose that there isanother different point x (cid:48) b ∈ V R ( h ) of minimum distance to b . Let v (cid:48) = b − x (cid:48) b and v = b − x b . Since J t is of constant rank k , its column vectors form a basis of thenormal space of dimension k . So v, v (cid:48) are normal vectors of x b , x (cid:48) b respectively and b = E ( x b , v ) = E ( x (cid:48) b , v (cid:48) ). Since b ∈ U , E is injective and it implies ( x b , v ) = ( x (cid:48) b , v (cid:48) ).It is a contradiction.On the other hand, because the Jacobian matrix of h is of full rank at any pointof V R ( h ), the value of λ = ( λ , . . . , λ k ) is uniquely defined by a zero x of h ( x ). So x b is the projection of a unique zero ( x b , λ b ) of Eqs. (6.2).Next we show that ( x b , λ b ) is an isolated simple zero of Eqs. (6.2). Otherwise,there must exist another point x (cid:48) b ∈ V C ( h ) near x b also satisfying (6.2), which isimpossible by applying the Tubular Neighborhood Theorem for the complex analyticmanifold V C ( h ) following a similar argument as we did for the real case. (cid:3) Corollary f = { f , ..., f k } be a set of polynomials in the ring R [ x ] and h = { f + w / √ β, ..., f k + w k / √ β } where β is a positive constant. Then V R ( h ) hasa tubular neighborhood U ⊂ R n + k . In addition, for any point ( a , b ) ∈ U , thereis a unique point ( x ∗ , w ∗ ) ∈ V R ( h ) to minimize the distance from ( a , b ) to V R ( h )and ( x ∗ , w ∗ ) is the projection of a unique isolated simple zero of the correspondingLagrangian system of min ( (cid:107) w − b (cid:107) + (cid:107) x − a (cid:107) ) / s.t. h = 0 . Proof.
Note that h attains full rank at any point of V R ( h ). Then the conclusionfollows directly from Lemma 6.2. (cid:3) Theorem
Let f = { f , ..., f k } be a set of polynomials in the ring R [ x , ..., x n ] with V R ( f ) (cid:54) = ∅ . For any β > , there exists an open set in R n containing V R ( f ) . Andfor any point a in this set there is a unique point on V R ( h ) to minimize the distancefrom ( a , ∈ R n + k to V R ( h ) , where h = { f + w / √ β, ..., f k + w k / √ β } .Moreover, if a moves along some piecewise smooth curve C in this set, the corre-sponding projection trajectory can be obtained by solving Equation (2.3) continuouslyfrom an initial point a ∈ C . Proof.
By Corollary 6.3, V R ( h ) has a tubular neighborhood U containing V R ( h ) in R n + k such that for any point ( a , b ) ∈ U , there is a unique point ( x ∗ , w ∗ ) ∈ V R ( h ) tominimize the distance from ( a , b ) to V R ( h ) and ( x ∗ , w ∗ ) is the projection of a uniqueisolated simple zero of the corresponding Lagrangian system of Equation (6.3):(6.4) x − a ... x n − a n w − b ... w k − b k = ∂f /∂x · · · ∂f k /∂x ... . . . ... ∂f /∂x n · · · ∂f k /∂x n / √ β . . . 1 / √ β ( n + k ) × k · λ ... λ k . with h = 0. Since ∅ (cid:54) = V R ( f ) = π x ( V R ( h ) ∩ V R ( w )) ⊂ π x ( U ∩ V R ( w )) and U is an openset, letting U (cid:48) = π x ( U ∩ V R ( w )), we know that U (cid:48) is an open set in R n containing V R ( f ). For any point a ∈ U (cid:48) , we have ( a , ⊂ U . Thus, by Corollary 6.3 there isa unique point ( x ∗ , w ∗ ) ∈ V R ( h ) to minimize the distance from ( a ,
0) to V R ( h ) and2 W. WU AND C. CHEN ( x ∗ , w ∗ ) is the projection of a unique isolated simple zero of Equation (6.4). Sincethere is a one-to-one correspondence between the zeros of Equation (6.4) and thoseof Equation (2.3), x ∗ is an isolated simple zero of Equation (2.3). Therefore, theprojection trajectory can be obtained by solving Equation (2.3) continuously from aninitial point a ∈ C . (cid:3) A natural question is how to move a into the tubular neighborhood. Recall theequation(6.5) x + β · J t · f = a . Suppose x , a and β satisfy this equation. So λ ( x − a ) + λβ · J t · f = 0 for any λ >
0. Let a = (1 − λ ) x + λ a and β = λβ . Then ( x − a ) + β · J t · f = 0.So x is a critical point of Equation (6.5) with a = a and β = β . When λ issmall, a is very close to x . By Theorem 2.3, x should be quite close to the curvefor sufficiently large β and consequently a is also close to the curve. Since β = λβ ,we can increase the value of β by the homotopy (5.2) and it gives a refined solution x even closer to the curve.In summary, we have the following algorithm. Algorithm 6.1
Move towards Tubular Neighborhood
Input:
A system f = { f , . . . , f k } ∈ R [ x , . . . , x n ]. A contraction factor 0 < λ < x , a ) satisfying Equation (6.5) with a fixed penalty factor β (cid:29) Output:
A point ( x , a ) satisfying Equation (6.5), where β = β , with a = (1 − λ ) x + λ a . begin let a = (1 − λ ) x + λ a and β = λβ let t = 1 /β and t = 1 /β construct H ( x, t ) = t ( x − a ) + J t · f ≡ x , t ) move the parameter t from t to t continuously yields H ( x , t ) = 0 return ( x , a ). end Remark a belongs to thetubular neighborhood of V R ( f ) since there is no information about the radius of thetubular neighborhood and the contraction factor is just chosen by the user.However, the tubular neighborhood condition is sufficient but unnecessary forregular Jacobian of Equation (6.5). Therefore, we can monitor the Jacobian duringchanging a along some direction. If the Jacobian is close to being singular, it indicatesthat a might be out of the tubular neighborhood and we can use Algorithm 6.1 tomove a towards the tubular neighborhood. Since the tubular neighborhood is an openset, a will be in the neighborhood after finitely many iterations. Thus, this restoresthe regularity of Equation (6.5). The next question to answer is how to choose thedirection to move a in Equation (6.5). By our assumption, V R ( f ) is a one dimensionalcurve and the best choice of moving a is along the tangent direction of the curve.However it will be difficult to obtain such a direction if Jacobian of f is rank-deficientalong the whole curve. Here we present a method for finding the tracing direction bycomputing eigenvalues. OMPANION CURVE a and x . That is there exists a w such that ( x, w ) is the projection of ( a ,
0) onto V R ( h ), where h = f + w/ √ β . When β is sufficiently large, by Theorem 2.3 we can consider x as the projection of a onto V R ( f ) approximately. If a is quite close to the curve V R ( f ), roughly speaking whenthe change of a is parallel to the tangent at x , then we should have ∆ a ≈ c ∆ x forsome constant c close to 1. This implies that(6.6) ( β ∂ J t · f∂x + I )∆ x = ∆ a ≈ c ∆ x. Therefore, ∆ a will be approximately equal to the eigenvector corresponding to theeigenvalue 1.Since V R ( f ) is of dimension one and this eigenvector is close to the tangent, othereigenvectors will approximately span the normal space of the curve at this point.Intuitively, a large change of a in the approximate normal space will not change theprojection x too much and consequently the eigenvalues of such eigenvectors will bevery large. Proposition β ∂ J t · f∂x + I is symmetric and it has n orthogonaleigenvectors. Proof.
We only need to show the matrix A = ∂ J t · f∂x is symmetric. Since f =( f , ..., f k ), it is straightforward to verify that(6.7) A ij = ∂ ( (cid:80) k(cid:96) =1 f (cid:96) ∂f (cid:96) ∂x i ) ∂x j = k (cid:88) (cid:96) =1 ( ∂f (cid:96) ∂x i · ∂f (cid:96) ∂x j + f (cid:96) ∂ f (cid:96) ∂x i ∂x j ) = A ji . The eigenvalues of a real symmetric matrix are real and their eigenvectors areorthogonal. (cid:3)
Example f := ( x − x ) in Section 3. With β = 10 , atthe witness point ( − . , − . β ∂ J t · f∂x + I = (cid:32) . − . − . . (cid:33) , which has two eigenvalues 233 . . . , − . . , . . , . We are now ready to present thecompanion curve tracing method in detail. The input is a finite set f of polynomialsin R [ x ]. The algorithm starts by applying the criterion in Section 4 to determineif V R ( f ) is empty. If V R ( f ) is determined to empty, the algorithm terminates andreturn ∅ . Otherwise, it will produce given number of points in V R ( h ), where h = { f + w / √ β, . . . , f k + w k / √ β } for a large β (cid:29)
1, such that (cid:107) f (cid:107) ≤ (cid:15) holds. Moreprecisely, we have the following algorithm. • Algorithm
CompanionCurveTracing • Input: – a finite set of polynomials f = { f , . . . , f k } .4 W. WU AND C. CHEN – a prescribed number N . • Output: return ∅ if V R ( f ) is determined to be empty by the algorithm; oth-erwise return at least N points for each connected component of V R ( f ) suchthat for each point p , the backward error (cid:107) f ( p ) (cid:107) ≤ (cid:15) . • Steps:1. choose a random point a in R n and a large number β (cid:29)
12. apply Proposition 4.1 to determine if V R ( f ) is empty; if true, then return ∅
3. obtain the real solution set S of the square system ( x − a ) + β J t · f = 0by homotopy continuation method4. let S (cid:48) = { x : (cid:107) f ( x ) (cid:107) < (cid:15) } for some tolerance (cid:15)
5. if S (cid:48) = ∅ , then return ∅
6. for each x ∈ S (cid:48) x (cid:48) , a (cid:48) ) where ( a (cid:48) ,
0) is presum-ably in the tubular neighborhood of V R ( f + w/ √ β )6.2 find the unit eigenvector ∆ x of β ∂ J t · f∂x + I at x (cid:48) with eigenvalue c close to 16.3 update a (cid:48) = a (cid:48) + hc ∆ x and x (cid:48) = x (cid:48) + h ∆ x , where h is the step size6.4 refine x (cid:48) by Newton iteration with fixed a (cid:48) by solving ( x (cid:48) − a (cid:48) ) + β J t · f = 0.6.5 if the smallest eigenvalue of β ∂ J t · f∂x + I at x (cid:48) is close to zero, itindicates that ( a (cid:48) ,
0) may be out of the tubular neighborhood, gotostep 6.16.6 goto step 6.2 with updated ( x (cid:48) , a (cid:48) ) until we have N points of thecurve starting from x
7. goto step 6 until we enumerate all points in S (cid:48)
7. Examples.
In this section, we illustrate the effectiveness of the companioncurve method on tracing some curves defined by rank-deficient systems. For all theexamples below, we simply choose β = 10 . Example f be the discriminant of the characteristic polynomialof the following matrix x x
11 1 x . It was proved that f is a sum of squares [15], which is necessarily rank-deficient at anypoint of V R ( f ). The companion curve method generates a straight line as illustratedby Fig. 4. Example f := x x + 100 x − x x − x x x + 36 x + 9 x + 3 x x − x x − x x x + 24 x x + 48 x x + 71 x − x x − x + 2 x x + 12 x − x − x − x + 3 , x − x x x + 3 x x + 36 x − x x + 100 x + 48 x x − x x + 24 x x − x x x − x x − x + 12 x + 2 x x − x − x x + 71 x + 7 x − x + 1 , which defines a real algebraic curve according to [4]. However, on the other hand,this system has two polynomials with four variables x , x , x , x . Hence, the nullityof J f is at least 2, which implies that f must be rank-deficient. The companioncurve method generates several curves, as illustrated by Fig. 5. OMPANION CURVE Fig. 4 . (Color online) Left: The traced curve (in red • ) of V R ( f ) in Example 7.1 and itscompanion curve (in blue (cid:3) ) overlap. Right: zooming in on the left figure to see the difference. Fig. 5 . (Color online) Left: Projections in ( x , x , x ) of the traced curve (in red • ) of V R ( f ) in Example 7.2 and its companion curve (in blue (cid:3) ) overlap. Right: zooming in part of the leftbottom branch of V R ( f ) (and its companion curve) to see the difference. Example f : x − x x + 100 x x + 160 x x x + 64 x x + 64 x x − x + 20 x x + 112 x x + 16 x x − x x − x x + 50 x − x + 1 ,x x − x x x + 100 x x + 160 x x x + 64 x x x + 64 x x x − x x + 20 x x + 112 x x x + 16 x x x − x x x − x x x − x + 50 x x + 16 x x − x − x x − x x − x − x + 14 x − x − x − x − , − x x + 2 x − x + 1 , − x x − x x + 3 x + 1 , which also defines a real algebraic curve according to [4]. This system has fourpolynomials with six variables x , . . . , x . Hence, the nullity of J f is at least 2,6 W. WU AND C. CHEN
Fig. 6 . (Color online) Left: Projections in ( x , x , x ) of the traced curve (in red • ) of V R ( f ) in Example 7.3 and its companion curve (in blue (cid:3) ) overlap. Right: zooming in on part of the rightbottom branch of V R ( f ) (and its companion curve) to see the difference. which implies that f must be rank-deficient. The companion curve method generatesseveral curves, as illustrated by Fig. 6. Example f = x y + x z + y z − xyz + 1has 4 real isolated solutions: { (1 , , , (1 , − , − , ( − , , − , ( − , − , } and isrank-deficient at all these points. The real zero set of f can be seen as a degen-erate curve. Fig. 7 illustrates the points and their corresponding companion curveproduced by the method. Note that the companion curve bounces back and fortharound the isolated solutions. This is due to the fact that in Algorithm Companion-CurveTracing
Step 6 . a away from an isolated point whileStep 6 . a back to a neighborhood of the isolated point.
8. Conclusion and future work.
In this paper, we proposed a companioncurve method for tracing real algebraic curves whose defining system has singularJacobian along the whole curve. The effectiveness of this method is illustrated byseveral non-trivial systems. We notice that, even the algebraic form of the real al-gebraic curve is singular, as long as the curve is geometrically regular or has onlysingletons, the method works very well. We do notice that the method may not beable to obtain a complete tracing of curves which are not only “algebraically singular”but also “geometrically singular”, simple as ( xy ) , without an algebraic preprocess.This is due to the fact that a slight perturbation of the curve may generate morethan one connected components although the curve itself is connected via the singu-lar point. A preliminary investigation suggests that the symmetry matrix ∂ J t · f∂x + I plays a vital role in suggesting tracing directions, which deserves a further in-depthinvestigation. Acknowledgements.
This work was funded by NSFC (No. 11771421), the KeyResearch Program of Frontier Sciences of CAS (No. QYZDB-SSW-SYS026), CAS“Light of West China” Program, National Key Research and Development Program
OMPANION CURVE Fig. 7 . (Color online) Left: The four points traced (in red • ) of V R ( f ) and their companioncurve (in blue + ). Right: zooming in on the left figure around point (1 , , . (No. 2020YFA0712300), and Chongqing Programs (No. cstc2018jcyj-yszxX0002, No.cstc2019yszx-jcyjX0003, No. cstc2020yszx-jcyjX0005). REFERENCES[1]
E. Allgower and K. Georg , Introduction to Numerical Continuation Methods , vol. 45, 012003.[2]
C. Beltr´an and A. Leykin , Robust certified numerical homotopy tracking , Foundations ofComputational Mathematics, 13 (2013), pp. 253–295.[3]
L. Blum, F. Cucker, M. Shub, and S. Smale , Complexity and Real Computation , Springer-Verlag New York, Inc., Secaucus, NJ, USA, 1998.[4]
C. Chen, W. Wu, and Y. Feng , Full rank representation of real algebraic sets and applications ,in CASC 2017, Springer, 2017, pp. 51–65.[5]
C. Chen, W. Wu, and Y. Feng , Visualizing planar and space implicit real algebraic curveswith singularities , J. Syst. Sci. Complex., 33 (2020), pp. 1252–1274.[6]
J. Cheng, S. Lazard, L. Pe˜naranda, M. Pouget, F. Rouillier, and E. Tsigaridas , Onthe topology of real algebraic plane curves , Mathematics in Computer Science, 4 (2010),pp. 113–137.[7]
M.-D. Choi and T. Lam , Extremal positive semidefinite forms , Mathematische Annalen, 231(1977), pp. 1–18.[8]
D. Daouda, B. Mourrain, and O. Ruatta , On the computation of the topology of a non-reduced implicit space curve , in ISSAC 2008, 2008, pp. 47–54.[9]
H. Engl, M. Hanke, and A. Neubauer , Regularization of Inverse Problem , vol. 375, 01 1996.[10]
J. Eriksson, P. Wedin, M. Gulliksson, and I. Soderkvist , Regularization methods foruniformly rank-deficient nonlinear least-squares problems , Journal of Optimization Theoryand Applications, 127 (2005), pp. 1–26.[11]
A. J. Gomes , A continuation algorithm for planar implicit curves with singularities , Computers& Graphics, 38 (2014), pp. 365 – 373.[12]
J. D. Hauenstein , Numerically computing real points on algebraic sets , Acta ApplicandaeMathematicae, 125 (2012), pp. 105–119.[13]
H. Hong , An efficient method for analyzing the topology of plane real algebraic curves , Math-ematics and Computers in Simulation, 42 (1996), pp. 571 – 582.[14]
K. Jin and J. Cheng , Isotopic meshing of a real algebraic space curve , J. Syst. Sci. Complex.,33 (2020), pp. 1275–1296.[15]
P. Lax , On the discriminant of real symmetric matrices , Communications on Pure and AppliedMathematics, 51 (2005), pp. 577–586.[16]
D. Lazard and F. Rouillier , Solving parametric polynomial systems , J. Symb. Comput., 42 W. WU AND C. CHEN(2007), pp. 636–667.[17]
J. M. Lee , Introduction to Smooth Manifolds , vol. 218 of Graduate Texts in Mathematics,Springer-Verlag, 2013.[18]
T. L. Lee, T. Y. Li, and C. H. Tsai , Hom4ps-2.0: a software package for solving polynomialsystems by the polyhedral homotopy continuation method , Computing, 83 (2008), p. 109.[19]
B. Martin, A. Goldsztejn, L. Granvilliers, and C. Jermann , Certified parallelotope con-tinuation for one-manifolds , SIAM Journal on Numerical Analysis, 51 (2013), pp. 3373–3401.[20]
B. Mishra , Algorithmic Algebra , Springer-Verlag, New York, 1993.[21]
D. Mumford , Algebraic Geometry I: Complex Projective Varieties , Springer-Verlag, BerlinHeidelberg, 1995.[22]
F. Rouillier, M.-F. Roy, and M. Safey El Din , Finding at least one point in each connectedcomponent of a real algebraic set defined by a single equation , Journal of Complexity, 16(2000), pp. 716 – 750.[23]
A. Tikhonov, A. Goncharskij, V. Stepanov, and A. Yagola , Numerical Methods for theSolution of Ill-Posed Problems , 01 1995.[24]
W. Wu, C. Chen, and G. Reid , Penalty function based critical point approach to compute realwitness solution points of polynomial systems , in CASC 2017, Springer, 2017, pp. 377–391.[25]
W. Wu and G. Reid , Finding points on real solution components and applications to differ-ential polynomial systems , in ISSAC 2013, 2013, pp. 339–346.[26]
W. Wu, G. Reid, and Y. Feng , Computing real witness points of positive dimensional poly-nomial systems , Theoretical Computer Science, 681 (2017), pp. 217 – 231.[27]
Y. Yu, B. Yu, and B. Dong , Robust continuation methods for tracing solution curves ofparameterized systems , Numerical Algorithms, 65 (2014), pp. 825–841.[28]