134 661 THESlS Illiflllljfllllflifljflllfllljlll LIBRARY Michigan State University This is to certify that the dissertation entitled rSIOLVING POLYNOMIAL SYSTEMS IN C BY POLYHEDRAL HOMOTOPIES presented by King L1 has been accepted towards fulfillment of the requirements for P1“ - n degree in Mathematics— ! l 7 Major professor August 1,2000 Date MS U is an Affirmative Action/Equal Opportunity Institution 0- 12771 PLACE IN RETURN BOX to remove this checkout from your record. To AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE woo mm.“ SOLVING POLYNOMIAL SYSTEMS IN C” BY POLYHEDRAL HOMOTOPIES By Xing Li A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Mathematics 2000 ABSTRACT SOLVING POLYNOMIAL SYSTEMS IN C" BY POLYHEDRAL HOMOTOPIES By Xing Li In the last two decades, the homotopy continuation method has been developed into a reliable and efficient numerical algorithm for solving all isolated zeros of poly- nomial systems. During the last few years, a major computational breakthrough has emerged in the area. Based on the Bernshtein theory on root count, the polyhedral h0- motopy is established to considerably reduce the number of homotopy paths that need to be traced to find all the isolated roots, making the method much more powerful. The main goal of this dissertation is to present a strategy which uses homotopy continuation method efficiently to solve polynomial systems via mixed cell calculation. ACKNOWLEDGMENTS I would like to express my sincere gratitude to Professor Tien-Yien Li, my disser- tation advisor, for his constant encouragement and kind guidance. His knowledge, insights and enthusiasm were invaluable. I would also like to thank Professor William C. Brown, Professor Chichia Chiu, Professor Ronald Fintushel, Professor Michael Frazier, Professor Richard Hill, Pro- fessor Wei-Eihn Kuan and Professor Jay Kurtz for their time and advice. iii TABLE OF CONTENTS INTRODUCTION 1 Polyhedral Homotopy 1.1 Bernshtein Theory .................. ~ ............ 1.2 Polyhedral Homotopy ............................ 1.3 Solve Binomial System ............................ 2 Linear Programming 3 Find Mixed Cells 3.1 Find all lower Edges of a lifted Lattice Set ................. 3.2 Extend level-k Subfaces ........................... 3.3 Find All Fine Mixed Cells .......................... 4 Numerical Implementation 4.1 Algorithm .................................. 4.2 Implementation ................................ 4.3 Numerical Results .............................. BIBLIOGRAPHY iv Introduction Polynomial systems arise quite commonly in many fields of science and engineering, such as formula construction, geometric intersection, inverse kinematics, power flow with PQ-specified bases, computation of equilibrium states, etc., see [10]. Elimination theory-based methods, most notably the Buchberger algorithm [5] for constructing Griibner bases, are the classical approach to solving multivariate polynomial systems, but their reliance on symbolic manipulation makes those methods somewhat unsuit- able for all but small problems. In 1977, Garcia and Zangwill [14] and Drexler [11] independently presented the- orems suggesting that homotopy continuation could be used to find the full set of isolated zeros of a polynomial system numerically. During the last two decades this method has been developed into a reliable and efficient numerical algorithm for ap- proximating all isolated zeros of polynomial systems. See [23] for a survey. Let P(x) = 0 be a system of n polynomial equations in n unknowns. Denoting P = (p1, ..., p"), we want to find all isolated solutions of p1($1,...,$n) = 0 (1) 1),,(231, ...,:cn) = 0, for x = ($1, ..., :12"). The classical homotopy continuation method for solving (1) is to define a system that is easy to solve Q(x) = (q1(x), ..., qn(x)) = 0 and then follow the curves in the real variable t which make up the solution set of 0 = H(x, t) = (1 — t)Q(x) + tP(x). (2) More precisely, if the system Q(x) , known as the start system, is chosen correctly, the following three properties hold: a Property 1 (Triviality). The solutions of Q(x) = 0 are known. 0 Property 2 (Smoothness). The solution set of H (x, t) = 0 for 0 S t S 1 consists of a finite number of smooth paths, each parametrized by t in [0, 1). 0 Property 3 (Accessibility). Every isolated solution of H (x, 1) :2 P(x) = 0 can be reached by some path originating at t = 0. It follows that this path starts at a solution of H (x, 0) = Q(x) = 0. When the three properties hold, the solution paths can be followed from the initial points (known because of Property 1) at t = 0 to all solutions of the original problem P(x) = 0 at t = 1 using standard numerical techniques [1, 2]. A homotopy H (x, t) = 0 with H(x,0) = Q(x) and H(x, l) = P(x), which may not be in the form of (2), is considered to be successful if it satisfies these three properties. A typical choice [8, 22, 24, 30, 46, 47] of the system Q(x) = (q1(x), ..., q,,(x)) which satisfies Properties 1-3 is, q1(:1:1,...,:c,,) = alx‘f‘ — b1 (3) q,,(:r1, ..., as") = anzf‘," — I)", where d1, ..., d,, are the degrees of p1(x), ..., pn(x) respectively, and a,, b,- are random complex numbers (and therefore nonzero with probability one). So in one sense, the original problem posed is solved. All solutions of P(x) = 0 are found at the end of d1 - - - (1,, paths that make up the solution set of H(x, t) = 0,0 5 t S 1. 2 Solutions to Solutions to Stan System P(x) = O > infinity Q(X) =0 \/\/// t=0 t=l l Figure 1: Solution curves of H (2:, t) = 0 The reason the problem is not satisfactorily solved by the above considerations is the existence of extraneous paths. Although the above method produces d = d1 - - - d" paths since Q(x) = 0 in (3) has d isolated nonsingular solutions, the system P(x) = 0 may have fewer than d solutions. We call such a system deficient. In this case, some of the paths produced by the above method will be extraneous paths. More precisely, even though Properties 1-3 imply that each solution of P(x) = 0 will lie at the end of a solution path, it is also consistent with these properties that some of the paths may diverge to infinity as the parameter t approaches 1 (the smoothness property rules this out for t —> to < 1). In other words, it is quite possible for Q(x) = 0 to have more solutions than P(x) = 0. In this case, some of the paths leading from roots of Q(x) = 0 are extraneous, and diverge to infinity when t —) 1 (See Figure 1). Empirically, we find that most systems arising in applications are deficient. A great majority of the systems have fewer than, and in some cases only a small fraction of, the expected number of solutions. For a typical example of this sort, let us look at the following Cassou-Nogues system p1 = 15b4cd'z + 6b4c3 + 21b4c2d — 144b’c — 8b2c2e —28b2cde — 648b2d + 36b2d28 + 9b40l3 -— 120, p2 = 30b4c3d -— 32cde2 — 720b2cd — 24b2c3e — 432b2c2 + 576ce — 576de +16b2cd2e + 16d2e2 + l6c2e2 + We + 39b4c2cz2 + 18b4cd3 —432b2d2 + 24b2d3e — 16b2c2de — 240c + 5184, (4) p3 = 216b2cd — 162b2d2 — 81b2c2 + 1008ce — 1008de + 15b2c2de —15b2c3e — 80cde2 + 40dze2 + 40c2e2 + 5184, p4 = 4b2cd — 3b2d'z — 4b2c2 + 22ce — 22de + 261. Since d1 = 7, d2 = 8,d3 = 6 and d4 = 4 for this system, the system Q(x) in (3) will produce d1 x d2 x d3 x d; = 7 x 8 x 6 x 4 = 1344 paths for the homotopy in (2). However, the system (4) has only 16 isolated zeros. Consequently, a major fraction of the paths are extraneous. Sending out 1344 paths in search of 16 solutions is a highly wasteful computation. The choice of Q(x) in (3) to solve the system P(x) = 0 requires an amount of computational effort proportional to dl - - - d”, known as the Bézout number, which bounds the number of isolated zeros, counting multiplicities, of P(x) in C" [39]. We wish to derive methods for solving deficient systems for which the computational effort is instead proportional to the actual number of solutions. In the last few years, a major computational breakthrough has emerged in the area. The new idea takes a great advantage of the Bernshtein theory [4] which provides a much tighter bound, compared to the Bézout bound, for the number of isolated zeros of P(x) in the algebraic tori (C‘ )", where C“ = C \ {0}. The so called polyhedral homotopy [18] is then established for the new method and the homotopy paths so produced is much fewer. Accordingly, the required computation effort is considerably reduced. The new algorithm is very promising. In particular, for polynomial systems without special structures, the new algorithm outperformed the existing methods by a big margin. The purpose of this dissertation is to present a strategy of solving polynomial sys- tems by polyhedral homotopy efficiently via newly developed mixed cell calculation. The polyhdreal homotopy and some necessary terminologies are introduced in Chap- ter 1. In Chapter 2, we give a basic linear programming algorithm which serves as a main tool for the mixed cell calculation presented in Chapter 3. Our algorithms have been implemented successfullly, the numerical results on substantial variety of examples are presented in Chapter 4. CHAPTER 1 Polyhedral Homotopy The Bernshtein theory on root count of polynomial systems is essential for our attempt to reduce the number of homotopy curves need to be traced when the homotopy continuation method is employed to find all isolated zeros of polynomial systems. In the first section of this chapter, the Bernshtein theory on root count in (0')”, where C“ = C\ {O} , as well as its extension to root count in C" are presented. In the second section, the polyhedral homotopy, based on the Bershtein theory, for finding all isolated zeros of a polynomial system is introduced. In the last section, we will disscuss how to solve a binomial system to obtain initial solutions of a polyhedral homotopy. 1.1 Bernshtein Theory Let the given polynomial system be P(x) = (p1(x), - -- ,pn(x)) E C[x], where x = (2:1,n- ,3“). With at" = 22‘," ”4:3," where a = (a1,--- ,an), write p1(X) : 26111238: aESI (1.1) pn(x) = 2 Clay, :6 Sn where 51,-” ,5" are fixed subsets of N" with cardinals k,- = #Sj, and c}... E C“ for a E S,-,j = 1,--- ,n. We call S,- the support of pJ-(x), denoted by supp(p,-), its convex hull K,- = conv(.S'j) in IR" the Newton polytope of p,-, and S = (5'1, - - - ,3”) the support of P(x), denoted by supp(P). We now embed the system (1.1) in the system P(c,x) = (p1(c,x), - -- ,pn(c,x)) where p1(c,X) = Z cum“. aESI (1.2) pn(C.X) = Z amar‘, lESn and the coemcients 01',- with a E Sj, for j = 1,-o- ,n in the system are taken to be a set of M E k1 + - - - + kn variables. Namely, the system P(x) in (1.1) is considered as a system in (1.2) corresponding to a set of specified values of coefficients 6 = (egg) or P(x) = P(c,x). We shall refer to the total number of isolated zeros, counting multiplicities, of a polynomial system as the root count of the system. Lemma 1.1 [17] For polynomial systems P(c, x) in (1.2), there exists a polynomial system G(c) = (gl(c),--- ,g,,(c)) in the variables c = (9,.) for a E S,- and j = 1, - - . ,n such that for those coefi‘icients E: = (c;,.) for which G (E) 75 0, the root count in (C‘ )u of the corresponding polynomial systems in (1.2) is a fixed number. And the root count in ((C" )" of any other polynomial systems in (1.2) is bounded above by this number. Remark 1.2 Since the zeros of the polynomial system G (c) in the above lemma form an algebraic set with dimension smaller than M, its complement is open and dense with full measure in CM . Therefore, with probability one, G(é) 75 0 for randomly chosen coefficients 6 = (ch) 6 (CM . Hence, polynomial systems P(é,x) in (1.2) with G(c) # 0 are said to be in general position. Theorem 1.3 ([4], Theorem A) The root count in (0')" of a polynomial system P(x) = (p1(x),...,p,,(x)) in general position equals to the mixed volume of its support. The terminology in this theorem needs explanation. For non-negative variables A1, - - . , A" and the Newton polytopes KJ- ofpj, for j = l, - -- ,n, let A1K1+- - -+/\,,K,, denote the Minkowski sum of A1K1,- -- , AnKn, that is, AlKl +"°+AnKn = {A1r1+---+Anr,,|rj E Kj,j = I,’” ,n}. It can be shown that the n-dimensional volume of this polytope Voln(A1K1 +- - ° + [\nKn) is a homogeneous polynomial of degree n in A1, - - - , A”. The coefficient of the term Al x - - - x A“ in this homogeneous polynomial is called the mixed volume of the polytopes K1, - - - , Kn, denoted by M(K1, ' - - ,Kn), or the mixed volume of the support of the system P(x) = (p1 (x), - -- , pn(x)), denoted by M(Sl, - - - ,3“) where 53- = supp(p,—) for j = 1, - -- ,n. Sometimes, when no ambiguities exist, it is called the mixed volume of P(x). In [6], this root count was nicknamed the BKK bound after its inventors, Bern- shtein [4], Kushnirenko [21] and Khovanskii [20]. In general, it provides a much tighter bound compared to variant Bézout bounds [32, 39]. An apparent limitation of the theorem is that it only counts the isolated zeros of polynomial systems in (C‘ )" rather than all the isolated zeros in the amne space C”. For the purpose of finding all the isolated zeros of a polynomial system in C", a generalized version of the theorem which counts the roots in C" is strongly desirable. This problem was first attempted in [36] where the notion of the shadowed sets was introduced and a bound for the root count in C" was obtained. Later, a significantly much tighter bound was discovered in the following theorem. Theorem 1.4 [27] The root count in C" of a polynomial system P(x) = (p1(x),--- ,pn(x)) with supports 5', = supp(p,~),j = 1,~-- ,n is bounded above by the mixed volume M(Sl U{0}, . . . , Sn U{0}). In other words, the theorem says that the root count in C“ of a polynomial system P(x) = (p1 (x), - - - , pn(x)) is bounded above by the root count in (C‘ )” of the poly- nomial system P(x) in general position obtained by augmenting the constant term to those p33 in P(x) in which the constant term is absent. As a corollary, when 0 E S,- for all j = 1, - - - ,n, namely, all p,- (x) in P(x) have constant terms, then the mixed volume of P(x) also serves as a bound for the root count of P(x) in C", rather than in (C‘)” as Theorem 1.3 asserts. This theorem was further extended in several different ways [19, 37]. 1.2 Polyhedral Homotopy In light of Theorem 1.4 given in the last section, to find all isolated zeros of a given polynomial system P(x) = (p1(x), - - - ,p,,(x)) in C" with support S = ($1, - - - , Sn), we first augment the monomial x° (=1) to those 12,- ’s which do not have constant terms. Followed by choosing coefficients of all the monomials in the system generically, a new system Q(x) = (q1(x), - - . ,q,,(x)) with support 5" = (Si, - - - ,S,’,) is obtained, where, of course, 8'; = Sj U {0} for j = l, - -- ,n. We will solve this system in the first place, and the details will be discussed in this section. Afterwards, in Chapter 4, we will present our algorithm to solve P(x) = 0. To begin, we write (11(X) = 251,96. 36$; Q(X) = 5 (1-3) Qn(x) = Z Emlxa' aES; Since all those coeficients 6,3,, for a E S}, j = 1, - - ~ ,n, are chosen generically, this system may be considered as a system in general position. Namely, there exists a polynomial system 0(0) = (91(C)w-- ,gm(<=)) (1-4) in the variables c = (cw), for a E S}, j = l,--- ,n, such that polynomial sys- tems with G(c) 74 0 reach the maximum root count in (C‘ )" for the support 5’ = (Si,"' ,5;) and we have 0(5) 79 0 for the set of coefficients 6 = (5,3,) in (1.3). Let t denote a new complex variable and consider the polynomial system Q(x, t) = (61(x, t), - - - ,cjn(x,t)) in the n + 1 variables (x, t) given by @109 t) _____ Z Elfixatwfla), .653 520“) = s (1.5) (Mint) = Z 5n..X‘t“"‘“’, 365% where each w,- : S; —-> R for j = 1, - .. ,n is chosen generically and known as a lifting on S}. For a fixed to, we rewrite the system in (1.5) as Z (51,.t301(8))xa : the S] 61(xa t0) Q(X,to) = qn(x, t0) = Z (5",.t3’n(l))xa. 1165;, 10 This system is in general position if for G (c) in (1.4), T(to) ..=. G(a,,.t;,""(") 7e 0, for a e 5;. and j = 1, . .- ,n. The system T(t) = 0 can have at most finitely many solutions, since T(t) is not identically 0 because T( 1) = G(Ej,.) 75 0. Let t1 = rlewli ' ' ' )tk : rkewh be the solutions of T(t) = 0. Then, for any 0 75 0,- for j = 1, - -- ,k, the systems Q(x,t) = (q1(x3 t): ' ' ' tq-fl(xit)) given by q'1(x,t) ___ 2(El’aeiw1(a)9)xatw1(a), 36$] Q(x, t) = q" (x, t) = Z (en’aeiuh; (3)0)th‘Wn (a) , :6 5;, are in general position for all t > 0 because Ej neiwj (a)0tw,- (a) : Ej a(tei9)tvj (a) and, G(Ej,.(te‘9)"’1'(‘)) = T(te‘a) 7s 0. Therefore, without loss of generality, (choose an angle 6 at random and change the coeficients 6,3. to éj,.e“"1'(“)9 if necessary) we may suppose the systems Q(x, t) in (1.5) are in general position for all t > 0. Together with Lemma 1.1 given in the last section, it follows that for all t > 0 the systems Q(x, t) in (1.5) have the same number of isolated zeros in (C‘ )“. This number, say It, should equal to the mixed volume of the support of Q(x) in (1.3) by Theorem 1.3. We shall skip this fact temporarily and will reach this assertion at the end of this section. Now, consider Q(x, t) = 0 as a homotopy, known as the polyhedral homotopy, defined on (C‘ )" x [0, 1]. We have Q(x, 1) = Q(x), and the zero set of this homotopy 11 is made up of k homotopy paths, say, x1(t),--- ,x’°(t), since for each 0 < t S 1, Q(x, t) has exactly I: isolated zeros from the argument given above. Since each q“,-(x, t) has nonzero constant term for all j = 1, - - - ,n, by a standard application of generalized Sard’s Theorem [7], all those homotopy paths are smooth with no bifurcations. Therefore, both Property 2 (Smoothness) and Property 3 (Accessibility) introduced earlier hold for this homotopy. However, at t = 0, Q(x, 0) E 0, see Figure 1.1. Consequently, the starting points x1(0), . - - ,x"(0) of those homotopy paths can not be identified, causing the breakdown of the standard homotopy continuation algorithm. This major obstacle can be overcome by the devise we describe below. A Q(x.0) =0 Q(X) =0 Figure 1.1: Solution curves of Q(x, t) = 0 For a = (0:1,- -- ,an) E R", consider the transformation y = t‘ax defined by 3]] = t—Ollxl, (1.6) yn = t'“"z,,. 12 For a: (a1,--- ,an) E N“, we have . — aloe. a" : (yltai )ai , , , (ynt°")“" — “1 . . . an aian+~~+a a yl n t n n = yat(a,a) . Here, (-, ) stands for the usual inner product in R". Substituting (1.7) into (1.5) yields, for j: 1,--- ,n My, t) E ij(yt°‘, t) = Z: Ej,ay‘t(°‘")t‘”"°’ :65;- = Z Ej.yat<(a.1).(a.w,-(a))> 365;. (1.8) = Z Ej,ay.t—aj :65;- = Z 52-..1"+ Z 5:..y‘t<“"""’l (1.12) ass; ass; (a,a)=s,- (61.5))[31' Evidently, for any path y(t) defined on [0, 1], we have, for all t > 0, H.(y(t),t) = 0 «=> Heart) = 0. Therefore, the zero set of Ha(y, t) = 0 consists of the same homotopy paths of the homotopy H (y, t) = 0 in (1.9). The difference is, the starting points of the homotopy paths considered in the homotopy Ha, (y, t) = 0 are solutions of the system f h?(}’, 0) = 2 51,83" = O: 1165] (ata)=Bl Ha(y,0) = < i (1.13) h:(yt 0) = Z grimy. : 0- aESQ l =a. As shown below, when this system is in certain desired form, its isolated nonsingular solutions that lie in (C‘ )" can be constructively identified. In those situations, Prop- erty 1 (Triviality) becomes partially valid for those homotopy paths of Ha(y, t) = 0 that emanate from those nonsingular solutions of (1.13) in (C‘)", and we may follow those paths to reach a partial set of isolated zeros of Q(y) at t = 1. The system (1.13) is known as the binomial system if each h?(y,0) consists of exactly two terms, that is, hi’(y,0) = ciy‘“ + c’w"1 = 0, (1.14) h:(ya0) : Cnyan + any“; = 0a 14 where aha; E S}, c,- = 51.41,- and c; = 51's;- for j = 1,--- ,n. And in this case, ({a1,a’1},--- ,{ama:,}) is called a mixed cell (of type (1,--- ,1)) of S’ = (Si, - - - , 5;) associated with inner normal 07. Proposition 1.5 The binomial system in (1.14) has al —a’l det f (1.15) a. D “I nonsingular solutions in (0')". The number kc is called the volume of the mixed cell ({a1,a’1},- .. , {a,,,a;}). The proof of this proposition is constructive and therefore provides an algorithm for solving the binomial system (1.14) in (C‘ )”. We will come back to this matter in the next section. In summary, for given a = (a1, - - - ,an) 6 R", by changing variables y = t‘ax, as in (1.6), in the homotopy Q(x, t) = (61(x, t), - -- ,(jn(x,t)) = 0 in (1.5), the homotopy H(y, t) = (h1(y,t), - -- ,hn(y,t)) = 0 in (1.9) is obtained, where hJ-(y,t) = q‘j(yt°‘,t). Followed by factoring out the lowest power tfif of t among all monomials in each individual h,-(y,t) = 0 for j = l,--- ,n we arrive at the homotopy Ha(y, t) = 0 in (1.11). When the start system Ha(y,0) = 0 of this homotopy is binomial, its nonsingular solutions in (C‘ )“, kc (as given in (1.15)) of them, become available. We may then follow those homotopy paths of Ho, (y, t) = 0 originated from those ka regular solutions of Ha(y,0) = 0 in (0')", and reach kc. isolated zeros of Q(y) at t = 1. Worth notifying here is the fact that the system Q(x), or Q(y), stays invariant at t = 1 during the process. Now, the existence of oz 6 IR“ for which the start system H, (y, 0) = 0 is binomial is warranted by the following 15 Proposition 1.6 For all the real functions wj : S; ——> R, j = 1, - -- ,n being gener- ically chosen, there must exist a E R" , for which the start system Ha(y,0) = 0 of the homotopy Ha(y, t) = 0 in (1.12) is binomial with a nonempty set of nonsingular solutions in (C’)", i.e., ha 79 0 in (1.15). The assertion of this proposition was proved implicitly in [18] with terminologies and machineries developed in combinatorial geometry, such as, random liftings, fine mixed subdivisions, lower facets of convex polytopes, etc., see also [23]. Here, we elect to reinterpret the result without those specialized terms. V Now, different a 6 IR" given in Proposition 1.6 leads to different homotopy Ha(y,t) = 0 in (1.11). Henceforth, following homotopy paths of those difl'erent homotopies will reach different sets of isolated zeros of Q(y). By taking the Puiseux series expansions of those homotopy paths of Ho, (y, t) = 0 originated at (0‘)" into consideration, it is not hard to see that those different sets of isolated zeros of Q(y) reached by different sets of homotopy paths actually disjoint from each other. Most importantly, it can be shown that every isolated zero of Q(y) can be obtained this way by following certain homotopy curve of the homotopy Ha(y, t) = 0 associated with certain a 6 IR“ given by Proposition 1.6. Thus the total number of isolated zeros of Q(y) must equal to the sum of those ha ’8 corresponding to all the possible a’s provided by Proposition 1.6, respectively. In [18], it was shown that this sum actually equal to the mixed volume of Q(y). This yields an alternative proof of Theorem 1.3, it is very different from Bernshtein’s original approach [4]. 1.3 Solve Binomial System Another major step in solving polynomial systems by using the polyhedral homo- topy method as we described in the previous section is finding the solutions of the 16 corresponding binomial system cly‘ll + c’lyall = 0, (1.16) any“ + 61.31“ = 0. produced by the mixed cell ({al, a’l}, - -- ,{am a“) as in (1.14). We now discuss the method for solving (1.16) in (C‘ )“. Let _ __’ J, 3:1:H'1n1 and, with y E (0)" in mind, we rewrite the system (1.16) as yvl = b1: (1.17) yv" : bm where b,- = ——j for j = 1,--- ,n. Let 1' V = vf v; v: (1-18) and for brevity, write yV___(yv1,.H,yun) and b=(bli°°'ibn)- Then, (1.17) becomes, yV = b. (1.19) With this notation, it is easy to verify that for an n X n integer matrix U, we have, M)” = W”). 17 Now, when the matrix V in (1.18) is an upper triangular matrix, i.e., F - v11 012 ' ' ° Uln 0 022 ' ' ' ‘Uzn V = , L O . . . 0 ”an d then the equation in (1.19) becomes yin] : b1, ymzyvzz : 02, 1 2 (1.20) slimy?" - - - 3.13.“ = b... By forward substitutions, all the solutions of the system (1.20) in (C‘ )“ can be found, and the total number of solutions is |v11| x - - - x |v,,,,| = [det V|. In general, we may upper triangularize V in (1.18) by the following process. Recall that the greatest common divisor d of two nonzero integers a and b, denoted by gcd(a, b), can be written as d = gcd(a, b) = ka + lb, for certain nonzero integers k and I. Let k l M = __9 9. d d We have det(M) = 1, and a k l a d M = :: b —g g b 0 Similar to using Givens rotation to produce zeros in a matrix for its QR factorization, the matrix M may be used to upper triangularize V as follows. For v E Z“, let a 18 and b be its i-th and the j-th (nonzero) components where i < j, that is, a -—> i-th v = b -—> j-th. With d = gcd(a, b), we let i-th j—th r 1 - l k l i-th 1 U(i,j) = (1.21) 1 -3 3 j-th 1 - 1 . 19 Evidently, U (i, j) is an integer matrix with |det(U (i, j))| = 1 and F' . T d i-th U (231)?) o j-th h d Thus a series of matrices in the form of U (i, j) in (1.21) may be used to successively produce zeros in the lower triangular part of the matrix V in (1.18), resulting in an upper triangular matrix. In simple terms, we may construct an integer matrix U , as a product of those U (i, j) ’s, with |det U | = 1 and UV is an upper triangular integer matrix. Now, as mentioned above, the solutions of the system (zU)V .—_ zUV = b (1.22) in (C‘ )" can be found by forward substitutions, since UV is an upper triangular integer matrix. And the total number of solutions in (0')" is |det(UV)| = |det(U)| . |det(V)[ = |det(V)[. By letting y = zU for each solution z of (1.22) in (C' )", we obtain all the solutions of the system (1.22) in (C‘ )", and hence, solve the system (1.16) in (C‘ )”. 20 CHAPTER 2 Linear Programming As outlined in the last chapter, when the polyhedral homotopy is employed to find all the isolated zeros of a polynomial system P(x) = (p1 (x), - - - , pn(x)) with supports 31,- -- ,5“, one major step is to indentify the mixed cells ({a1,a’1},--- ,{an,a:,}) induced by generic liftings w,- : Sj —> R for j = 1, - - - ,n. As a point of departure in deve10ping our algorithms for finding all mixed cells in the next chapter, we introduce in this chapter some basic terminologies and algorithms in linear programming [3] that will be used in the method. Consider the model problem min (f, x) (2.1) s.t. (c,,x) S b,-, i = 1,--- ,m where f 6 IR", c, 6 1R“,b = (b1,--- ,bm)T 6 R“, x = (x1,--- ,xn), m > n. The feasible region of (2.1), denoted by R, defines a polyhedral set. By a nondegenerate extreme point of R we mean a vertex point of R with exactly n active constraints. Let x0 be a nondegenerate extreme point of R and J = {j1, - -- , jn} be the set of indices of currently active constraints at x°, that is, (c,,x°) = b,-, ifi E J (c,,x°) < 0;, ifz ¢ J. 21 Let DT = [c,-,,- -- .cjn]. Since x0 is a vertex point, D must be nonsingular. Let D"1 =[u1,-~- ,un]. Then for any a > 0 and 1 S k S n, we have (cj,,x° - em.) = (c,-,,x°) — a(c,-,,u),) = (cinxo) = b,-, ifi at k, (2 2) (Cijo — em.) = (cj,,x°) — o(c,,,u,,) = bit. - 0’ < bit. and for small a > 0 , (c,-,x0 — em.) = (c,,x°) — o(c,-,u,.) < b,-, for i ¢ J. 0 Thus the n edges of the feasible region R emanating from x can be represented in theform xo—ouh, 0>0, k=1,--- ,n. These edges provide possible search directions to minimize the objective function (f, x). Let x1 = x0 — on.- with or > 0. Then the value of the cost function at x1 is (f,x1) = (f,x°) — a(f,u,-), and it decreases when (f, u.) > 0. It can be easily shown that x0 is an optimal solution of (2.1) if (f, u.) S 0 for all i = 1, - - - ,n. If some of the (f, 11,-) ’s are positive, then the greatest rate of decrease of the cost function is obtained by choosing k such that (f,u,,) = max{(f,u,-) | 1 S i g n}. Let s = u), be the next search direction. From (2.2), for all positive a, the i- th constraint is still active at x1 = x° - as for every i 6 J\{jk}, and the jk-th constraint becomes inactive but stays feasible. To make xl feasible, we must choose a > 0 such that (c,,x° -— as) = (c,,x°) -— o(c,~,s) g b,, for i 5! J. (2.3) 22 If (c,-,s) Z 0 for alli ¢ J, then the inequalities in (2.3) are valid for all a > 0 and problem (2.1) is unbounded from below with no solution. Otherwise, from (2.3), the largest possible a for x1 to stay feasible is (cit X0) '— bi (c,,s) all i 5! J with (c,,s) < 0}. oo=min{ Let I be the smallest integer such that (cl, X0) — bl 0° = ° — cos is a new extreme point of the feasible region R in (2.1) with Then x1 = x reduced value of the objective function. This procedure can be continued until either an optimal solution is reached or the problem is determined to be unbounded from below. We summarize the above discussion in Algorithm 1 below [3]. Algorithm 1 Solving the model problem (2.1) . Step 0: Initialization. Start with an extreme point x0 of (2.1), J = {i1,--- ,in}, And D‘1 = [u,,-] = [u1,--- ,ufl], where DT = [d1,--- ,d,,] = [c,-,,--- ,c,,,] is nonsingu- lar. Step 1: Computation of the search direction 3. Determine the smallest index k such that (f,u;,) = max{(f,u,-) |i= 1,--- ,n}. If (f, uk) 5 0, stop with optimal solution x0. Otherwise, set s = u), and go to Step 2. Step 2: Compute the maximum feasible step size a. If (c,,s) Z 0 for all i = l,--- ,m, print the message “problem is un- bounded from below” and stop. Otherwise, compute the smallest index I 23 Step 3: and a such that (m. X") - b1 (c..X°) - b.- a : (C115) = min { (Ci, 3) and go to Step 3. alli 6! J with (c,,s) < 0}. Update. Set x0 := x0 — as. Replace k-th column of DT by c; and update the inverse D’l . Replace the k-th element of J by l . Go to Step 1. The process of obtaining next feasible solution from a given feasible solution with one execution round of Step 1, 2 and 3 is called a pivot operation in Linear Program- ming. 24 CHAPTER 3 Find Mixed Cells In this chapter we will elaborate our algorithms for finding all mixed cells by solving a series of linear programming problems. For i = 1,- - - , n, let S,- be the support of p,-(x) in the polynomial system P(x) = (p1(x), - -- ,pn(x)) and w,- : S, ——> IR be a generically chosen function. Let A S,-={a=(a,w,-(a)) |8€S§}, f01”l=1,°°° ,n and for a = (ah-u ,an) E IR", write 6: = ((1,1). Recall that a mixed cell of S = (51,- .. ,S,,) induced by the lifting w = (w1,~- ,m,.) is a collection of pairs {a1,a'1},--- ,{an,a:,}, with aha] E S,, i = 1,--- ,n such that there exists an a = (011,-.- ,an) 6 IR“ for which (find) = (52:51), i: 1r” in and (as) > (“,a) for as s,\{a.-,a;}, 2': 1,... ,n. The geometric meaning of finding those mixed cells is that with generic lifting w, on lattice points S, C N“ for each i = 1, - . - ,n, we are looking for hyperplanes with 25 V Figure 3.1: A lifting on lattice points normal 0‘: = ((1,1) where a 6 IR", and each hyperplane supports the convex hull of S,- at exactly two points {5.35;} of S5, for each i = 1, - - - ,n, as shown in Figure 3.1. For 1 S i S n, £2 = {as} C S,- is called a lower edge of S, if there is a vector (3: = (a, 1) with a 6 R“ such that (as) = (axes (as) 5 (Ba), 6 e s.\{a,a'}. For 1 _<_ k _<_ n, E). =(é1,--- ,ék) where e,- = {are} C S}, for i = 1,-~- ,k, is called a level-k subface of S = (S1, - - - ,Sn) if there is a vector 61 = (a, 1) with 0: 6 R" such that for alli=l,--- ,lc, Obviously, a level- 1 subface of S is just a lower edge of S1 and a level-n subface of S induces a mixed cell of S. Thus, to find the mixed cells of S, one may proceed 26 by finding all the lower edges of S for i = l, - - - ,n in the first place, followed by extending the level-k subfaces of S from k = 1 to k = n. It can be shown that the mixed volume of S M(S) = (—1)"‘1 ivomm) + (—1)"‘2 ZvoMK. + K,) i=1 i —b. 2 min { -b. ((31, 8) (Ci: 8) Go to Step 4. alli ¢ J with (c,,s) < 0}. Set Xo == Xo — as and update J = {i1,--- ,iv} and D4. Set 5(8) := as) 11 (Pm {{a.,a,}|k,z e J}), and ’P := r\{{a,.,a,}|k,z 6 J} Go to Step 2. 3.2 Extend level-k Subfaces For a level-h subface E). = (é1,--- ,ék) of S =(S1,o~ ,Sn) with l g k < n where éi = {5:35;} 6 £(S'i) for i = 13 ° ' ' 1k) we say élc-i-l = {ék+ltélc+l} E £(§k+1) extends 33 E), if EH1 = (é1,~- ,ék+1) is a level-(k + 1) subface of S. Let 6(a) = {strain} E as“) {attain} extends a}. E). is called extendible if 8(1331.) 96 (b, it is nonextendible otherwise. To find all mixed cells of S = (Sh-u ,5“), we will start from k = 1 and extend E). step by step. If E). is nonextendible, there is no mixed cell of S which contains ({a1,a’1},- - - ,{ah afl) , and extension attempt will be repeated on the next Eh. Ob- viously, when k = n -— 1, an extendible E). yields mixed cells of S with elements in 8 (BE) (possibly several). In this section, we describe our algorithm to calculate £(Ek) efficiently for a given level-k subface E), = (é1,--- ,él.) where e,- = {54,53 C S,- for i = 1, - -- ,k. Now, consider the following system in the n + 1 unknowns 00,011, - - - ,an (a, a) 2 010, a e s... (and) S (é,d), 563i) fOIi=1,°'°,k (37) (éhé) : ($16!): 7:: 11'” ,k, where 61 = (011,- -- ,an, 1) E Rn“. The following lemma is obvious. Lemma 3.5 For a level-k subface E). = (é1,--- ,ék), where 6,- = {$4,512} C S,, if system (3. 7) has a solution (ao,a1, - .. ,0“) such that (as) = (6,307) = 0:0 for s,,s, 6 31.4.1, then {aha} extends 137).. With S,- = {a,-,1,--~ ,n,,mi} for i = 1,... ,n and a = (011,”- ,an) we may rewrite system (3.7) as |/\ (ah+1,j,0‘> — 00 —’wk+1(3k+1,j) j = 11' " amlc+1 (at—aijia) S wi(ai,j)—wi(ai) j: 1"” 1mi1 aid 6 Si\{aiia~:}a z=l,--- 3k (as-aha) = w.-(a$)-w.-(aa) i=1.--- .k- 34 By using the last lc equality constraints to eliminate k variables of a, the above system can be reduced to the following general inequality system : c’ldlaj, + c’uaaj, + + Cid-way", 5 b1 62.13091 + 62¢th + + C’2,j,,.‘1j.,: S 02 (3.8) c’,,,J-loz,-1 + cimaj, + + c’p’jncjn, S by k+1 wherep= Emi—Zk and n’=n—k+1. As beforbfby a coordinate transformation (x1, ~ -- ,x,,:) = (ah, - - - ,a,-fl,)L, where L 6 IR”, x", is nonsingular, the system can be further reduced to the following inequal- ity system: 01,1331 S b1 02,131 + 02,232 S 02 (3.9) ens-’01 + 672.2932 + + comma S 5» cmlxl + Cp,2$2 + + c,,,,,x,7 g by. Lemma 3.6 System (3.9) has a solution (ao,al,--- ,an) satisfying (infirm-,3) = (5h+1,j,&) = 0:0 for i S i, j g ml.“ if and only if system (3.9) has a solution x1, - - . ,x,, such that c,-,1x1+c.-,2x2+- - ~+c.,,,x,, = b,- and cj,1x1+c,-,2x2+- ' ~+c,-,,,x,, = b,-. Inequalities in (3.9) defines a polyhedron R in R” , and for an extreme point, or a vertex, of R, there are at least 17 active constraints. From Lemma 3.5 and Lemma 3.6, it follows that Lemma 3.7 If x0 is an extreme point of R and J = {i1, - .. ,it} is the indices of active constraints of the system at x0 with t Z n, then {a.+,,,p,a,.+,,-,} extends E), forany e Jn{1,---.m...}. 35 Similar to the discussion following Lemma 3.3, if {a.+1,,-,s,.+1,,.} E £(Sk+1) for {i,j} C I E {1, - - - ,m,.+1}, it will lead to a corresponding point x0 of R, its indices of active constraints includes {i, j}. Hence, to construct £(Ek) C [.(SkH), we may look for all those extreme points of R whose indices of active constraints contain at least a pair of {i, j} in I. To achieve this, we may certainly apply the Two-Point Test introduced in the last section and confine to I the indices of the “two points” to be tested. However, it is very likely that most of the 61,111,,- ’s that appeare in the pairs in £(SH1) fail to extend E). with their associated pairs in £(S).+1). Namely, those points do not exist in any of the pairs in 8 (Eh). This phenomenon never occurs when we compute the set of lower edges £(B) of B in the last section since all the points in B are extreme points. Consequently, every point of 3 appears in certain pairs of 5(3). From this important observation, we introduce the following One-Point Test to be used in additional to the Two-Point Test in our algorithm. One-Point Test Problem: min —C.-.,la:1 — 610.2932 — ---— 6.0.01», 01,131 S 01 02,1501 + 62,2502 S 02 (3.10) cmlxl + (3,7,2232 + ---+ cmxn S b,, cmlxl + cmxg + + cmx" S by where 1 3 i0 < my“. Lemma 3.8 Given 1 5 i0 < my“, if the optimal value of system (3.10) is greater than -b,-o, then {a.+1,,-,,a,.+l,,-} does not extend E), for all i E {1, - -- ,mk+1}\{io}. PROOF: Suppose there exists 1 S jo 3 ml,“ for which {ak+1,,o,a,.+1 .10} extends Eh. 36 By Lemmas 3.5 and 3.6, system (3.9) has a solution (x1, . -- ,xfl) satisfying 610,131 + 010,232 ' ° ' "r' 0&0,an = bio: 010,131 + 011.3932 +Cjom$n = bjo~ Hence the objective function value at (x1, - - - :30) is _ci01131 — 63-01232 — . o o — €10,713" = _biO’ which contradicts the fact that the optimal value of the system (3.10) is greater than —bio . D From the above lemma, points appeared in the pairs in £(S1.+1) may be tested systematically by using One-Point Test to check the possibilities of their appearances in the pairs in £(Ek) . When the optimal value obtained is not as desired for a particular point £119,130 , all the pairs associated with aptly-0 in £(Sk+1) should be deleted from further considerations. In the meantime, in the process of reaching the Optimal value of the problem, newly obtained extreme points of R provide a collection of new pairs of £(Ek) as long as their active constraints contain a pair of {i, j} in I = {1,--- ,m,.+1}. Furthermore, we no longer test points am,- in £(Sk+1) whose index i have appeared in any of the indices of the active constraints of the extreme points of R being obtained. The system of constraints in problem (3.10) is the same inequality system in (3.9) which defines the polyhedron R in IR". To find an initial extreme point of R to start Algorithm 1 on the problem, we may employ the same strategy by augmenting a new variable 6 Z 0 as in calculating 16(8) of B in the last section. Two-Point Tests will be used only after One-Point Tests have exhausted all the testings on possible candidates. Our experiences show that the Two-Point Test only plays a minor role in constructing 8(Ek) , namely, when we finish the One-Point Tests, most of the pairs in £(Ek) have been found. 37 Combining the One-Point Test and the Two-Point Test, we list the following al- gorithm for constructing 8 (173),) Algorithm 3 Given Eh, construct E (1:31,) Step 0: Initialization. Set up the inequality system (3.9). Start from an extreme point x0 with J = {i1,--~ ,in} and D“1 = [u1,--- ,u,,], where DT = [c,',,~- ,c,-,,], and set PHI := [.(SHI). Step 1: One-Point Test Problems. Step 1.0 Set i0 := 0, go to Step 1.1. Step 1.1 Set up objective function Find 1' = min {j I j > i0 and {ak+,,,-,a,.+l,j.} C F1,“ for some j’}. If such 7' does not exists, go to Step 2. Otherwise set i0 := 1' and f = (—c,-0,1,--- ,—c,-0,,,), go to Step 1.2. Step 1.2 Determine the smallest index k such that (1311),) : max{(f:ui> Ii: 13' ' '17)}- If (f,u),) S 0, go to Step 1.5. Otherwise, set s = u)c and go to Step 1.3. Step 1.3 Compute the smallest index 1 and a such that 0 _ . 0 ..., _(cbx) b1=min{(c,,x) b, _ (cbs) (c,,s) alli ¢ J with (c,,s) < 0}. Go to Step 1.4. Step 1.4 Set x0 := x0 — as and update J = {i1,--- ,in} and D‘l. If I < mkH, check if any lower edge {5k+1,1,5k+1 J} in 13),“ 38 extends 13).“. Collect these lower edges, if they exist, and delete them from EH1. Go to Step 1.2. Step 1.5 If the current value of objective function is not equal to -b,o, delete all lower edges containing point 6H1,“ from 13),“. Go to Step 1.1. Step 2: Two-point Test Problems. Step 2.1 Set up objective function. If R1.“ = 0, stop. Otherwise select a lower edge {ék+1.iovél¢+1.jo} 6 13hr Set f == (_C‘ioJ — €30.11“ 1-CioJI — 010.17): and 13”,“ :2 13“,,1\{a,.+1,,-0,a,.+1,,-,}, go to Step 2.2. Step 2.2 Determine the smallest index k such that (f,u;,) = max{(f,u,~) I i: 1,--- ,n}. If (f, uk) 3 0, go to Step 2.1. Otherwise, set s = u), and go to Step 2.3. Step 2.3 Compute the smallest index I and a such that 0' : (chx0) - bl = min { (c,,x°) _ bi (Cl, 8) (Cg, S) alli a! J with (c,,s) < 0}. Go to Step 2.4. Step 2.4 Set x0 := x0 - as and update J = {i1,-~- ,in} and D‘l. If I < my“, check if any lower edge {a,,1,,,a,.+,,,-} in F1,“ extends 1:1“. Collect those lower edges, if they exist, and delete them from PHI. Go to Step 2.2. 39 Remark 1 Numerical testing shows that setting up the inequality system (3.9) is very time consuming. One strategy we employ is to save the inequality systems at all previous levels. Thus the inequality system in the current level can be set up by using the inequality system that already exist. 3.3 Find All Fine Mixed Cells For S =(S1,---,S,,) with S,- = {a.-,1,~- ,a,,m,} C N", i = 1,--- ,n and generically chosen w = (w1,--- ,wn) with .63.}, i=1,...,.., s,- = {a = (a, ...,-(1.)) we combine our algorithms described in the last two sections in the following algorithm for finding all the mixed cells in S = (S1, - - - ,Sn). Algorithm 4 Find all mixed cells in S = (51,- -- ,5"). Step 0: Initialization. Find us.) for alli = 1,--- ,n. Set .71 1: 16(31), ’6 I: 1. Step 1: Backtracking. If k = 0 Stop. Iffih=0,set k:=k—1 andgo toStep 1. Otherwise go to Step 2. Step 2: Select next level-k subface to extend. Select e). E .73)., and set P). := P).\{é,.}. Let E), =(é1, - - - ,ék) and go to Step 3. Step 3: Extending the level-k subface. Find £(Ek). 40 If 8(3),) = 0, go to Step 1, otherwise set .73).“ = 8(Ek), k := k + 1 then go to Step 4. Step 4: Collect mixed cells. Ifk = n, all C' = (e1, - - - ,en_1,e),é E f}, are fine mixed cells, pick up all these mixed cells, then set k := k — 1, go to Step 1. Otherwise go to Step 2. Remark 2 In finding £(Ek) at Step 3, inequalities associated with the points in S,- which never appear in E (13),) for i = 1, - - - ,k - 1 should not be considered as constraints, since these points will never enter the level-k subface. 41 CHAPTER 4 Numerical Implementation 4.1 Algorithm For finding all isolated zeros of a polynomial system P(x) = (p1(x), - -- , p,,(x)) in C", where pi(x) = Z €5,3X‘, for i = 1: ' ' ' in: 365; we outline the major steps in brief terms as follows: (A) Set up the polyhedral homotopy Q(x, t) : C” x [0,1] ——> C" as q,(x, t) = Z (c,-,. + (1 — t)e,,,)t""(‘)x‘, for i = 1, ~ -~ ,n, a€S¢U{O} where w,- : S,- U {0} —> R are chosen generically and 6,3,. are randomly chosen complex numbers. (B) Find all mixed cells of extended support 51 U {0}, - - - , 3,, U {0}. For each inner normal (1 associated with a mixed cell, define the homotopy Ha(y,t) E t-5Q(yt°‘,t) = 0, where 6 = (61,- -- ,6“) and fl,- is the lowest order in t among all the terms in qi (yta: t) ° 42 (C) Solve the binomial system Ha(y, 0) = 0 in C", then follow homotopy paths of Ha(y,t) = 0 to find all the isolated zeros of P(x). We set up our initial system Q(x) = (q1(x), - -- ,q'n(x)) by perturbating the coef- ficients of P(x), that is 6.0:) = Z (61,. +e.-,.)X‘. 2': 1.--° .71 aESiU{0} where 5,3. 6 IR" are randomly chosen small complex numbers, and let Q(XJ) = (1 - t)Q(X) + tP(X)- The homotopy Q(x, t) in (A) is obtained by setting up the polyhedral homotopy for Q(x, t) instead with the same variable t. Apparently, we have Q(x, 1) = P(x) and for each t 6 [0,1], Q(x, t) and Q(x) have the same support. Let {1,-(x, t) = Z (c,-,. + c,,.)x‘t‘”‘(‘), i = 1, - -- ,n. .es.u{o} It is clear that, for any a 6 1R“, the lowest order terms in t of both qJ-(xt°‘,t) and q“,- (yta, t) are the same. Hence the mixed cells and their associated inner normals stay invariant, and the start system of the homotopy H.(y, t) = t""Q(yt°.t) = 0 (4.1) is the same as that of the Ra(y, t) E t’5Q(yt°‘,t) = 0 with Q(yt°,t) = (61(yt°.t),°-- .(in(yt°‘.t))- Here. again. 0 = (fi1,---,fi.) and for j = l.--- ,n. 6,- is the lowest order in t among all the terms in q,- (yta, t). Thus, when nonsingular solutions of Ho, (y, t) in (C‘ )” are available, we may follow those homotopy paths of Ha(y, t) in (4.1) instead with those starting points. 43 4.2 Implementation In this section, we will briefly describe our software environment that has been de- veloped. Currently there are several publically available software packages dedicated to solv- ing polynomial systems by homotopy continuations. HOMPACK [45] and CONSOL [31] are written in FORTRAN77, pss [28] and Pelican [16] are written in C and PHC [42] is written in Ada. Some of these software packages are integrated multi—purpose packages. Our package is designed as a high-performance polynomial system solver, focused on better eficiency, portability and simplicity. Our program is written in C++, a standard programming language which provides excellent support for ob ject-oriented programming, abstraction, and encapsulation. The following UML diagrams illustrate the structure of our polynomial system solver system. Diagram 4.2 shows that our polynomial system solver relys on two packages. One of them provides utility tools such as linear programming and lin- ear system solvers. The other one may not be critical, it is mainly used for parsing polynomial expressions and performing simple polynomial manipulation to provide certain degree of user friendly interface. The class diagram 4.2 shows our Polyno- mialSolver uses four major components: PolyhedralHomotopy, MixedCell, Binomial- SystemSolver, and Continuation (which uses Newton’s iteration). The state diagram 4.2 shows the transitions among states in the execution of our program. 4.3 Numerical Results Our software package has solved many well-known polynomial systems successfully. In this section we present some of our numerical results. 44 _—l Polynomial system ___] Linear Algebra Linear Programming T Polynomial system Solver Figure 4.1: Package Diagram PolynomialSolver load() run() 1 PolyhedralHomotopy BinomialSystemSoTver MixedCellFinder COBUTIMUOD getLiftedSupport() loadti first() operator()() qetBinamralSystem() first() next() evaluateValue() next() aetCe11() evaluateJocobianl) 1 . noralize(); aetSo utron() Figure 4.2: Class Diagram 45 1 Newton operator()() ? Setting Up Homotopy System next() Transforming next Homotopy [get next cell] FM Finding Mixed Cell] [no more cell] J 69 ‘ next r Solving [get the next solution] Binomial System Following Curve next() J L Figure 4.3: State Diagram 1. System of 'Ih'inks from the PoSSo test suite[40] f 45y + 35a — 165v — 36, 353/ + 252 + 40t — 27u, 25yu — 165v2 + 15x — 182 + 30t, P(xlyiztutvtt) =< l5yz + 20tu — 9x, —11v3 + xy + Zzt, —11uv + 31)2 + 99x, 46 2. Generalized eigenvalue problem[9] r L —10x1x§ + 2x2x§ -— x3x§ + x4x§ + 3x5x§ + xlxe + 2x2x6 + x3x6 +2x4x6 + x5x3 + 10x1 + 2x; — x3 + 2x; — 2x5, 2x1x§ - 11222:: + 2x3x§ -— 2x4x§ + 2:52: + 2x1x6 + xgxe + 2x3x6 +x4x5 + 3x5x6 + 2x1 + 9x2 + 3x3 - x4 — 2x5, —x1x§ + 2x2x§ — 12x3x§ — x4x§ + x5x§ + xlxa + 2x2x6 — 2x4x6 —2x5x5 — x1 + 3x; + 10x3 + 2x.; — x5, 2:133, — 2x2x§ — 23x: — 10x4x§ + 2x5x§ + 2x1x5 + xzxs — 2x3x6 +2x4x5 + 3x5x5 + 2x1 — x2 + 2x3 + 12x4 + x5, 3x1x§ + $225; + x3x§ + 2x4x§ — llxsxfi + $12.35 + 3x2x5 — 2x3x6 +3x4x5 + 3x5x5 — 2x1 — 2x2 — x3 + x4 + 10x5, $1+$2+$3+$4+$5—1, 3. A system from economics modeling [31] The following formulates the system for general dimension n. n—k-l xk+ xix”). x,,—c,,,k=1,---,n—1. P(X)= ( Z ) l=l 2?: $1 + I. The constant c), can be chosen at random. This problem has solved for dimension up 12. 47 4. Totally mixed Nash equilibria for 5 players with two pure strategies[29] r 1390350657193 + 0.4641393136113192 - 0.6266605268p2 — 1.400171891p4 —2.090683800p4p3p2 + 4.089263882p4p3 + 1.129827638124192 +1.881614464p5p4p3p2 + 0.4716169661p5p3p2 — 0.8625849122125154};2 —1.398871056P5P2 + 09599693844125 + 0.071402539712515. +0.1073802376p5p3 — 0.0259664538125114153 — 0.2067814278, —1.196136754p3 — 09249804195114 + 0.3188761009p4p3 - 1.045301323193451 —0.0306661782p1 + 0.5987012929p4p3p1 - 0.4448182692p4p1 —0.3908068031p5p4p3p1 - 1.212939725p5 + 2.586129779p5p4 4.1180169224125193 — 1.0515195071251941»... + 2.13497937515512315l 4.337061849125124). — 0.2961272671P5P1 + 0.7316111016, 22729431631;2 —- 04131564265124 — 1.920446680104122 + 000805092341;1 +1.342851102p4p1 — 2.979502184p2p1 + 3.3915718341541521;1 4.5975693742195114)». + 0.3002794716125192 — 0.7893445350p5 +1.276948001p5p4 — 4.601376311P5P4P1 + 23568043221551»1 +3.498840190p5p4p2p1 — 1.3553750151251521;1 — 1.231070236, —2.206336116p3 + 2.318673689p3p2 — 1.267478048p2 + 1.110654516p3p1 +1.533592098p1 — 1.872504375p2P1 + 0.3299103675p3p2p1 4.4007504721951931». + 2.0936745161551112 — 1.772874182p5 +2.993821915p5p3 — 1.356762392p5P3P1 + 0.0637534233p5p1 +0.5870371377P5P2P1 + 1.018269743P5P3P2P1 + 1.400431557, —2.522718869p3 + 0.8323646978p3p2 — 1.375030881p2 — 0.3055448755p4 +0.6760632172p4p3p2 -— 0.4262974456134193 + 1.268255245p3p1 +0.5352674901P4P2 — 1024495558121 + 1.8182754041941531;1 4.354832512194121 — 1.595112039122111 + 2.2379562421241521), +3.370102170p3p2p1 — 3.465040669p4p3p2p1 + 2.132631128, 48 5. Benchmark i1 from the Interval Arithmetics Benchmarks [15] x1 — 0.25428722 - 0.18324757x4x3x9, x2 — 0.37842197 - 0.16275449x1x10x5, x3 — 0.27162577 — 0.16955071x1x2x10, x4 —- 0.19807914 — 0.15585316x7x1x5, x5 — 0.44166728 — 0.19950920x7x5x3, x5 — 0.14654113 — 0.18922793x3x5x10, x7 — 0.4293716] — 0.21180484x2x5x3, x8 — 0.07056438 — 0.17081208x1x7x5, x9 — 0.34504906 — 0.19612740x10x6x3, [ x10 — 0.42651102 - 0.21466544x4x8x1, This system is very sparse, the mixed volume is much smaller than the total degree. 6. Cyclic-n problem [13] The general formulation goes as follows: a I: ZH$(i+j)modn, k: 1,... ,n-1 P(X) ____ i=1 j=l n 11%“ " 1, i=1 This system is widely considered as a benchmark problem for polynomial system solving. We have solved this system up to dimension 11. 49 7. The construction of Virasoro algebras[38] The test results are summarized in the following table. In the table, M(A°) denotes the mixed volume of the extended support A0 = (A1 U {0}, . . . , .A" U {0}). The number of zeros in the examples are obtained from the numerical results of our algorithm. It is well-known that homotopy curves may converge to solutions in an algebraic variety with nonzero dimension, i.e., they may lead to non-isolated zeros of the target polynomial systems. In our root count, we do not exclude those numerical solutions at which the J acobian matrices of the corresponding polynomial systems are almost singular. And our computation were carried out on a 400Mhz Intel Pentium f 8x? + 8x1x2 + 8x1x3 + 2x1x4 + 2x1x5 + 2x1x6 + 2x1x7 —8x2x3 — 2x4x7 — 2x5x6 — x1, 8x1x2 — 8x1x3 + 8x3 + 8x2x3 + 2x2x4 + 2x2x5 + 2x2x5 +2x2x7 — 2x4x3 — 2x5x7 — x2, —8x1x2 + 8x1x3 + 8x2x3 + 8x3 + 2x3x4 + 2x3x5 + 2x3x6 +2x3x7 — 2x4x5 — 2x6x7 - x3, 2x1x4 — 2x1x7 + 2x2x4 — 2x2x5 + 2x3x4 — 2x3x5 + 8x2 +8x4x5 + 2x4x3 + 2x4x7 + 6x4x3 — 6x5x3 — x4, 2x1x5 — 2x1x5 + 2x2x5 — 2x2x7 — 2x3x4 + 2x3x5 + 8x4x5 —6x4x3 + 8x3 + 2x5x5 + 2x5x7 + 6x5x3 — x5, —2x1x5 + 2x1x5 — 2x2x4 + 2x2x5 + 2x3x6 —- 2x3x7 + 2x4x6 +2x5x5 + 8x: + 8x6x7 + 6x5x3 - 6x7x3 - x5, —2x1x4 + 2x1x7 — 2x2x5 + 2x2x7 — 2x3x5 + 2x3x7 + 2x4x7 +2x5x7 + 8x5x7 — 6x6x3 + 8x-2, + 6x7x3 — x7, —6x4x5 + 6x4x8 + 6x5x3 — 6x5x7 + 6x3x8 + 6x7x3 + 8x3 — x3, II CPU with 256 MB of RAM, running SunOS 5.6. 50 System Total Degree M(A°) # of Zeros in C" CPU Time Trinks 24 10 10 290ms Eigenvalue 243 10 10 430ms Economics-13 354294 2048 2048 8m28200ms Economics-14 1062882 4096 4096 22m35s620ms Nash equilibia 1024 44 44 48240ms Benchmark i1 59049 66 50 48100ms Cyclic-10 3628800 35940 34940 lh33m46s Cyclic-11 39916800 184756 184756 8h55m36s Virasoro 256 256 256 378 690 ms Table 4.1: Numerical Results 51 BIBLIOGRAPHY 52 BIBLIOGRAPHY [1] E. L. Allgower and K. Georg (1990), Numerical Continuation Methods, an In- troduction, Springer Series in Comput. Math., Vol. 13, Springer-Verlag(Berlin Herdelberg, New York). [2] E. L. Allgower and K. Georg (1993), “Continuation and path following”, Acta [3] [4] [5] [6] [7] [8] [9] Numerica, 1-64. M. J. Best and K.Ritter (1985), Linear Programming: Active Set Analysis and Computer Programs, Prentice-Hall, Inc. New Jersey. D. N. Bernshtein (1975), “The number of roots of a system of equations”, Func- tional Analysis and Appl., 9(3), 183-185. Translated from Funktsional. Anal. i Prilozhen., 9(3), 1-4. B. Buchberger ( 1985), “Grfibner basis: An algorithmic method in polynomial ideal theory”, In Multidimensional System Theory (N.K. Bose, ed.), D. Reidel Publishing Company (Dordrecht Boston Lancaster), 184-232. J. Canny and J. M. Rojas (1991), “An optimal condition for determining the exact number of roots of a polynomial system”, In Proceedings of the 1991 In- ternational Symposium on Symbolic and Algebraic Computation, ACM, 96-101. S. N. Chow, J. Mallet-Paret, and J. A. Yorke (1978), “Finding zeros of maps: homotopy methods that are constructive with probability one”, Math. Comput., 32, 887-899. S. N. Chow, J. Mallet-Paret, and J. A. Yorke (1979), “Homotopy method for 10- cating all zeros of a system of polynomials”, In Functional difl'erential equations and approximation of fixed points (H. O. Peitgen and H. O. Walther, eds.), Lec- ture Notes in Mathematics Vol. 730, Springer-Verlag (Berlin Heidelberg), 77-88. M. Chu, T.-Y. Li and T. Sauer (1988) ”Homotopy method for general lambda- matrix problems”, SIAM J. Matrix Anal. Appl., vol. 9, No. 4, pp 528-536. 53 [10] H.H. Chung (1998), Polyhedral homotopy and its applications to polynomial sys- tem solving, Ph.D. dissertation, Michigan State University. [11] F. J. Drexler (1977), “Eine Methode zur Berechnung siimtlicher Léisungen von Polynomgleichungssystemen”, Numer. Math. 29, 45-58. [12] I. Emiris (1994), Sparse Elimination and Applications in Kinematics, Ph.D. the- sis, Computer Science Division, Dept. of Electrical Engineering and Computer Science, University of California (Berkeley). [13] I. Emiris and J. Canny (1995), “Efficient incremental algorithms for the sparse resultant and the mixed volume”, J. Symbolic Computation 20, 117-149. [14] C. B. Garcia and W. I. Zangwill (1979), “Finding all solutions to polynomial systems and other systems of equations”, Math. Programming, 16, 159-176. [15] P. Van Hentemyck, D. McAllester, and D. Kapur (1995), ”Solving polynomial systems using a branch and prune approach”. Technical Report No. CS-95-01, Department of Computer Science, Brown University. [16] B. Huber, “Pelican manual”, available via the author’s Web page. [17] B. Huber (1996), Solving sparse polynomial systems, Ph.D. dissertation, Cornell University. [18] B. Huber and B. Sturmfels (1995), “A polyhedral method for solving sparse polynomial systems”, Math. Camp, 64, 1541-1555. [19] B. Huber and B. Sturmfels (1997), “Bernstein’s theorem in afine space”, Discrete Comput. Geom., 7, no. 2, 137-141. [20] A. G. Khovanskii (1978), “Newton polyhedra and the genus of complete intersec- tions”, Functional Anal. Appl., 12(1), 38-46. Translated from Funktsional. Anal. i Prilozhen., 12(1), 51-61. [21] A. G. Kushnirenko (1976), “Newton Polytopes and the Bézout Theorem”, Functional Anal. Appl., 10(3), 233-235. Translated from Funktsional. Anal. i Prilozhen., 10(3), 82-83. [22] T. Y. Li (1983), “On Chow, Mallet-Paret and Yorke homotopy for solving systems of polynomials”, Bulletin of the Institute of Mathematics, Acad. Sin., 11, 433- 437. 54 [23] T. Y. Li (1997), “Numerical solution of multivariate polynomial systems by ho- motopy continuation methods”, Acta Numerica, 6, 399-436. [24] T. Y. Li and T. Sauer (1989), “A simple homotopy for solving deficient polyno- mial systems”, Japan J. Appl. Math., 6, 409-419. [25] T. Y. Li, T. Sauer and J. A. Yorke (1989), “The cheater’s homotopy: an eficient procedure for solving systems of polynomial equations”, SIAM J. Numer. Anal, 26, 1241-1251. [26] T. Y. Li and X. Wang (1992), “Nonlinear homotopies for solving deficient poly- nomial systems with parameters”, SIAM J. Numer. Anal, 29, 1104-1118. [27] T. Y. Li and X. Wang (1996), “The BKK root count in C"”, Math. Camp, 65, no. 216, 1477-1484. [28] G. L. Malajovich. pss 2. alpha, polynomial system solver, version 2. alpha. Dis- tributed by the author through gopher. [29] Richard D. McKelvey and Andrew McLennan (1997), ”The maximal number of regular totally mixed Nash equilibria”, Journal of Economic Theory, Volume 72, pages 411-425. [30] A. P. Morgan (1986), “A homotopy for solving polynomial systems”, Appl. Math. Comput., 18, 173-177. [31] A. P. Morgan (1987), Solving polynomial systems using continuation for engi- neering and scientific problems, Prentice-Hall (Englewood Cliffs, N. J) [32] A. P. Morgan and A. J. Sommese (1987), “A homotopy for solving general poly- nomial systems that respect m-homogeneous structures” , Appl. Math. Comput., 24, 101-113. [33] A. P. Morgan and A. J. Sommese (1989:1), “Coefficient-parameter polynomial continuation”, Appl. Math. Comput., 29, 123-160. Errata: Appl. Math. Comput., 51, 207 (1992). [34] AP. Morgan and CW. Wampler (1990), “Solving a planar four-bar design prob- lem using continuation”, ASME J. of Mechanical Design, 112, 544-550. [35] C. H. Papadimitriou and K. Steiglitz (1982) Combinatorial optimization: algo- rithms and complexity, Prentice-Hall, Inc. New Jersey. 55 [36] J. M. Rojas (1994), “A convex geometric approach to counting the roots of a polynomial system”, Theoret. Comput. Sci., 133, 105-140. [37] J. M. Rojas and X. Wang (1996), “Counting affine roots of polynomial systems via pointed Newton polytopes”, J. of Complexity, 12, no. 2, 116-133. [38] S. Schrans and and W. 'Droost (1990) ”Generalized Virasoro Constructions for SU(3)”, Nuclear Phys. B, Vol. 345, No. 2—3, pp. 584—606. [39] I. R. Shafarevich (1977), Basic Algebraic Geometry, Springer-Verlag (New York). [40] C. Traverso (1997), “The PoSSo test suite examples”, [Online] Available at http: / / www.inria.fr / safir / POL / index.html. [41] J. Verschelde and K. Gatermann (1995), “Symmetric Newton polytopes for solv- ing sparse polynomial systems”, Adv. Appl. Math., 16, 95-127. [42] J. Verschelde (1995), “PHC and MVC: two programs for solving polynomial systems by homotopy continuation”. Presented at the PoSSo workshop on soft- ware, Paris. Available by anonymous ftp to ftp.cs.kuleuven.ac.be in the directory / pub / NumAnal-Apleath / PHC. [43] J. Verschelde (1996), “Homotopy continuation methods for solving polynomial systems, Ph. D. thesis, Department of computer S cience, Katholieke Universiteit Leuven (Leuven, Belgium). [44] J. Verschelde, K. Gatermann, and R. Cools (1996), “Mixed-volume computation by dynamic lifting applied to polynomial system solving” , Discrete Comput. Geom., 16, no. 1, 69-112. [45] L. T. Watson, S. C. Billups, and A. P. Morgan. ”Algorithm 653: HOMPACK: a suite of codes for globbally convergent homotopy algorithms.” ACM Tiuns. Math. Softw., 13(3): 281-310, 1987. [46] A. H. Wright (1985), “Finding all solutions to a system of polynomial equations”, Math. Comput., 44, 125-133. [47] W. Zulener (1988), “A simple homotopy method for determining all isolated solutions to polynomial systems”, Math. Comput., 50, 167-177. 56 "I[illllflllllilllllllli