TORIC NEWTON METHOD FOR POLYNOMIAL HOMOTOPIES

JAN VERSCHELDE

Abstract. This paper defines a generalization of Newton’s method to deal with solution paths defined by polynomial homotopies that lead to extremal values. Embedding the solutions in a toric variety leads to explicit scaling relations between coefficients and solutions. Toric Newton is a symbolic-numeric algorithm where the symbolic pre-processing exploits the polyhedral structures. The numerical stage uses the additional variables introduced by the homogenization to scale the components of the solution vectors to the complex unit circle. Toric Newton generates appropriate affine charts and enables to approximate the magnitude of large solutions of polynomial systems. AMS Subject Classification. 14M25, 52B55, 65H10, 68Q40. Keywords. polynomial system, toric varieties, Newton’s method, homotopy continuation.

1. Introduction

If a polynomial system has a large root, we first want to estimate its magnitude before the approximation of its actual value. In the context of polyhedral homotopies we can certify [15] diverging solution paths and separate those from paths leading to extremal values. The purpose of this paper is to show how polyhedral methods can approximate the magnitude of large roots. In [14] and [26] polyhedral homotopies were introduced to exploit sparsity. The number of solution paths equals the bound of Bernshteǐn’s theorem. First we give some notations to state Bernshteǐn’s theorems.

Denote x^{a} = x_{1}^{a_{1}}x_{2}^{a_{2}} ⋯ x_{n}^{a_{n}}. The support of a polynomial f(x) = ∑_{a∈A} c_{a}x^{a} is the set A = { a ∈ ℤ^{n} | c_{a} ≠ 0 }. The Newton polytope P of f is the convex hull of its support A. The structure of a polynomial system F (x) = 0 is represented by the Newton polytopes P = (P_{1}, P_{2}, …, P_{n}) spanned by the supports A = (A_{1}, A_{2}, …, A_{n}). Denote ℂ^{∗} = ℂ \ {0}. In [1] Bernshteǐn proved that the number of isolated solutions in (ℂ^{∗})^{n} of a polynomial system with generic coefficients equals the mixed volume V_{n}(P) of its Newton polytopes P.

In the nongeneric case we examine faces of systems. A face ∂_{v}P of a polytope P is characterized by a direction v ∈ ℝ^{n} \{0} and spanned by ∂_{v}A = { a ∈ A | ∀b ∈ A : ⟨a, v⟩ ≤ ⟨b, v⟩ }. A face polynomial of f is ∂_{v}f(x) = ∑_{a∈∂_{v}A} c_{a}x^{a}. A face system of F is ∂_{v}F (x) = 0 with equations ∂_{v}f_{i}(x) = 0, i = 1, 2, …, n. For V_{n}(P) > 0, Bernshteǐn’s second theorem [1] states that F (x) = 0 has fewer than V_{n}(P) isolated solutions in (ℂ^{∗})^{n} if and only if there is some ∂_{v}F (x) = 0 that has a solution in (ℂ^{∗})^{n}.

A toric variety X is a solution set of a polynomial system that contains the torus (ℂ^{∗})^{n} as a dense open subset with (ℂ^{∗})^{N} as a group acting on X. In section 2 we make this definition from [11] concrete in applying the homogenization procedure described in [5]. The group action is used in the third section to scale the solution components to the complex unit circle, narrowing the torus (ℂ^{∗})^{n} to a numerically favorable region. See [6], [7], [11], [13], and [20] for more on toric varieties.

Research at MSRI is supported by a post-doctoral fellowship. Supported in part by NSF under Grant DMS - 9804846 at Michigan State University.

Representing diverging and badly-scaled solutions is a common struggle in the daily life of a numerical analysist. In [19, Chapter 5], [18], [27, 28] we find coefficient and equation scaling, combined with centering, an effective and useful tool. The classical homogenization [2], [17],

[19] in projective space easily leads to clustering and fails to recognize the deficiency behavior when applied to sparse polynomial systems. In the third section of this paper we exemplify this situation.

In the fourth section we define Toric Newton method and analyze the cost of its execution. This algorithm is a symbolic-numeric approach. The Toric Newton method consists of a symbolic pre-processing stage, where the homogeneous forms are determined and the normal fans are refined to be simplicial. This stage is performed only once. The numerical stage is executed iteratively while following the solution paths and takes place in the original space of dimension n.

Toric Newton augments the robustness of polyhedral end games [15] as used in numerical homotopies [17]. Convergence properties of Toric Newton are addressed in the fifth section. The complexity of Newton in the dense case has been analyzed in [2] and has recently been extended to sparse complexity models [8], multi-homogeneous [9] and over-determined systems [10]. The Toric Newton method parallels similar developments in providing algorithms using torus embeddings to compute sparse resultants [21], residues [3], [4], and Bezoutians [22]. All these methods rely on polytopes [29]. Here we give practical evidence for the performance of the polyhedral constructions.

2. Newton Polytopes and Homogeneous Coordinates

If the Newton polytope lies in a hyperplane, then the coordinates of the polynomial are homogeneous. The exponents of the scaler λ are given by the components of the direction that is perpendicular to the hyperplane, as expressed formally in the following proposition.

Proposition 2.1. ([13, page 193]) Consider f(x) = ∑_{a∈A} c_{a}x^{a}.

If ∃v ∈ ℤ^{n} \ {0}, ∃d ∈ ℤ: ∀a ∈ A, ⟨a, v⟩ = d, then f(λ^{v_{1}}x_{1}, λ^{v_{2}}x_{2}, …, λ^{v_{n}}x_{n}) = λ^{d}f(x).

Proof. Denote λ^{v}x = λ^{v_{1}}x_{1}λ^{v_{2}}x_{2} ⋯ λ^{v_{n}}x_{n}, then (λ^{v}x)^{a} = (λ^{v_{1}}x_{1})^{a_{1}}(λ^{v_{2}}x_{2})^{a_{2}} ⋯ (λ^{v_{n}}x_{n})^{a_{n}} = λ^{⟨a,v⟩}. So, f(λ^{v}x) = ∑_{a∈A} c_{a}(λ^{v}x)^{a} = ∑_{a∈A} c_{a}λ^{⟨a,v⟩}x^{a} = λ^{d}f(x). □

The multilinearity of the inner product provides us with a whole range of homogeneous coordinates when the polytope lies in the intersection of many hyperplanes.

Example 2.1. (projective space) Let (x, y, z) ∈ ℙ^{2} and consider f(x, y, z) = x + y + z. The support of f is the set A = {(1, 0, 0), (0, 1, 0), (0, 0, 1)} and lies in the hyperplane ⟨a, v⟩ = d with v = (1, 1, 1), d = 1. The homogeneous coordinates are (λ^{1}x, λ^{1}y, λ^{1}z).

Example 2.2. (multi-projective space) Let (x_{0}, x_{1}, y_{0}, y_{1}) ∈ ℙ^{1} × ℙ^{1} and consider f(x_{0}, x_{1}, y_{0}, y_{1}) = x_{1}y_{1} − x_{0}y_{0}. The support of f is the set A = {(0, 1, 0, 1), (1, 0, 1, 0)} and lies in the hyperplanes ⟨a, v^{i}⟩ = d with v^{1} = (1, 1, 0, 0), v^{2} = (0, 0, 1, 1), v^{3} = (1, 0, 0, 1), and d = 1. First we have the two usual bi-homogeneous coordinates (λ^{1}x_{0}, λ^{1}x_{1}, y_{0}, y_{1}) and (x_{0}, x_{1}, λ^{1}y_{0}, λ^{1}y_{1}), which allow to scale the x-variables independently from the y’s. In addition, we have the third relation, because the Newton polytope is an edge and not of dimension two as one would expect. The corresponding homogeneous coordinates are (λ^{1}x_{0}, x_{1}, y_{0}, λ^{1}y_{1}).

Examples 2.1 and 2.2 are particular cases of homogeneous coordinates. The general definition is given below. In applying the homogenization method in [5] we copy the original monomials x^{a}. By this adaption we recognize the structure of the polytope and have explicit scaling relations between the coefficients and solutions, as needed in the toric version of Newton’s method.

Definition 2.1. Consider f(x) = ∑_{a∈A} c_{a}x^{a}, with P = conv(A). Assume vol(P) > 0. Let

P have N facets given by N inward pointing normals v^{i} ∈ ℤ^{n} \ {0} with all components relatively prime. For every v^{i}, define d_{i} = min_{a∈A}⟨a, v^{i}⟩. The facets are spanned by ∂_{v^{i}}A = {a ∈ A|⟨a, v^{i}⟩ = d_{i}}. The homogenization map ϕ_{A} introduces for every v^{i} a new variable z_{i}, as follows:

ϕ_{A}(x^{a}) = x^{a} ∏_{i=1}^{N} z_{i}^{⟨a,v^{i}⟩−d_{i}} and ϕ_{A} (∑_{a∈A} c_{a}x^{a}) = ∑_{a∈A} c_{a}ϕ_{A}(x^{a}).  (1)

If vol(P) = 0, then A lies in some hyperplanes ⟨x, v^{i}⟩ = c_{i}, and we add ±v^{i} to the normals of the k-faces of P, where k = dim(A) − 1.

The result of the homogenization is the map ϕ_{A}(f) : (ℂ^{∗})^{n} × ℂ^{N} → ℂ : (x, z) ↦ ϕ_{A}(f)(x, z). We may allow zero values for the z-variables, because of the following property.

Proposition 2.2. Let deg(f, x) denote the maximal exponent of the variable x in the polynomial f. Then,

deg (ϕ_{A} (∑_{a∈A} c_{a}x^{a}) , z_{i}) ≥ 0 and deg(ϕ_{A}(x^{a}), z_{i}) \{ \begin{matrix}= & 0 & if a ∈ ∂_{v^{i}}A \\ > & 0 & if a ∉ ∂_{v^{i}}A \end{matrix} . (2)

Proof. The directions perpendicular to the facets are given by inward pointing normals, i.e.: for a facet normal v^{i}: d_{i} ≤ ⟨a, v^{i}⟩, ∀a ∈ A. Therefore all exponents of z_{i} of ϕ_{A} are nonnegative. □

This homogenization corresponds geometrically to assigning as exponents of z_{i} the lattice distances relative to the ith facet of P. If we collect for every monomial the exponents of the new z-variables and store these vectors rowwise in a matrix, then we obtain the vertex-facet incident matrix [29] of the polytope P. The normal fan Σ_{P} collects all cones that are normal to the faces of P. The generators of those normal cones are the normals to the facets.

Example 2.3. (running example) Consider f(x_{1}, x_{2}) = c_{10}x_{1} + c_{01}x_{2} + c_{20}x_{1}^{2} + c_{12}x_{1}x_{2}^{2}. The Newton polytope has four facets, with inward pointing normals given by v^{1} = (1, 1), v^{2} = (0, 1), v^{3} = (−2, −1), and v^{4} = (1, −1). In Figure 1, the homogenizing along v^{4} is illustrated.

The homogenization is ϕ_{A}(f)(x, z) = c_{10}z_{3}^{2}z_{4}^{2}x_{1} + c_{01}z_{2}z_{3}^{2}x_{2} + c_{20}z_{1}z_{4}^{3}x_{1}^{2} + c_{12}z_{1}^{2}z_{2}^{2}x_{1}x_{2}^{2}.

The homogeneous coordinates correspond to the inward pointing normals to the facets.

[figure: z_{4}; z_{1}; z_{2}; z_{3}; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; 0; 1; 2; 3; P; v^{1}; v^{2}; v^{3}; v^{4}; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; Σ_{P}]

Figure 1. On the left we see that the exponents of z_{4} are lattice distances relative to the fourth facet of P with normal v^{4}. On the right the normal fan Σ_{P} of P is displayed.

Proposition 2.3. Consider the facet normal v^{i} with corresponding d_{i}. Denote λ^{v^{i}}xz = (λ^{v_{1}^{i}}x_{1}, λ^{v_{2}^{i}}x_{2}, …, λ^{v_{n}^{i}}x_{n}, z_{1}, …, z_{i−1}, λ^{−1}z_{i}, z_{i+1}, …, z_{n}). Then ϕ_{A}(f)(λ^{v^{i}}xz) = λ^{d_{i}}ϕ_{A}(f)(x, z).

Proof. First consider one monomial x^{a}:

ϕ_{A}(x^{a}) |_{(x,z)=λ^{v^{i}}xz} = x^{a} ∏_{j=1}^{N} z_{j}^{⟨a,v^{j}⟩−d_{j}} |_{(x,z)=λ^{v^{i}}xz}  (3)

= (λ^{v_{1}^{i}}x_{1}λ^{v_{2}^{i}}x_{2} ⋯ λ^{v_{n}^{i}}x_{n})^{a} ( ∏_{j=1,≠i}^{N} z_{j}^{⟨a,v^{j}⟩−d_{j}}) (λ^{−1}z_{i})^{⟨a,v^{i}⟩−d_{i}} (4)

= λ^{⟨a,v^{i}⟩}x^{a}λ^{−⟨a,v^{i}⟩+d_{i}} ∏_{j=1}^{N} z_{j}^{⟨a,v^{j}⟩−d_{j}}  (5)
= λ^{d_{i}}ϕ_{A}(x^{a})  (6)

The multilinearity of ϕ_{A} completes the proof. □

The transition from toric to affine space is immediate by taking λ = z_{i} in every λ^{v^{i}}xz.

Example 2.4. (Example 2.3 continued.) The homogeneous coordinates corresponding to the facets are given by

\begin{matrix}( & λ^{1}x_{1} & , & λ^{1}x_{2} & , & λ^{−1}z_{1} & , & z_{2} & , & z_{3} & , & z_{4} & ) \\ ( & λ^{0}x_{1} & , & λ^{1}x_{2} & , & z_{1} & , & λ^{−1}z_{2} & , & z_{3} & , & z_{4} & ) \\ ( & λ^{−2}x_{1} & , & λ^{−1}x_{2} & , & z_{1} & , & z_{2} & , & λ^{−1}z_{3} & , & z_{4} & ) \\ ( & λ^{1}x_{1} & , & λ^{−1}x_{2} & , & z_{1} & , & z_{2} & , & z_{3} & , & λ^{−1}z_{4} & ) \end{matrix}  (7)

for the respective normals v^{i}, i = 1, 2, 3, 4. Observe that as λ → 0, the x-values in affine space go either to zero or infinity. We have a multitude of homogeneous coordinates by taking combinations, see Σ_{P} in Figure 1. Positive combinations of any two vectors normal to adjacent facets correspond to constructing a vector inside the cone σ ∈ Σ_{P} that has these vectors among its generators. The relations (7) illustrate how the group (ℂ^{∗})^{N} acts on the variety.

For general polynomial systems not all supports are equal. Suppose there are r different supports which we group into A = (A_{1}, A_{2}, …, A_{r}) spanning the Newton polytopes P = (P_{1}, P_{2}, …, P_{r}). The idea is to embed the system using r distinct sets of homogenization variables, which provides a more cost-effective generalization than taking the Minkowski sum of the polytopes in P. This idea works because any normal to the Minkowski sum is a sum of normals to the components in the sum, whence we can capture all faces in this way. Definition 2.1 is extended to systems with different supports as follows, taking r = n for notational convenience.

Definition 2.2. Let A be the supports of F (x) = 0. Suppose the polytope spanned by A_{i} has N_{i} facet normals. Abbreviate z = (z^{1}, z^{2}, …, z^{n}). The mixed homogenization map ϕ_{A} is

ϕ_{A}(F) : (ℂ^{∗})^{n} × ℂ^{N_{1}} × ℂ^{N_{2}} × ⋯ × ℂ^{N_{n}} → ℂ^{n} : (x, z) ↦ ϕ_{A}(F)(x, z),  (8)

with ϕ_{A}(F)(x, z) = (ϕ_{A_{1}}(f_{1})(x, z^{1}), ϕ_{A_{2}}(f_{2})(x, z^{2}), …, ϕ_{A_{n}}(f_{n})(x, z^{n})) and ϕ_{A_{i}}(f_{i}) as in (1). For a homotopy H(x, t) = 0, the continuation parameter t is treated as a constant in setting up the homogenized homotopy ϕ_{A}(H)(x, z, t) = 0.

For several polytopes sharing the same facet normal, we can use the same corresponding z-variable. With this economical use of z-variables, Definition 2.2 yields the same homogenization map ϕ_{A} as in Definition 2.1 when applied to A = (A, A, …, A). Since the general homogenization is also determined by inward pointing normals v^{i}, Propositions 2.2 and 2.3 have a direct generalization to the case of different supports and are omitted.

Proposition 2.4. Suppose A is the support of f(x). Let (x^{∗}, z^{∗}) ∈ (ℂ^{∗})^{n} × ℂ^{N} : ϕ_{A}(f)(x^{∗}, z^{∗}) = 0.

1. If z_{i}^{∗} = 0 and z_{j}^{∗} ≠ 0, ∀j ≠ i, then ∂_{v^{i}}f(x) = 0 has a solution in (ℂ^{∗})^{n}.

2. Let I = { i | z_{i}^{∗} = 0, #(∩∂_{v^{i}}A) ≥ 2 }, z_{j}^{∗} ≠ 0, ∀j ∉ I, and v = ∑_{i∈I} v^{i}.

Then ∂_{v}f(x) = 0 has a solution in (ℂ^{∗})^{n}. The solutions x^{∗} of ∂_{v}f correspond to solutions of f with components equal to 0 and/or ∞.

Proof. Although the second statement generalizes the first, we distinguish two cases:

1. Suppose there is only one z_{i}^{∗} = 0 with corresponding facet normal v^{i}. Then all monomials in ϕ_{A}(f) whose exponent vector does not lie on the ith facet vanish when z_{i} = z_{i}^{∗} = 0. For z_{j}^{∗} ≠ 0, ∀j ≠ i, the homogeneous coordinates λ^{v^{j}}x^{∗}z^{∗} (as in Proposition 2.3) are used successively, each time setting λ = z_{j}^{∗}, ∀j ≠ i, to obtain a solution in (ℂ^{∗})^{n} of ∂_{v^{i}}f(x) = 0.

2. When several z_{i}^{∗}’s are zero, more monomials in the system vanish. The condition on I guarantees that there are still enough monomials left in ∂_{v}f to admit a solution in (ℂ^{∗})^{n}. Successive application of Proposition 2.3 setting λ = z_{j}^{∗}, ∀j ∉ I, yields a solution in (ℂ^{∗})^{n} of ∂_{v}f(x) = 0.

The homogeneous coordinates contain as exponents for λ the components of the normals. Any normal v has at least one component different from zero, call it v_{k}. Setting λ = z_{i} → 0 to transform to affine space will lead to x_{k} → ∞, if v_{k} < 0, or x_{k} → 0, in case v_{k} > 0. □

In [5], solutions that have no z-variables equal to zero, or only in cases covered by Proposition 2.4 are called nontrivial. Indeed, without the conditions on I in Proposition 2.4, all monomials could vanish. The generalization of Proposition 2.4 to polynomial systems is immediate.

Proposition 2.5. Let F (x) = 0 have supports A. The nontrivial solutions of ϕ_{A}(F)(x, z) = 0 with some z-variables equal to zero correspond to solutions of face systems.

We close this section with some remarks on systems with different supports. In the homogenization we care for all coefficients in the system, therefore we distinguish between supports, rather than between polytopes. Secondly, for sparse polynomials it often happens that the supports lie in lower-dimensional planes. In that case we treat the normal vectors that define the hyperplanes that contain the supports as facet normals in addition to the normal to the faces that define the lower-dimensional polytopes. Finally, aside from this complication, working with different supports does not alter the homogenization method.

3. Scaling Coefficients and Solutions

The variables introduced by the homogenization allow to scale ϕ_{A}(f)(x, z) = 0, so that every component of x lies on the unit circle in ℂ^{∗}, while keeping the z-variables bounded.

Proposition 3.1. Suppose x ∈ (ℂ^{∗})^{n} satisfies ϕ_{A}(f)(x, z) = 0 for fixed z with 0 < z_{j} ≤ 1, ∀j = 1, 2, …, N. Let v = (− log |x_{1}|, − log |x_{2}|, …, − log |x_{n}|), with log of base 10. Then

1. ∃c ∈ ℝ_{+}^{n}, ∃J = (j_{1}, j_{2}, …, j_{n}), j_{i} ∈ {1, 2, …, N} : v = ∑_{i=1}^{n} c_{i}v^{j_{i}}.

2. The homogeneous coordinates

(λ^{(− ∑_{i=1}^{n} c_{i}v^{j_{i}})^{k}}x_{k}|_{k=1}^{n} , …, λ^{c_{i}}z_{j_{i}}|_{i=1}^{n} , …)  (9)

with λ = 0.1 scale (x, z) to (x̃, z̃) with |x̃_{k}| = 1 and 0 < z̃_{j} ≤ 1.

Proof. The above statements are true because:

1. The normal fan Σ_{P} of P is complete. Whence every vector v can be written as a positive combination of the generators of a cone σ ∈ Σ_{P}.

2. (− ∑_{i=1}^{n} c_{i}v^{j_{i}})_{k} = log |x_{k}|, so x̃_{k} = (0.1)^{log |x_{k}|}x_{k} = \frac{x_{k}}{|x_{k}|}, whence |x̃_{k}| = 1.

0 < z_{j_{i}} ≤ 1 ⇒ 0 < z̃_{j_{i}} = (0.1)^{c_{i}}z_{j_{i}} ≤ 1, as c_{i} ≥ 0. Thus 0 < z̃_{j} ≤ 1. □

Note that any base b can be used in the logarithm, it suffices to replace 0.1 in Proposition 3.1 by b^{−1}.

Example 3.1. (toric scaling analogue to projective space) Suppose the Newton polytope is a simplex. The homogeneous coordinates are

\begin{matrix}( & λ^{1}x_{1} & , & λ^{0}x_{2} & , & λ^{−1}z_{1} & , & z_{2} & , & z_{3} & ) \\ ( & λ^{0}x_{1} & , & λ^{1}x_{2} & , & z_{1} & , & λ^{−1}z_{2} & , & z_{3} & ) \\ ( & λ^{−1}x_{1} & , & λ^{−1}x_{2} & , & z_{1} & , & z_{2} & , & λ^{−1}z_{3} & ) \end{matrix}  (10)

for the respective facet normals v^{1} = (1, 0), v^{2} = (0, 1), and v^{3} = (−1, −1). Since this is a toric analogue, we are concerned about solutions with zero components and have more z-variables than in the case of ordinary homogenization. We consider two situations:

1. (x_{1}, x_{2}) = (2, 2), v = (− log(2), − log(2)), J = (3, 3), c = (log(2), 0)

⇒ z_{3} = (0.1)^{log(2)}z_{3} = \frac{1}{2}z_{3}.

2. (x_{1}, x_{2}) = (2, 1), v = (− log(2), − log(1)), J = (2, 3), c = (log(2), log(2))

⇒ z_{2} = (0.1)^{log(2)}z_{2} = \frac{1}{2}z_{2} and z_{3} = (0.1)^{log(2)}z_{3} = \frac{1}{2}z_{3}. In the first case it suffices to use one of the homogeneous coordinates, while in the second case a combination is constructed.

Example 3.2. (Example 2.3 continued.) Here we also consider two situations:

1. (x_{1}, x_{2}) = (\frac{1}{2}, 2), v = (− log(\frac{1}{2}), − log(2)), J = (4, 4), c = (log(2), 0)

⇒ z_{4} = (0.1)^{log(2)}z_{4} = \frac{1}{2}z_{4}.

2. (x_{1}, x_{2}) = (2, 2), v = (− log(2), − log(1)), J = (2, 3), c = (log(2), log(2))

⇒ z_{3} = (0.1)^{2 log(2)}z_{3} = \frac{1}{4}z_{3} and z_{4} = (0.1)^{3 log(2)}z_{4} = \frac{1}{8}z_{4}. The scaling relation is constructed by decomposing v in Σ_{P}, see Figure 1.

In the general case when the system has an r-tuple of different supports A, we decompose v = − log |x| in every component of Σ_{P} = (Σ_{P_{1}}, Σ_{P_{2}}, …, Σ_{P_{r}}) to compute new values of the z-variables which are distinct for every different support.

Propositions 2.5 and 3.1 allow us to give a concise proof of Bernshteǐn’s second theorem formulated below as Theorem 3.1 using a homotopy argument.

Definition 3.1. We say that the homotopy H(x, t) = 0 is in generic position with respect to the supports A that span the Newton polytopes P, if ∀t ∈ [0, 1): H(x, t) = 0 has supports A and has exactly V_{n}(P) regular solutions in (ℂ^{∗})^{n}.

Given a system F (x) = 0, the canonical example of a homotopy in generic position is

H(x, t) = (1 − t)F^{(0)}(x) + tF (x) = 0, t ∈ [0, 1],  (11)

where the system F^{(0)}(x) = 0 has random complex coefficients and the same supports as F.

Theorem 3.1. Consider F (x) = 0, with supports A, Newton polytopes P and V_{n}(P) ≠ 0.

#{ m_{x} | F (x) = 0, x isolated } < V_{n}(P) ⇔ ∃v ∈ ℝ^{n} \ {0}, ∃x ∈ (ℂ^{∗})^{n} : ∂_{v}F (x) = 0

where m_{x} denotes the multiplicity of the isolated solution x of F (x) = 0.

Proof. Construct any homotopy H(x, t) = 0 that is in generic position with respect to P and for which H(x, t = 1) = F (x). Consider the homogenized homotopy ϕ_{A}(H)(x, z, t) = 0. ⇒ If F (x) = 0 has fewer than V_{n}(P) isolated solutions in (ℂ^{∗})^{n}, then, as t → 1, either

some solution paths must diverge outside (ℂ^{∗})^{n} or some solutions paths must go to a component of solutions, because the homotopy is in generic position. In the first case, by Proposition 3.1, the only way to diverge is that some z_{i} → 0. Then, according to Proposition 2.5, we have a solution in (ℂ^{∗})^{n} of some face system. In the second case, the solution component allows at least one free variable to choose in moving on that component. So we can move a free variable outside (ℂ^{∗})^{n} for which some z_{i} is zero. Then we can again apply Proposition 2.5.

⇐ If there is a face system that has a solution x ∈ (ℂ^{∗})^{n}, then in the corresponding tuple of z-variables for which ϕ_{A}(H)(x, z, 1) = 0, ∃i: z_{i} = 0. Since the homotopy is in generic position, there is at least one solution path that ended at (x, z). By Proposition 2.5, the transition to affine coordinates shows that the solution path that ended at this couple (x, z) diverged to a solution that has components at infinity and/or zero components. So F (x) = 0 must have fewer than V_{n}(P) isolated solutions in (ℂ^{∗})^{n}. □

4. The Static Scaling in Projective Newton

Let us recall Newton’s method in projective space (see [17]). We assume we have a system H = (h_{1}, h_{2}, …, h_{n}) of homogeneous polynomials: h_{i}(λz) = λ^{deg(h_{i})}h_{i}(z) for i = 1, 2, …, n with λz = (λz_{0}, λz_{1}, …, λz_{n}) the usual homogeneous coordinates. For a solution at infinity: z_{0} = 0.

Algorithm 4.1. (Newton’s method in projective space)

Input: (H, z).  (Homogeneous system, approximate root)
Output: z.  refined approximation
i : |z_{i}| = max_{j=1}^{n} |z_{j}|;  maximal component
z := (\underbrace{\frac{z_{0}}{z_{i}} , …, \frac{z_{i−1}}{z_{i}}}, 1, \underbrace{\frac{z_{i+1}}{z_{i}} , …, \frac{z_{n}}{z_{i}}});  scale the approximation
x := ( \overbrace{z_{1}, …, z_{i−1}} , \overbrace{z_{i+1}, …, z_{n}} );  switch to an affine chart
F := H|_{z_{i}=1};  select evaluation map
x := x − D_{x}^{−1}F (x)F (x);  affine Newton method
z := (x_{1}, …, x_{i}, 1, x_{i+1}, …, x_{n}).  refined approximation

In the context of path following, H is a parametrized family of systems.

Example 4.1. (projective Newton on a sparse system) Consider the application of projective Newton on a sparse system of polynomials. Suppose the supports of the polynomials are all equal. We see the Newton polytope displayed in Figure 2. In the homogenization the simplex is implicitly taken as Newton polytope and to every facet we associate a new variable z_{i}. Since in projective space one does not distinguish solutions with zero components, the original monomials are not copied.

[figure: .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; 0; 1; 2; 3; 1; 2; 3; z_{0}; z_{2}; z_{1}; P]  f(x) = c_{10}x_{1} + c_{01}x_{2} + c_{20}x_{1}^{2} + c_{12}x_{1}x_{2}^{2} f̃(z) = c_{10}z_{0}^{2}z_{1} + c_{01}z_{0}^{2}z_{2} + c_{20}z_{0}z_{1}^{2} + c_{12}z_{1}z_{2}^{2}

Figure 2. On the left we see how the polytope P of f is approximated by a simplex which corresponds to the usual projective homogenization f̃ given on the right.

In case of divergence, z_{1} and/or z_{2} go to infinity. The scaling step in Algorithm 4.1 selects i = 1 or i = 2 and z_{0} strives to zero. In that case we see in Figure 2 that the monomial z_{1}z_{2}^{2} dominates. However, a diverging path cannot go to a monomial and goes to one of the other face systems. The multi-homogeneous version of Newton’s method will make the same monomial dominant.

Example 4.2. (projective Newton leads to clustering) Consider the following simple system:

\{ \begin{matrix}x_{1} − 10^{8} = 0 \\ x_{2}^{2} − 1 = 0 \end{matrix}  (12)

In affine space, the solutions (x_{1}, x_{2}) = (10^{8}, ±1) are well-separated. In projective space, they are represented as (z_{0}, z_{1}, z_{2}) = (1, 10^{8}, ±1) and scaled to (z_{0}, z_{1}, z_{2}) = (10^{−8}, 1, ±10^{−8}), which is a pair of numerically clustered solutions. Although this example is cooked up for this special occasion, this phenomenon is observed as typical in many applications. As the sum of the two Newton polytopes of the two polynomials is a rectangle, the toric homogenization corresponds to the bi-homogenization, which allows to scale the variables independently: (x_{1}, x_{2}) = (10^{8}, ±1) ↦ (z_{10}, z_{11}, z_{20}, z_{21}) = (10^{−8}, 1, 1, ±1), avoiding the clustering of solutions.

The most important difference with multi-projective Newton is the dynamic use of the scaling procedure of Proposition 3.1 which guarantees distinct paths, as formulated below in Theorem 6.1.

5. Toric Newton scales in a Dynamic Way

The set up of Toric Newton enables an efficient numerical scaling while path following.

Algorithm 5.1. (Toric Newton – symbolic set-up)

Input: (A, △).  (support , triangulation)
Output: H;  homogeneous coordinates matrix
Σ_{P}.  simplicial normal fan
1. V := Facet Normals(A, △);  extract facet normals
2. H := Homogenize(A, V);  relative lattice distances
3. Σ_{P} := Cone Subdivision(H, A, V).  refine normal fan

Next we analyze the cost of this algorithm. The dimension n and the cardinality #A measure the input size. The complexity of the output is bounded from below by the number of cells #△ in any triangulation △ of A. Recall that we are interested in all isolated solutions and to count the number of roots we compute volumes. So #△ is a lower bound for the number of isolated solutions.

Proposition 5.1. (Cost of Step 1 in Algorithm 5.1.) Consider a triangulation △ of the polytope P spanned by A. Then all facets of P can be enumerated in O(#△n(n^{3} + #A)) operations.

Proof. The candidate facets are the facets of the simplices in the triangulation. Every simplex has n + 1 facets. The computation of each candidate requires O(n^{3}) arithmetical operations. To test whether this candidate normal constitutes a facet normal we have to compute inner products with the points in A and look whether the inner product is minimal for points on the facet of the simplex. This gives the bound O(#△n(n^{3} + #A)). □

The number of arithmetic operations of Step 2 of Algorithm 5.1 equals O(n#V). Step 3 is only needed when the polytope is not simple, i.e.: when some vertices are contained in more than n facets, see the top of the polytope in Figure 3. In Step 3 we make the fan simplicial subdividing the normal cones. One method to subdivide normal cones is to cut corners (see [13, page 190] and [29, pages 54–55]) at vertices that have more than n facets, as illustrated in Figure 3. We refer to [11, Section 2.6] for more on refinements of fans in connection with toric varieties.

For a normal cone σ ∈ Σ_{P} at the vertex a we take the average \overline{v} of the facet normals that generate σ. Determine c such that the hyperplane ⟨x, \overline{v}⟩ = c separates a from the rest of the vertices that span the polytope. The normal cones τ to the intersection points of the hyperplane with the edges of the polytope refine the normal cone σ. The recursive structure of the normal fan Σ_{P} is illustrated at the right of Figure 3. In the decomposition of a vector v in a subdivision of σ, we replace the average vector \overline{v} in terms of the other generators of σ, so that we obtain a decomposition that uses only facet normals.

[figure: a_{3}; a_{2}; a_{1}; 000; 020; 200; 220; 111; b^{1}; b^{2}; b^{3}; b^{4}; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; .; A = [ \begin{matrix}1 1 1 \\ 0 0 0 \\ 2 0 0 \\ 0 2 0 \\ 2 2 0 \end{matrix} ]; V = [ \begin{matrix}0 0 1 \\ 0 1 − 1 \\ −1 0 − 1 \\ 1 0 − 1 \\ 0 − 1 − 1 \end{matrix} ]; v^{1}; v^{2}; v^{3}; v^{4}; v^{5}; H = [ \begin{matrix}1 0 0 0 0 \\ 0 0 2 0 2 \\ 0 0 0 2 2 \\ 0 2 2 0 0 \\ 0 2 0 2 0 \end{matrix} ]; Σ_{P} = [ \begin{matrix}1 1 1 \\ 0 0 0 \\ 2 0 0 \\ 0 2 0 \\ 2 2 0 \end{matrix} ]; (v^{1}, v^{3}, v^{5}); (v^{1}, v^{4}, v^{5}); (v^{1}, v^{2}, v^{3}); (v^{1}, v^{2}, v^{4}); ^{[ \begin{matrix}b^{1} \\ b^{2} \\ b^{3} \\ b^{4} \end{matrix} ]}; (\overline{v}, v^{2}, v^{4}); (\overline{v}, v^{2}, v^{3}); (\overline{v}, v^{4}, v^{5}); (\overline{v}, v^{3}, v^{5})]

Figure 3. Refinement of a normal cone of a nonsimple polytope by cutting the corner at (1, 1, 1). The points b^{i} lie on the intersection of the hyperplane cut with the edges that contain (1, 1, 1).

In three dimensions we use only one separating hyperplane, because every edge is enclosed in exactly two facets. This is no longer the case in higher dimensions, so we proceed recursively as the points in the intersection of the polytope with the separating hyperplane belong to a space of lower dimension. By Proposition 2.4 from [29, page 55], edges are in oneto-one correspondence with vertices of the polytope cut out by the separating hyperplane. Recursive application proves the following.

Proposition 5.2. (Cost of Step 3 in Algorithm 5.1.) To triangulate an n-dimensional cone at most n − 2 separating hyperplanes are needed.

Although the number of cones grows exponentially, one new refining cone requires only elementary operations to compute the intersection points of the edges with the separating hyperplane. Therefore we can perform a lazy refinement of Σ_{P}, decomposing only those normal cones when needed in the scaling procedure.

Just as projective Newton method, we switch to an appropriate affine chart of dimension n.

Algorithm 5.2. (Toric Newton – number crunching)

Input: ϕ = ϕ_{A}(F);  system in homogeneous coordinates
(x, z);  approximate root with |x_{i}| = 1 and |z_{j}| ≤ 1
Σ_{P}.  simplicial normal fan
Output: (x, z).  refined approximation with |x_{i}| = 1 and |z_{j}| ≤ 1

F := ϕ|_{F ix z};  evaluation map
x := x − D_{x}^{−1}F (x)F (x);  affine Newton method
v := − log |x|;  direction of correction
J = (j_{1}, j_{2}, …, j_{n})  find σ ∈ Σ_{P} : v ∈ σ
:= Localize Position(Σ_{P} , v);  σ = span(v^{j_{1}}, v^{j_{2}}, …, v^{j_{n}})
c = (c_{1}, c_{2}, …, c_{n}) : v = ∑_{i=1}^{n} c_{i}v^{j_{i}};  decompose v

for all i ∈ {1, 2, …, n} do

x_{i} := \frac{x_{i}}{|x_{i}|}; z_{j_{i}} := (0.1)^{c_{i}}z_{j_{i}};  scale solution and coefficients

end for.

In analyzing the cost of localizing the position of the direction vector v we assume that the fan Σ_{P} has a recursive simplicial structure, as illustrated in Figure 3.

Proposition 5.3. (Cost of Localize Position in Algorithm 5.2.) For v ∈ ℝ^{n} and Σ_{P}, P = conv(A), finding σ ∈ Σ_{P} so that v ∈ σ takes at least #A and at most #A(n − 1) inner products.

Proof. If P is simple, then σ is the normal cone at the point a ∈ A for which ⟨a, v⟩ = min_{b∈A}⟨b, v⟩. Computing this point requires #A inner products. If P is not simple, there are at most n − 2 refinement levels, according to Proposition 5.2. At every level, there are at most #A points to consider. □

For r-tuples of supports A, we construct Σ_{P} = (Σ_{P_{1}}, Σ_{P_{2}}, …, Σ_{P_{r}}). The cost of the geometric constructions is scaled by a factor r. From our experiences with mixed-volume computations [25] we find linear upscaling in the cost rather remarkable.

6. Numerical Convergence Properties

In the absence of singularities, Newton’s method converges quadratically, doubling the precision of the approximate root in each step. With a homotopy in generic position, we will not hit a solution component, except perhaps only at the very end of the paths. What may happen however, is that paths become clustered, hereby not only limiting the convergence but possibly also terminating the numerical path tracking with a failure to reach the end.

Affine path following methods deal well with paths diverging outside (ℂ^{∗})^{n}, for those paths exhibit a steady direction leading to the face systems that cause the deficiency. To obtain a certificate of divergence, polyhedral end games [15] can effectively trace the paths in affine space, because to compute the direction of divergence one does not need to follow the paths up to infinity. Reaching badly-scaled though probably still meaningful solutions by affine path following methods is often impossible due to the finiteness of numbers a computer can represent. Projective Newton brings infinity to our working space, but roots that are clearly separated in affine space can cluster by a projective transformation.

The dynamic of the scaling procedure of Proposition 3.1 provides the proof of the numerical convergence properties of Toric Newton.

Theorem 6.1. The solution paths defined by a homotopy in generic position compactified torically remain separated, except possibly at the end when converging to multiple solutions.

Proof. At the start of the paths the solutions (x, z) are all mutually distinct tuples. We will show that they stay separated, as the continuation parameter t progresses. By Proposition 3.1, the components of x remain on the unit circle and the components of z may possibly go from 1 to 0.

Either a path stays inside (ℂ^{∗})^{n} or diverges, as t → 1. Only in the latter case, according to Theorem 3.1, the solution paths have the same direction. So if the path stays inside (ℂ^{∗})^{n}, its direction v in Proposition 3.1 is different from the other paths yielding distinct values for z. Thus all paths that stay inside (ℂ^{∗})^{n} remain separated, unless they converge to a multiple solution.

Still according to Theorem 3.1, only paths that diverge to solutions of the same ∂_{v}F (x) = 0 can have the same direction v. However, this direction becomes only apparent at the end of the paths, when the monomials of ∂_{v}F (x) = 0 are dominating the rest of the system. In that case, even when two paths would converge to the same z-vector, they would converge to the same x-vector only if that vector would be a multiple solution of ∂_{v}F (x) = 0. □

The following is an immediate reformulation of Theorem 6.1.

Corollary 6.1. For a homotopy in generic position, Toric Newton method converges quadratically in correcting the solutions along the paths, except when approaching singularities at the end. Regular solutions are approximated by a quadratically converging Toric Newton method.

We close this section with an important numerical observation. One may object that as z-variables can go to zero we could have to rely on very tiny values to distinguish the solutions. However, we must not look at the actual values of the z-variables, but at their exponents. These exponents have been obtained in a numerically stable way, by adding consecutive coefficients in the decomposition of the directions of the paths, see the scaling procedure of Proposition 3.1. To deal with extremal values of solutions, we usually work with a logarithmic scale, as the magnitude of the solutions is most relevant.

7. Computational Experiences

Here we report on the geometric computations showing the practical feasibility of the method.

7.1. Implementation Issues. The experiments are conducted with the polynomial homotopy continuation environment PHC [23] created, developed and maintained by the author. On the occasion of reported calculations with SAGBI homotopies [24], a new version of the package was recently released.

The construction of the triangulation is performed with the dynamic lifting algorithm [25]. This algorithm is efficient for sparse inputs. To overcome the space complexity, one should refine the fan in a lazy way subdividing only those parts of the cones that are strictly needed while path following. The implementation of Toric Newton produces values (θ, z) ∈ [0, 2π)^{n} × [0, 1]^{N}, where x_{i} = e^{iθ_{i}}. To minimize numerical round off, the affine Newton’s method must be executed in [0, 2π)^{n} × [0, 1]^{N}. This stage of implementation in PHC is still under construction.

Multi-precision arithmetic (as available in the latest release of PHC) is recommended for transforming from toric to affine coordinates at the end of the paths and for evaluation of the residuals. Only and only at that stage, we fine multi-precision arithmetic useful. Alternatively and numerically, one could treat the z_{i}’s one after the other and apply Newton’s method intermediately to correct the rounding errors. For well-conditioned affine roots, this method works well.

7.2. Geometric Samples. Here we give an impression on the complexity of the geometry involved on a 15-point configuration randomly generated for several dimensions. Table 1 summarizes the characteristics and timings. The computations were done on a 166 MHz Pentium II processor with 64 Mb internal memory running Linux.

15 points in {0, 1, 2, 3}^{n} for n = 3, 4, 5, 6, 7  cpu min, sec, millisec
n #△  vol #V  max ω deg(ϕ_{A}, z_{i}) #σ  #Σ_{P}  construct decompose
3 16  64 17  57  15 6  63  20 ms  380 ms
4 21  172 25  1,227  36 15  341  90 ms  500 ms
5 57 1,081 59  5,904  70 27 17,299  5 s 40 ms  850 ms
6 110 4,159 140 3,619,556  660 78 55,725  20 s 30 ms 1 s 450 ms
7 157 29,041 248 12,567,604  1,364 135 664,135 5 m 31 s 470 ms 3 s 30 ms

Table 1. Characteristics and timings for random samples. The first part lists the number of cells #△, the volume of the polytope and the number of facet normals #V. In the second part we compare the maximal height (max ω) to construct △ to the highest degree of a z-variable in the homogenization. The sixth column displays the maximum number of generating facet normals #σ of a cone and the total number of cones in the refinement of the fan Σ_{P}. The last part shows timings for constructing this refinement and for decomposing 1000 random vectors.

From Table 1 we observe the growth in complexity as the dimension increases. The maximal lifting value ω is generated by the dynamic lifting algorithm [25]. Although the integer lifting values generated by the implementation of the dynamic lifting algorithm can be improved by balanced floating-point liftings developed in [12], it is noteworthy that the degrees of the z-variables are of a much lesser magnitude.

Newton polytopes are in general not simple, as experienced from the #σ. The exhaustive refinement of the normal fan is quite memory exhaustive. For larger dimensions we must use a lazy decomposition of normal cones. The fast timings in constructing and consulting this geometric data structure are promising and encouraging.

7.3. A Numerical Example. The following system is inspired on an example in [16].

\{ \begin{matrix}x_{1} − 10 = 0 \\ x_{2} − x_{1}^{2} = 0 \\ x_{3} − x_{2}^{2} = 0 \\ x_{4} − x_{3}^{2} = 0 \end{matrix} ⇔ \{ \begin{matrix}γ_{11}(x_{1} − 10) + γ_{12}(x_{2} − x_{1}^{2}) + γ_{13}(x_{3} − x_{2}^{2}) + γ_{14}(x_{4} − x_{3}^{2}) = 0 \\ γ_{21}(x_{1} − 10) + γ_{22}(x_{2} − x_{1}^{2}) + γ_{23}(x_{3} − x_{2}^{2}) + γ_{24}(x_{4} − x_{3}^{2}) = 0 \\ γ_{31}(x_{1} − 10) + γ_{32}(x_{2} − x_{1}^{2}) + γ_{33}(x_{3} − x_{2}^{2}) + γ_{34}(x_{4} − x_{3}^{2}) = 0 \\ γ_{41}(x_{1} − 10) + γ_{42}(x_{2} − x_{1}^{2}) + γ_{43}(x_{3} − x_{2}^{2}) + γ_{44}(x_{4} − x_{3}^{2}) = 0 \end{matrix}

where the γ-coefficients are randomly generated complex numbers. Both systems are equivalent. The embedding at the right is only meant to illustrate how four nice quadratic polynomial equations can hide one root with large components. This problem is not recognized by conventional numerical scaling methods and is the kind of target application for Toric Newton.

8. Closing Remarks

To solve a wide class of polynomial systems in a reliable and efficient manner by homotopy continuation, a variety of different methods has been developed. Coefficient and equation scaling, often used in combination with centering, is a remedy to deal with systems that have extremal coefficients. Projective transformations compactify our working space. Polyhedral end games used in affine path following methods certify diverging paths. Toric Newton generates an optimal affine chart to approximate the magnitude of extremal roots.

Much is left undone. At a practical level, path trackers have to be developed and implemented that recognize the toric situation. Instead of complex condition numbers there is a need for structured condition numbers that measure the accuracy of the computations.

Acknowledgments. This work could not have been started without Birk Huber with whom I first learned in [15] to compute mathematical certificates of divergence. The practical experiences on path following in projective space shared to me by Tangan Gao were very inspiring. I thank T.Y. Li for a hearty conversation that led to the inclusion of assumption that homotopies need to be in generic position. I am grateful for the conversations with Alicia Dickenstein and Eduardo Cattani.

References

[1] D.N. Bernshteǐn. The number of roots of a system of equations. Functional Anal. Appl., 9(3):183–185, 1975. Translated from Funktsional. Anal. i Prilozhen., 9(3):1–4,1975.

[2] L. Blum, F. Cucker, M. Shub, and S. Smale. Complexity and Real Computation. Springer–Verlag, New York, 1997.

[3] E. Cattani and A. Dickenstein. A global view of residues in the torus. Journal of Pure and Applied Algebra, 117–118:119–144, 1997. Proceedings of MEGA’96.

[4] E. Cattani, A. Dickenstein, and B. Sturmfels. Computing multidimensional residues. In L. Gonzalez-Vega and T. Recio, editors, Algorithms in Algebraic Geometry and Applications, volume 143 of Progress in Mathematics, pages 135–164. Birkhäuser, Basel, 1996. Proceedings of MEGA’94.

[5] D. Cox, J. Little, and D. O’Shea. Using Algebraic Geometry, volume 185 of Graduate Texts in Mathematics. Springer–Verlag, New York, 1998.

[6] D.A. Cox. The homogeneous coordinate ring of a toric variety. J. Algebraic Geometry, 4(1):17–50, 1995.

[7] V. Danilov. The geometry of toric varieties. Russ. Math. Surveys, 33:97–154, 1978.

[8] J.P. Dedieu. Condition number analysis for sparse polynomial systems. In F. Cucker and M. Shub, editors, Foundations of Computational Mathematics. Selected Papers of a Conference, Held at IMPA in Rio de Janeiro, January 1997, pages 75–101. Springer–Verlag, 1997.

[9] J.P. Dedieu and M. Shub. Multihomogeneous Newton methods. IBM Technical Report, IBM Research Division, 1997.

[10] J.P. Dedieu and M. Shub. Newton and predictor-corrector methods for over-determined systems of equations. IBM Technical Report, IBM Research Division, 1998.

[11] W. Fulton. Introduction to Toric Varieties, volume 131 of Annals of Mathematics Studies. Princeton University Press, Princeton, New Jersey, 1993.

[12] T. Gao, T.Y. Li, J Verschelde, and M. Wu. Balancing the lifting values to improve the numerical stability of polyhedral homotopy continuation methods. Manuscript, 1997.

[13] I.M. Gel’fand, M.M. Kapranov, and A.V. Zelevinsky. Discriminants, Resultants and Multidimensional Determinants. Birkhäuser, Basel, 1994.

[14] B. Huber and B. Sturmfels. A polyhedral method for solving sparse polynomial systems. Math. Comp., 64(212):1541–1555, 1995.

[15] B. Huber and J. Verschelde. Polyhedral end games for polynomial continuation. Numerical Algorithms, 18(1):91–108, 1998.

[16] T. Krick and L.M. Pardo. A computational method for diophantine approximation. In L. Gonzalez-Vega and T. Recio, editors, Algorithms in Algebraic Geometry and Applications, volume 143 of Progress in Mathematics, pages 193–253. Birkhäuser, Basel, 1996. Proceedings of MEGA’94.

[17] T.Y. Li. Numerical solution of multivariate polynomial systems by homotopy continuation methods. Acta Numerica, 6:399–436, 1997.

[18] K. Meintjes and A.P. Morgan. A methodology for solving chemical equilibrium systems. Appl. Math. Comput., 22(4):333–361, 1987.

[19] A. Morgan. Solving polynomial systems using continuation for engineering and scientific problems. Prentice-Hall, Englewood Cliffs, N.J., 1987.

[20] T. Oda. Convex Bodies and Algebraic Geometry: an Introduction to the Theory of Toric Varieties. Springer–Verlag, Berlin Heidelberg New York, 1988.

[21] J.M. Rojas. Solving degenerate sparse polynomial systems faster. Manuscript available at http://www.cityu.edu.hk/ma/staff/rojas/.

[22] M.F. Roy and A. Szpirglas. Bezoutiens et résidus, affines, projectifs et dans le tore. Manuscript, 1998.

[23] J. Verschelde. PHCPACK: A general-purpose solver for polynomial systems by homotopy continuation. Paper and software available at the author’s web-pages.

[24] J. Verschelde. Numerical evidence for a conjecture in real algebraic geometry. Preprint #1998-064, MSRI, 1998. Paper and software available at the author’s web-pages.

[25] J. Verschelde, K. Gatermann, and R. Cools. Mixed-volume computation by dynamic lifting applied to polynomial system solving. Discrete Comput. Geom., 16(1):69–112, 1996.

[26] J. Verschelde, P. Verlinden, and R. Cools. Homotopies exploiting Newton polytopes for solving sparse polynomial systems. SIAM J. Numer. Anal., 31(3):915–930, 1994.

[27] L.T. Watson, S.C. Billups, and A.P. Morgan. Algorithm 652: HOMPACK: a suite of codes for globally convergent homotopy algorithms. ACM Trans. Math. Softw., 13(3):281–310, 1987.

[28] L.T. Watson, M. Sosonkina, R.C. Melville, A.P. Morgan, and H.F. Walker. HOMPACK90: A suite of Fortran 90 codes for globally convergent homotopy algorithms. ACM Trans. Math. Softw., 23(4):514– 549, 1997. Available at http://www.cs.vt.edu/˜ ltw/.

[29] G.M. Ziegler. Lectures on Polytopes, volume 152 of Graduate Texts in Mathematics. Springer–Verlag, New York, 1995.

Mathematical Sciences Research Institute, 1000 Centennial Drive, Berkeley, CA 94720-5070, USA.

Department of Mathematics, Michigan State University, East Lansing, MI 48824-1027, USA.

  E-mail address: jan@msri.org, jan@math.msu.edu or jan.verschelde@na-net.ornl.gov, and urls:
http://www.msri.org/people/members/jan http://www.mth.msu.edu/~jan