\documentstyle[12pt,amssymb]{amsart} % \headheight=8pt \topmargin=0pt % \textheight=624pt \textwidth=432pt % \oddsidemargin=18pt \evensidemargin=18pt % other margins : \headheight=8pt \topmargin=-15pt \textheight=634pt \textwidth=470pt \oddsidemargin=2pt \evensidemargin=2pt % using PicTeX : \input{prepictex} \input{pictex} \input{postpictex} % some environments : \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}{Lemma}[section] \newtheorem{corollary}{Corollary}[section] \newtheorem{definition}{Definition}[section] \newtheorem{example}{Example}[section] \newtheorem{proposition}{Proposition}[section] \newtheorem{algorithm}{Algorithm}[section] % using AMS mathematical symbols \newcommand{\nn}{{\Bbb N}} \newcommand{\ee}{{\Bbb E}} \newcommand{\zz}{{\Bbb Z}} \newcommand{\rr}{{\Bbb R}} \newcommand{\cc}{{\Bbb C}} \newcommand{\pp}{{\Bbb P}} \def\qed{\hfill $\Box$} % end of proof \def\qex{\hfill \quad\vrule height1.2ex width0.5em depth 0pt} % end of example \begin{document} \title{Toric Newton Method for Polynomial Homotopies} \author{Jan Verschelde} \address{\hskip-\parindent Mathematical Sciences Research Institute \\ 1000 Centennial Drive \\ Berkeley, CA 94720-5070, USA.} \address{\hskip-\parindent Department of Mathematics \\ Michigan State University \\ East Lansing, MI 48824-1027, USA.} \email{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} \thanks{Research at MSRI is supported by a post-doctoral fellowship. Supported in part by NSF under Grant DMS - 9804846 at Michigan State University} \bibliographystyle{plain} \begin{abstract} \noindent 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. \medskip \noindent {\bf AMS Subject Classification.} 14M25, 52B55, 65H10, 68Q40. % 14M25 Toric varieties % 52B55 Computational aspects related to convexity % 65H10 Systems of nonlinear equations % 68Q40 Symbolic Computation, algebraic computation \medskip \noindent {\bf Keywords}. polynomial system, toric varieties, Newton's method, homotopy continuation. \end{abstract} \maketitle \section{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~\cite{hv98} 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~\cite{hs95} and~\cite{vvc94} polyhedral homotopies were introduced to exploit sparsity. The number of solution paths equals the bound of Bernshte\v{\i}n's theorem. First we give some notations to state Bernshte\v{\i}n's theorems. \smallskip Denote ${\bf x}^{\bf a} = x_1^{a_1} x_2^{a_2} \cdots x_n^{a_n}$. The {\em support} of a polynomial $f({\bf x}) = \sum_{{\bf a} \in A} c_{\bf a} {\bf x}^{\bf a}$ is the set $A = \{ \ {\bf a} \in \zz^n \ | \ c_{\bf a} \not= 0 \ \}$. The {\em Newton polytope} $P$ of $f$ is the convex hull of its support~$A$. The structure of a polynomial system $F({\bf x}) = {\bf 0}$ is represented by the Newton polytopes ${\cal P} = (P_1,P_2,\ldots,P_n)$ spanned by the supports ${\cal A} = (A_1,A_2,\ldots,A_n)$. Denote $\cc^* = \cc \setminus \{ 0 \}$. In~\cite{ber75} Bernshte\v{\i}n proved that the number of isolated solutions in~$(\cc^*)^n$ of a polynomial system with generic coefficients equals the mixed volume $V_n({\cal P})$ of its Newton polytopes~$\cal P$. \smallskip In the nongeneric case we examine faces of systems. A {\em face}~$\partial_{\bf v} P$ of a polytope~$P$ is characterized by a direction~${\bf v} \in \rr^n \setminus \{ {\bf 0} \}$ and spanned by $\partial_{\bf v} A = \{ \ {\bf a} \in A \ | \ \forall {\bf b} \in A: \langle {\bf a} , {\bf v} \rangle \leq \langle {\bf b} , {\bf v} \rangle \ \}$. A {\em face polynomial of~$f$} is~$\partial_{\bf v} f({\bf x}) = \sum_{{\bf a} \in \partial_{\bf v} A} c_{\bf a} {\bf x}^{\bf a}$. A {\em face system of~$F$} is $\partial_{\bf v} F({\bf x}) = {\bf 0}$ with equations $\partial_{\bf v} f_i({\bf x}) = 0$, $i=1,2,\ldots,n$. For~$V_n({\cal P}) > 0$, Bernshte\v{\i}n's second theorem~\cite{ber75} states that $F({\bf x}) = {\bf 0}$ has fewer than~$V_n({\cal P})$ isolated solutions in~$(\cc^*)^n$ if and only if there is some $\partial_{\bf v} F({\bf x}) = {\bf 0}$ that has a solution in~$(\cc^*)^n$. \smallskip A toric variety $X$ is a solution set of a polynomial system that contains the torus $(\cc^*)^n$ as a dense open subset with $(\cc^*)^N$ as a group acting on~$X$. In section~2 we make this definition from~\cite{ful93} concrete in applying the homogenization procedure described in~\cite{cls98}. The group action is used in the third section to scale the solution components to the complex unit circle, narrowing the torus $(\cc^*)^n$ to a numerically favorable region. See~\cite{cox95}, \cite{dan78}, \cite{ful93}, \cite{gkz94}, and~\cite{oda88} for more on toric varieties. \smallskip Representing diverging and badly-scaled solutions is a common struggle in the daily life of a numerical analysist. In~\cite[Chapter 5]{mor87}, \cite{mm87}, \cite{wbm87,wsmmw97} we find coefficient and equation scaling, combined with centering, an effective and useful tool. The classical homogenization~\cite{bcss97}, \cite{li97}, \cite{mor87} 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. \smallskip 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$. \smallskip Toric Newton augments the robustness of polyhedral end games~\cite{hv98} as used in numerical homotopies~\cite{li97}. Convergence properties of Toric Newton are addressed in the fifth section. The complexity of Newton in the dense case has been analyzed in~\cite{bcss97} and has recently been extended to sparse complexity models~\cite{ded97}, multi-homogeneous~\cite{ds97} and over-determined systems~\cite{ds98}. The Toric Newton method parallels similar developments in providing algorithms using torus embeddings to compute sparse resultants~\cite{roj98}, residues~\cite{cd97},~\cite{cds96}, and Bezoutians~\cite{rs98}. All these methods rely on polytopes~\cite{zie95}. Here we give practical evidence for the performance of the polyhedral constructions. \section{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 $\lambda$ are given by the components of the direction that is perpendicular to the hyperplane, as expressed formally in the following proposition. \begin{proposition} {\rm (\cite[page~193]{gkz94})} Consider ${\displaystyle f({\bf x}) = \sum_{{\bf a} \in A} c_{\bf a} {\bf x}^{\bf a}}$. \newline If $\exists {\bf v} \in \zz^n \setminus \{ {\bf 0} \}$, $\exists d \in \zz$: $\forall {\bf a} \in A, \langle {\bf a} , {\bf v} \rangle = d$, then $f( \lambda^{v_1} x_1 , \lambda^{v_2} x_2 , \ldots , \lambda^{v_n} x_n ) = \lambda^d f({\bf x})$. \end{proposition} {\em Proof.} Denote $\lambda^{\bf v} {\bf x} = \lambda^{v_1} x_1 \lambda^{v_2} x_2 \cdots \lambda^{v_n} x_n$, then $(\lambda^{\bf v} {\bf x})^{\bf a} = (\lambda^{v_1} x_1)^{a_1} (\lambda^{v_2} x_2)^{a_2} \cdots ( \lambda^{v_n} x_n)^{a_n} = \lambda^{\langle {\bf a} , {\bf v} \rangle}$. So, ${\displaystyle f(\lambda^{\bf v} {\bf x}) = \sum_{{\bf a} \in A} c_{\bf a} (\lambda^{\bf v} {\bf x})^{\bf a} = \sum_{{\bf a} \in A} c_{\bf a} \lambda^{\langle {\bf a}, {\bf v} \rangle} {\bf x}^{\bf a} = \lambda^d f({\bf x})}$.~\qed \medskip 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. \begin{example} \label{exproj} (projective space) {\rm Let $(x,y,z) \in \pp^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 $\langle {\bf a} , {\bf v} \rangle = d$ with ${\bf v} = (1,1,1)$, $d = 1$. The homogeneous coordinates are $(\lambda^1 x , \lambda^1 y , \lambda^1 z )$.~\qex } \end{example} \begin{example} \label{exmulti} (multi-projective space) {\rm Let $(x_0,x_1,y_0,y_1) \in \pp^1 \times \pp^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 $\langle {\bf a} , {\bf v}^i \rangle = d$ with ${\bf v}^1 = (1,1,0,0)$, ${\bf v}^2 = (0,0,1,1)$, ${\bf v}^3 = (1,0,0,1)$, and $d = 1$. First we have the two usual bi-homogeneous coordinates $(\lambda^1 x_0 , \lambda^1 x_1 , y_0,y_1)$ and $(x_0,x_1,\lambda^1 y_0 , \lambda^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 $(\lambda^1 x_0, x_1 , y_0 , \lambda^1 y_1)$.~\qex } \end{example} Examples~\ref{exproj} and~\ref{exmulti} are particular cases of homogeneous coordinates. The general definition is given below. In applying the homogenization method in~\cite{cls98} we copy the original monomials~${\bf x}^{\bf 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. \begin{definition} \label{deftorhom} {\rm Consider ${\displaystyle f({\bf x}) = \sum_{{\bf a} \in A} c_{\bf a} {\bf x}^{\bf a}}$, with $P = {\rm conv}(A)$. Assume ${\rm vol}(P) > 0$. Let $P$ have $N$ facets given by $N$ inward pointing normals ${\bf v}^i \in \zz^n \setminus \{ {\bf 0} \}$ with all components relatively prime. For every ${\bf v}^i$, define $d_i = \min_{{\bf a} \in A} \langle {\bf a} , {\bf v}^i \rangle$. The facets are spanned by $\partial_{{\bf v}^i} A = \{ {\bf a} \in A | \langle {\bf a} , {\bf v}^i \rangle = d_i \}$. The {\em homogenization map $\phi_A$} introduces for every ${\bf v}^i$ a new variable $z_i$, as follows: \begin{equation} \label{eqdeftorhom} \phi_A ( {\bf x}^{\bf a} ) = {\bf x}^{\bf a} \prod_{i=1}^N z_i^{\langle {\bf a} , {\bf v}^i \rangle - d_i} \quad {\rm and} \quad \phi_A \left( \sum_{{\bf a} \in A} c_{\bf a} {\bf x}^{\bf a} \right) = \sum_{{\bf a} \in A} c_{\bf a} \phi_A ( {\bf x}^{\bf a} ). \end{equation} If ${\rm vol}(P) = 0$, then $A$ lies in some hyperplanes $\langle {\bf x} , {\bf v}^i \rangle = c_i$, and we add $\pm {\bf v}^i$ to the normals of the $k$-faces of $P$, where $k = {\rm dim}(A) - 1$. } \end{definition} The result of the homogenization is the map $\phi_A(f) : (\cc^*)^n \times \cc^N \rightarrow \cc : ({\bf x}, {\bf z}) \mapsto \phi_A(f)({\bf x}, {\bf z})$. We may allow zero values for the $z$-variables, because of the following property. \begin{proposition} \label{prophomcoord1} Let $\deg(f,x)$ denote the maximal exponent of the variable $x$ in the polynomial $f$. Then, \begin{equation} \deg\left(\phi_A \left( \sum_{{\bf a} \in A} c_{\bf a} {\bf x}^{\bf a} \right) , z_i \right) \geq 0 \quad {\rm and} \quad \deg(\phi_A ( {\bf x}^{\bf a} ) , z_i ) \left\{ \begin{array}{lll} = & 0 & {\rm if} \ {\bf a} \in \partial_{{\bf v}^i} A \\ > & 0 & {\rm if} \ {\bf a} \notin \partial_{{\bf v}^i} A \end{array} \right. . \end{equation} \end{proposition} {\em Proof.} The directions perpendicular to the facets are given by inward pointing normals, i.e.: for a facet normal ${\bf v}^i$: $d_i \leq \langle {\bf a} , {\bf v}^i \rangle$, $\forall {\bf a} \in A$. Therefore all exponents of $z_i$ of~$\phi_A$ are nonnegative.~\qed \medskip This homogenization corresponds geometrically to assigning as exponents of $z_i$ the lattice distances relative to the $i$th 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~\cite{zie95} of the polytope~$P$. The normal fan $\Sigma_P$ collects all cones that are normal to the faces of~$P$. The generators of those normal cones are the normals to the facets. \begin{example} \label{runningex} (running example) {\rm 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 ${\bf v}^1 = (1,1)$, ${\bf v}^2 = (0,1)$, ${\bf v}^3 = (-2,-1)$, and ${\bf v}^4 = (1,-1)$. In Figure~\ref{figlatt}, the homogenizing along~${\bf v}^4$ is illustrated. \begin{figure}[hbt] \begin{center} \begin{picture}(360,120)(0,0) \put(0,0){ \begin{picture}(120,120)(0,0) \put(30,20){\vector(0,1){80}} \put(20,30){\vector(1,0){120}} \put(30,60){\circle*{3}} \put(60,30){\circle*{3}} \put(90,30){\circle*{3}} \put(60,90){\circle*{3}} \thicklines \put(30,60){\line(1,-1){30}} \put(32,77){$z_4$} \put(30,60){\line(1,1){30}} \put(37,37){$z_1$} \put(60,30){\line(1,0){30}} \put(68,20){$z_2$} \put(60,90){\line(1,-2){30}} \put(78,60){$z_3$} \put(45,75){\vector(1,-1){15}} \setplotarea x from 0 to 60, y from 0 to 60 \setdashes <4pt> \setlinear \plot -20 10 80 110 / \plot 10 10 100 100 / \plot 40 10 120 90 / \plot 70 10 140 80 / \put(-23,4){${}_0$} \put(7,4){${}_1$} \put(37,4){${}_2$} \put(67,4){${}_3$} \put(15,60){$P$} \end{picture} } \put(200,0){ \begin{picture}(120,120)(0,0) \put(60,60){\circle*{3}} \put(60,60){\vector(1,1){40}} \put(102,102){${\bf v}^1$} \put(60,60){\vector(0,1){50}} \put(62,112){${\bf v}^2$} \put(60,60){\vector(-2,-1){50}} \put(0,35){${\bf v}^3$} \put(60,60){\vector(1,-1){50}} \put(115,8){${\bf v}^4$} \setplotarea x from 0 to 60, y from 0 to 60 \setdashes <4pt> \setlinear \plot 60 68 65 65 / % cone between v1 and v2 \plot 60 73 70 70 / \plot 60 80 75 75 / \plot 60 87 80 80 / \plot 60 94 85 85 / \plot 65 65 68 52 / % cone between v1 and v4 \plot 70 70 76 44 / \plot 75 75 84 36 / \plot 80 80 92 28 / \plot 85 85 100 20 / \plot 52 56 60 68 / % cone between v3 and v2 \plot 44 52 60 73 / \plot 36 48 60 80 / \plot 28 44 60 87 / \plot 20 40 60 94 / \plot 52 56 68 52 / % cone between v3 and v4 \plot 44 52 76 44 / \plot 36 48 84 36 / \plot 28 44 92 28 / \plot 20 40 100 20 / \put(105,65){$\Sigma_P$} \end{picture} } \end{picture} \caption{On the left we see that the exponents of $z_4$ are lattice distances relative to the fourth facet of~$P$ with normal~${\bf v}^4$. On the right the normal fan $\Sigma_P$ of $P$ is displayed.} \label{figlatt} \end{center} \end{figure} The homogenization is $\phi_A(f)({\bf x},{\bf 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$.~\qex } \end{example} The homogeneous coordinates correspond to the inward pointing normals to the facets. \begin{proposition} \label{prophomcoord2} Consider the facet normal ${\bf v}^i$ with corresponding $d_i$. \newline Denote $\lambda^{{\bf v}^i} {\bf x} {\bf z} = (\lambda^{v_1^i} x_1 , \lambda^{v_2^i} x_2 , \ldots , \lambda^{v_n^i} x_n , z_1,\ldots,z_{i-1}, \lambda^{-1} z_i, z_{i+1}, \ldots, z_n)$. \newline Then $\phi_A(f)(\lambda^{{\bf v}^i} {\bf x} {\bf z}) = \lambda^{d_i} \phi_A(f)({\bf x},{\bf z})$. \end{proposition} {\em Proof.} First consider one monomial~${\bf x}^{\bf a}$: \begin{eqnarray} \phi_A({\bf x}^{\bf a}) \left. \right|_{({\bf x},{\bf z}) = \lambda^{{\bf v}^i} {\bf x} {\bf z}} & = & {\bf x}^{\bf a} \prod_{j=1}^N z_j^{\langle {\bf a} , {\bf v}^j \rangle - d_j} \left. \right|_{({\bf x},{\bf z}) = \lambda^{{\bf v}^i} {\bf x} {\bf z}} \\ & = & (\lambda^{v_1^i} x_1 \lambda^{v_2^i} x_2 \cdots \lambda^{v_n^i} x_n)^{\bf a} \left( \prod_{j=1,\not=i}^N z_j^{\langle {\bf a} , {\bf v}^j \rangle - d_j} \right) (\lambda^{-1} z_i)^{\langle {\bf a} , {\bf v}^i \rangle - d_i} \\ & = & \lambda^{ \langle {\bf a} , {\bf v}^i \rangle } {\bf x}^{\bf a} \lambda^{ - \langle {\bf a} , {\bf v}^i \rangle + d_i } \prod_{j=1}^N z_j^{\langle {\bf a} , {\bf v}^j \rangle - d_j} \\ & = & \lambda^{d_i} \phi_A({\bf x}^{\bf a}) \end{eqnarray} The multilinearity of $\phi_A$ completes the proof.~\qed \medskip The transition from toric to affine space is immediate by taking $\lambda = z_i$ in every $\lambda^{{\bf v}^i} {\bf x} {\bf z}$. \begin{example} (Example~{\rm \ref{runningex}} continued.) {\rm The homogeneous coordinates corresponding to the facets are given by \begin{equation} \label{eqhomcoord} \begin{tabular}{ccccccccccccc} ( & $\lambda^1 x_1$ & , & $\lambda^1 x_2$ & , & $\lambda^{-1} z_1$ & , & $z_2$ & , & $z_3$ & , & $z_4$ & ) \\ ( & $\lambda^0 x_1$ & , & $\lambda^1 x_2$ & , & $z_1$ & , & $\lambda^{-1} z_2$ & , & $z_3$ & , & $z_4$ & ) \\ ( & $\lambda^{-2} x_1$ & , & $\lambda^{-1} x_2$ & , & $z_1$ & , & $z_2$ & , & $\lambda^{-1}z_3$ & , & $z_4$ & ) \\ ( & $\lambda^1 x_1$ & , & $\lambda^{-1} x_2$ & , & $z_1$ & , & $z_2$ & , & $z_3$ & , & $\lambda^{-1}z_4$ & ) \\ \end{tabular} \end{equation} for the respective normals ${\bf v}^i$, $i=1,2,3,4$. Observe that as $\lambda \rightarrow 0$, the $x$-values in affine space go either to zero or infinity. We have a multitude of homogeneous coordinates by taking combinations, see $\Sigma_P$ in Figure~\ref{figlatt}. Positive combinations of any two vectors normal to adjacent facets correspond to constructing a vector inside the cone $\sigma \in \Sigma_P$ that has these vectors among its generators. The relations~(\ref{eqhomcoord}) illustrate how the group~$(\cc^*)^N$ acts on the variety.~\qex } \end{example} For general polynomial systems not all supports are equal. Suppose there are $r$ different supports which we group into ${\cal A} = (A_1,A_2,\ldots,A_r)$ spanning the Newton polytopes ${\cal P} = (P_1,P_2,\ldots,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~$\cal 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~\ref{deftorhom} is extended to systems with different supports as follows, taking $r=n$ for notational convenience. \begin{definition} \label{defmixtorhom} {\rm Let $\cal A$ be the supports of $F({\bf x}) = {\bf 0}$. Suppose the polytope spanned by $A_i$ has $N_i$ facet normals. Abbreviate ${\bf z} = ({\bf z}^1, {\bf z}^2, \ldots, {\bf z}^n)$. The {\em mixed homogenization map $\phi_{\cal A}$} is \begin{equation} \phi_{\cal A}(F) : (\cc^*)^n \times \cc^{N_1} \times \cc^{N_2} \times \cdots \times \cc^{N_n} \rightarrow \cc^n : ({\bf x}, {\bf z}) \mapsto \phi_{\cal A}(F)({\bf x}, {\bf z}), \end{equation} with $\phi_{\cal A}(F)({\bf x}, {\bf z}) = (\phi_{A_1}(f_1)({\bf x},{\bf z}^1), \phi_{A_2}(f_2)({\bf x},{\bf z}^2)$, $\ldots$, $\phi_{A_n}(f_n)({\bf x},{\bf z}^n))$ and $\phi_{A_i}(f_i)$ as in~{\rm (\ref{eqdeftorhom})}. For a homotopy $H({\bf x},t) = {\bf 0}$, the continuation parameter $t$ is treated as a constant in setting up the homogenized homotopy~$\phi_{\cal A}(H)({\bf x},{\bf z},t) = {\bf 0}$. } \end{definition} For several polytopes sharing the same facet normal, we can use the same corresponding $z$-variable. With this economical use of $z$-variables, Definition~\ref{defmixtorhom} yields the same homogenization map $\phi_A$ as in Definition~\ref{deftorhom} when applied to ${\cal A} = (A,A,\ldots,A)$. Since the general homogenization is also determined by inward pointing normals ${\bf v}^i$, Propositions~\ref{prophomcoord1} and~\ref{prophomcoord2} have a direct generalization to the case of different supports and are omitted. \begin{proposition} \label{propfacesols1} Suppose $A$ is the support of $f({\bf x})$. \newline Let $({\bf x}^*,{\bf z}^*) \in (\cc^*)^n \times \cc^N$ $:$ $\phi_A(f)({\bf x}^*,{\bf z}^*) = 0$. \begin{enumerate} \item If $z^*_i = 0$ and $z^*_j \not= 0$, $\forall j \not= i$, then $\partial_{{\bf v}^i}f({\bf x}) = 0$ has a solution in~$(\cc^*)^n$. \item Let $I = \{ \ i \ | \ z^*_i = 0, \# (\cap \partial_{{\bf v}^i} A) \geq 2 \ \}$, $z^*_j \not= 0$, $\forall j \not\in I$, and ${\bf v} = {\displaystyle \sum_{i \in I} {\bf v}^i}$. \newline Then $\partial_{\bf v}f({\bf x}) = 0$ has a solution in~$(\cc^*)^n$. \end{enumerate} The solutions ${\bf x}^*$ of~$\partial_{\bf v} f$ correspond to solutions of~$f$ with components equal to~$0$ and/or~$\infty$. \end{proposition} {\em Proof}. Although the second statement generalizes the first, we distinguish two cases: \begin{enumerate} \item Suppose there is only one $z^*_i = 0$ with corresponding facet normal~${\bf v}^i$. Then all monomials in $\phi_A(f)$ whose exponent vector does not lie on the $i$th facet vanish when $z_i = z^*_i = 0$. For $z^*_j \not= 0$, $\forall j \not= i$, the homogeneous coordinates $\lambda^{{\bf v}^j} {\bf x}^*{\bf z}^*$ (as in Proposition~\ref{prophomcoord2}) are used successively, each time setting $\lambda = z^*_j$, $\forall j \not= i$, to obtain a solution in~$(\cc^*)^n$ of $\partial_{{\bf v}^i} f({\bf x}) = {\bf 0}$. \item 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~$\partial_{\bf v} f$ to admit a solution in $(\cc^*)^n$. Successive application of Proposition~\ref{prophomcoord2} setting $\lambda = z^*_j$, $\forall j \notin I$, yields a solution in $(\cc^*)^n$ of $\partial_{\bf v} f({\bf x}) = {\bf 0}$. \end{enumerate} The homogeneous coordinates contain as exponents for $\lambda$ the components of the normals. Any normal $\bf v$ has at least one component different from zero, call it $v_k$. Setting $\lambda = z_i \rightarrow 0$ to transform to affine space will lead to $x_k \rightarrow \infty$, if $v_k < 0$, or $x_k \rightarrow 0$, in case $v_k > 0$.~\qed \medskip In~\cite{cls98}, solutions that have no $z$-variables equal to zero, or only in cases covered by Proposition~\ref{propfacesols1} are called~{\em nontrivial}. Indeed, without the conditions on~$I$ in Proposition~\ref{propfacesols1}, all monomials could vanish. The generalization of Proposition~\ref{propfacesols1} to polynomial systems is immediate. \begin{proposition} \label{propfacesols} Let $F({\bf x}) = {\bf 0}$ have supports~$\cal A$. The nontrivial solutions of $\phi_{\cal A}(F)({\bf x},{\bf z}) = {\bf 0}$ with some $z$-variables equal to zero correspond to solutions of face systems. \end{proposition} 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. \section{Scaling Coefficients and Solutions} The variables introduced by the homogenization allow to scale $\phi_A(f)({\bf x},{\bf z}) = {\bf 0}$, so that every component of~$\bf x$ lies on the unit circle in~$\cc^*$, while keeping the $z$-variables bounded. \begin{proposition} \label{propscaling} Suppose ${\bf x} \in (\cc^*)^n$ satisfies~$\phi_A(f)({\bf x},{\bf z}) = {\bf 0}$ for fixed $\bf z$ with $0 < z_j \leq 1$, $\forall j = 1,2,\ldots,N$. Let ${\bf v} = (-\log|x_1|,-\log|x_2|,\ldots,-\log|x_n|)$, with $\log$ of base 10. Then \begin{enumerate} \item $\exists {\bf c} \in \rr_+^n$, $\exists J = (j_1,j_2,\ldots,j_n)$, $j_i \in \{ 1,2,\ldots,N \}:$ ${\displaystyle {\bf v} = \sum_{i=1}^n c_i {\bf v}^{j_i}}$. \item The homogeneous coordinates \begin{equation} \left( \left. \lambda^{\left( - \sum_{i=1}^n c_i {\bf v}^{j_i} \right)_k} x_k \right|_{k=1}^n, \ldots, \left. \lambda^{c_i} z_{j_i} \right|_{i=1}^n , \ldots \right) \end{equation} with $\lambda = 0.1$ scale $({\bf x},{\bf z})$ to $({\widetilde {\bf x}},{\widetilde {\bf z}})$ with $|{\widetilde x}_k| = 1$ and $0 < {\widetilde z}_j \leq 1$. \end{enumerate} \end{proposition} {\em Proof}. The above statements are true because: \begin{enumerate} \item The normal fan $\Sigma_P$ of $P$ is complete. Whence every vector ${\bf v}$ can be written as a positive combination of the generators of a cone $\sigma \in \Sigma_P$. \item $\left( - \sum_{i=1}^n c_i {\bf v}^{j_i} \right)_k = \log |x_k|$, so ${\widetilde x}_k = (0.1)^{\log|x_k|} x_k = \frac{x_k}{|x_k|}$, whence $|{\widetilde x}_k| = 1$. $0 < z_{j_i} \leq 1 \Rightarrow 0 < {\widetilde z}_{j_i} = (0.1)^{c_i} z_{j_i} \leq 1$, as $c_i \geq 0$. Thus $0 < {\widetilde z}_j \leq 1$.~\qed \end{enumerate} \medskip \noindent Note that any base $b$ can be used in the logarithm, it suffices to replace 0.1 in Proposition~\ref{propscaling} by~$b^{-1}$. \begin{example} (toric scaling analogue to projective space) {\rm Suppose the Newton polytope is a simplex. The homogeneous coordinates are \begin{equation} \begin{tabular}{ccccccccccc} ( & $\lambda^1 x_1$ & , & $\lambda^0 x_2$ & , & $\lambda^{-1} z_1$ & , & $z_2$ & , & $z_3$ & ) \\ ( & $\lambda^0 x_1$ & , & $\lambda^1 x_2$ & , & $z_1$ & , & $\lambda^{-1} z_2$ & , & $z_3$ & ) \\ ( & $\lambda^{-1} x_1$ & , & $\lambda^{-1} x_2$ & , & $z_1$ & , & $z_2$ & , & $\lambda^{-1} z_3$ & ) \\ \end{tabular} \end{equation} for the respective facet normals ${\bf v}^1 = (1,0)$, ${\bf v}^2 = (0,1)$, and ${\bf v}^3 = (-1,-1)$. Since this is a {\em 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: \begin{enumerate} \item $(x_1,x_2) = (2,2)$, ${\bf v} = (-\log(2),-\log(2))$, $J=(3,3)$, ${\bf c} = (\log(2),0)$ \newline $\Rightarrow$ $z_3 = (0.1)^{\log(2)} z_3 = \frac{1}{2} z_3$. \item $(x_1,x_2) = (2,1)$, ${\bf v} = (-\log(2),-\log(1))$, $J=(2,3)$, ${\bf c} = (\log(2),\log(2))$ \newline $\Rightarrow$ $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$. \end{enumerate} In the first case it suffices to use one of the homogeneous coordinates, while in the second case a combination is constructed.~\qex } \end{example} \begin{example} (Example~{\rm \ref{runningex}} continued.) {\rm Here we also consider two situations: \begin{enumerate} \item $(x_1,x_2) = (\frac{1}{2},2)$, ${\bf v} = (-\log(\frac{1}{2}),-\log(2))$, $J=(4,4)$, ${\bf c} = (\log(2),0)$ \newline $\Rightarrow$ $z_4 = (0.1)^{\log(2)} z_4 = \frac{1}{2} z_4$. \item $(x_1,x_2) = (2,2)$, ${\bf v} = (-\log(2),-\log(1))$, $J=(2,3)$, ${\bf c} = (\log(2),\log(2))$ \newline $\Rightarrow$ $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$. \end{enumerate} The scaling relation is constructed by decomposing $\bf v$ in $\Sigma_P$, see Figure~\ref{figlatt}.~\qex} \end{example} In the general case when the system has an $r$-tuple of different supports~$\cal A$, we decompose ${\bf v} = - \log|{\bf x}|$ in every component of $\Sigma_{\cal P} = (\Sigma_{P_1},\Sigma_{P_2}, \ldots, \Sigma_{P_r})$ to compute new values of the $z$-variables which are distinct for every different support. \smallskip Propositions~\ref{propfacesols} and~\ref{propscaling} allow us to give a concise proof of Bernshte\v{\i}n's second theorem formulated below as Theorem~\ref{theoberntwo} using a homotopy argument. \begin{definition} {\rm We say that the homotopy $H({\bf x},t) = {\bf 0}$ is {\em in generic position} with respect to the supports $\cal A$ that span the Newton polytopes $\cal P$, if $\forall t \in [0,1)$: $H({\bf x},t) = {\bf 0}$ has supports~$\cal A$ and has exactly~$V_n({\cal P})$ regular solutions in~$(\cc^*)^n$. } \end{definition} Given a system~$F({\bf x}) = {\bf 0}$, the canonical example of a homotopy in generic position is \begin{equation} H({\bf x},t) = (1-t) F^{(0)}({\bf x}) + t F({\bf x}) = {\bf 0}, \quad t \in [0,1], \end{equation} where the system~$F^{(0)}({\bf x}) = {\bf 0}$ has random complex coefficients and the same supports as~$F$. \begin{theorem} \label{theoberntwo} Consider~$F({\bf x}) = {\bf 0}$, with supports~$\cal A$, Newton polytopes~$\cal P$ and~$V_n({\cal P}) \not= 0$. \begin{displaymath} \# \{ \ m_{\bf x} \ | \ F({\bf x}) = {\bf 0}, {\bf x}~{\rm isolated} \ \} < V_n({\cal P}) ~ \Leftrightarrow ~ \exists {\bf v} \in \rr^n \setminus \{ {\bf 0} \}, \exists {\bf x} \in (\cc^*)^n: \partial_{\bf v} F({\bf x}) = {\bf 0} \end{displaymath} where $m_{\bf x}$ denotes the multiplicity of the isolated solution $\bf x$ of $F({\bf x}) = {\bf 0}$. \end{theorem} {\em Proof}. Construct any homotopy $H({\bf x},t) = {\bf 0}$ that is in generic position with respect to~$\cal P$ and for which $H({\bf x},t=1) = F({\bf x})$. Consider the homogenized homotopy $\phi_{\cal A}(H)({\bf x},{\bf z},t) = {\bf 0}$. \begin{itemize} \item[$\Rightarrow$] If $F({\bf x}) = {\bf 0}$ has fewer than~$V_n({\cal P})$ isolated solutions in~$(\cc^*)^n$, then, as $t \rightarrow 1$, either some solution paths must diverge outside~$(\cc^*)^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~\ref{propscaling}, the only way to diverge is that some $z_i \rightarrow 0$. Then, according to Proposition~\ref{propfacesols}, we have a solution in $(\cc^*)^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 $(\cc^*)^n$ for which some $z_i$ is zero. Then we can again apply Proposition~\ref{propfacesols}. \item[$\Leftarrow$] If there is a face system that has a solution ${\bf x} \in (\cc^*)^n$, then in the corresponding tuple of $z$-variables for which $\phi_{\cal A}(H)({\bf x},{\bf z},1) = {\bf 0}$, $\exists i$: $z_i = 0$. Since the homotopy is in generic position, there is at least one solution path that ended at~$({\bf x},{\bf z})$. By Proposition~\ref{propfacesols}, the transition to affine coordinates shows that the solution path that ended at this couple~$({\bf x},{\bf z})$ diverged to a solution that has components at infinity and/or zero components. So $F({\bf x}) = {\bf 0}$ must have fewer than~$V_n({\cal P})$ isolated solutions in~$(\cc^*)^n$.~\qed \end{itemize} \section{The Static Scaling in Projective Newton} Let us recall Newton's method in projective space (see~\cite{li97}). We assume we have a system $H = (h_1,h_2,\ldots,h_n)$ of homogeneous polynomials: $h_i(\lambda {\bf z}) = \lambda^{\deg(h_i)} h_i({\bf z})$ for $i= 1,2,\ldots,n$ with $\lambda {\bf z} = (\lambda z_0, \lambda z_1, \ldots , \lambda z_n)$ the usual homogeneous coordinates. For a solution at infinity: $z_0 = 0$. \begin{algorithm} \label{algprojnewt} (Newton's method in projective space) {\rm % \newline \begin{tabular}{lcr} Input: $(H,{\bf z})$. & \ \ \ \ & {\em (Homogeneous system, approximate root)} \\ Output: $\bf z$. & & {\em refined approximation} \\ \\ ${\displaystyle i : |z_i| = \max_{j=1}^n |z_j|}$; & & {\em maximal component} \\ ${\bf z} := \left( \underbrace{\frac{z_0}{z_i}, \ldots, \frac{z_{i-1}}{z_i}} , 1 , \underbrace{\frac{z_{i+1}}{z_i} , \ldots , \frac{z_n}{z_i}} \right)$; & & {\em scale the approximation} \\ ${\bf x} := \left( ~ \overbrace{z_1,\ldots,z_{i-1}} ~ , ~ \overbrace{z_{i+1},\ldots,z_n} ~ \right)$; & & {\em switch to an affine chart} \\ $F := \left. H \right|_{z_i = 1}$; & & {\em select evaluation map} \\ ${\bf x} := {\bf x} - D_{\bf x}^{-1} F({\bf x}) F({\bf x})$; & & {\em affine Newton method} \\ ${\bf z} := (x_1,\ldots,x_i,1,x_{i+1},\ldots,x_n)$. & & {\em refined approximation} \end{tabular} } \end{algorithm} In the context of path following, $H$ is a parametrized family of systems. \begin{example} (projective Newton on a sparse system) {\rm 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~\ref{fignewt}. 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. \begin{figure}[hbt] \begin{center} \begin{picture}(420,100)(0,0) \put(0,0){ \begin{picture}(120,100)(0,0) \put(30,40){\circle*{3}} \put(60,10){\circle*{3}} \put(90,10){\circle*{3}} \put(60,70){\circle*{3}} \put(27,70){\line(1,0){6}} \thicklines \put(30,40){\line(1,-1){30}} \put(30,40){\line(1,1){30}} \put(60,10){\line(1,0){30}} \put(60,70){\line(1,-2){30}} \setplotarea x from 0 to 60, y from 0 to 60 \setdashes <4pt> \setlinear \plot 30 10 30 100 / \plot 30 10 120 10 / \plot 30 100 120 10 / \put(25,4){${}_0$} \put(58,2){${}_1$} \put(88,2){${}_2$} \put(118,2){${}_3$} \put(22,39){${}_1$} \put(22,69){${}_2$} \put(22,99){${}_3$} \put(80,60){$z_0$} \put(70,-3){$z_2$} \put(13,55){$z_1$} \put(50,35){$P$} \end{picture} } \put(200,60){ $f({\bf x}) = c_{10} x_1 + c_{01} x_2 + c_{20} x_1^2 + c_{12} x_1 x_2^2$} \put(200,30){ ${\widetilde f}({\bf 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$} \end{picture} \caption{On the left we see how the polytope $P$ of $f$ is approximated by a simplex which corresponds to the usual projective homogenization ${\widetilde f}$ given on the right.} \label{fignewt} \end{center} \end{figure} In case of divergence, $z_1$ and/or $z_2$ go to infinity. The scaling step in Algorithm~\ref{algprojnewt} selects $i=1$ or $i=2$ and $z_0$ strives to zero. In that case we see in Figure~\ref{fignewt} 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.~\qex } \end{example} \begin{example} (projective Newton leads to clustering) {\rm Consider the following simple system: \begin{equation} \left\{ \begin{array}{r} x_1 - 10^8 = 0 \\ x_2^2 - 1 = 0 \end{array} \right. \end{equation} In affine space, the solutions $(x_1,x_2) = (10^8,\pm 1)$ are well-separated. In projective space, they are represented as $(z_0,z_1,z_2) = (1,10^8,\pm 1)$ and scaled to $(z_0,z_1,z_2) = (10^{-8},1,\pm 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,\pm 1) \mapsto (z_{10},z_{11},z_{20},z_{21}) = (10^{-8},1,1,\pm 1)$, avoiding the clustering of solutions.~\qex } \end{example} The most important difference with multi-projective Newton is the dynamic use of the scaling procedure of Proposition~\ref{propscaling} which guarantees distinct paths, as formulated below in Theorem~\ref{theoseparate}. \section{Toric Newton scales in a Dynamic Way} The set up of Toric Newton enables an efficient numerical scaling while path following. \begin{algorithm} \label{algtorsymnewt} (Toric Newton -- symbolic set-up) {\rm % \newline \begin{tabular}{lcr} Input: $(A,\triangle)$. & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ & {\em (support , triangulation)} \\ Output: $H$; & & {\em homogeneous coordinates matrix} \\ ~~~~~~~~~~~ $\Sigma_P$. & & {\em simplicial normal fan} \\ \\ 1. $V$ := Facet\_Normals$(A,\triangle)$; & & {\em extract facet normals} \\ 2. $H$ := Homogenize$(A,V)$; & & {\em relative lattice distances} \\ 3. $\Sigma_P$ := Cone\_Subdivision$(H,A,V)$. & & {\em refine normal fan} \end{tabular} } \end{algorithm} 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 $\# \triangle$ in any triangulation $\triangle$ of~$A$. Recall that we are interested in {\em all} isolated solutions and to count the number of roots we compute volumes. So $\# \triangle$ is a lower bound for the number of isolated solutions. \begin{proposition} (Cost of Step 1 in Algorithm~\ref{algtorsymnewt}.) Consider a triangulation $\triangle$ of the polytope $P$ spanned by~$A$. Then all facets of~$P$ can be enumerated in $O(\# \triangle n ( n^3 + \#A ))$ operations. \end{proposition} {\em 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(\# \triangle n ( n^3 + \#A ))$.~\qed \medskip The number of arithmetic operations of Step 2 of Algorithm~\ref{algtorsymnewt} 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~\ref{figconeref}. In Step~3 we make the fan simplicial subdividing the normal cones. One method to subdivide normal cones is to cut corners (see~\cite[page~190]{gkz94} and~\cite[pages~54--55]{zie95}) at vertices that have more than $n$ facets, as illustrated in Figure~\ref{figconeref}. We refer to~\cite[Section~2.6]{ful93} for more on refinements of fans in connection with toric varieties. \smallskip For a normal cone $\sigma \in \Sigma_P$ at the vertex ${\bf a}$ we take the average~${\overline {\bf v}}$ of the facet normals that generate~$\sigma$. Determine $c$ such that the hyperplane $\langle {\bf x} , {\overline {\bf v}} \rangle = c$ separates ${\bf a}$ from the rest of the vertices that span the polytope. The normal cones $\tau$ to the intersection points of the hyperplane with the edges of the polytope refine the normal cone~$\sigma$. The recursive structure of the normal fan $\Sigma_P$ is illustrated at the right of Figure~\ref{figconeref}. In the decomposition of a vector ${\bf v}$ in a subdivision of $\sigma$, we replace the average vector ${\overline {\bf v}}$ in terms of the other generators of $\sigma$, so that we obtain a decomposition that uses only facet normals. \begin{figure}[hbt] \begin{center} \begin{picture}(400,120)(0,0) \put(-30,0){ \begin{picture}(120,120)(0,0) \put(30,60){\vector(0,1){50}} \put(18,110){$a_3$} \put(30,60){\vector(1,0){100}} \put(130,55){$a_2$} \put(30,60){\vector(1,-4){15}} \put(30,0){$a_1$} \put(30,60){\circle*{3}} \put(15,63){${}_{000}$} \put(90,60){\circle*{3}} \put(92,63){${}_{020}$} \put(40,20){\circle*{3}} \put(25,18){${}_{200}$} \put(100,20){\circle*{3}} \put(102,18){${}_{220}$} % \put(65,40){\circle*{3}} % 110 \put(65,95){\circle*{3}} \put(59,100){${}_{111}$} \put(45,85){{\scriptsize ${\bf b}^1$}} \put(47,70){{\scriptsize ${\bf b}^2$}} \put(75,85){{\scriptsize ${\bf b}^3$}} \put(65,65){{\scriptsize ${\bf b}^4$}} \thicklines \put(30,60){\line(1,-4){10}} \put(90,60){\line(1,-4){10}} \put(30,60){\line(1,0){60}} \put(40,20){\line(1,0){60}} \setplotarea x from 0 to 120, y from 0 to 120 \setlinear \plot 30 60 65 95 / \plot 90 60 65 95 / \plot 40 20 65 95 / \plot 100 20 65 95 / \put(55,85){\circle*{3}} \put(58,73){\circle*{3}} \put(75,73){\circle*{3}} \put(72,85){\circle*{3}} \plot 55 85 58 73 / \plot 55 85 72 85 / \plot 58 73 75 73 / \plot 72 85 75 73 / \end{picture} } \put(160,0){ \begin{picture}(200,120)(0,0) \put(0,90){{\scriptsize $A = \left[ \begin{array}{c} 1~1~1 \\ 0~0~0 \\ 2~0~0 \\ 0~2~0 \\ 2~2~0 \\ \end{array} \right]$}} \put(70,90){{\scriptsize $V = \left[ \begin{array}{r} 0~\!~~~~0~\!~~~~1 \\ 0~\!~~~~1~-1 \\ -1~\!~~~~0~-1 \\ 1~\!~~~~0~-1 \\ 0~-1~-1 \\ \end{array} \right]$}} \put(155,110){${\bf v}^1$} \put(155,100){${\bf v}^2$} \put(155,90){${\bf v}^3$} \put(155,80){${\bf v}^4$} \put(155,70){${\bf v}^5$} \put(190,90){{\scriptsize $H = \left[ \begin{array}{c} 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{array} \right]$}} \put(0,30){{\scriptsize $\Sigma_P = \left[ \begin{array}{c} 1~1~1 \\ 0~0~0 \\ 2~0~0 \\ 0~2~0 \\ 2~2~0 \\ \end{array} \right]$}} \put(57,50){\vector(1,0){100}} \put(57,41){\vector(1,0){20}} \put(57,32){\vector(1,0){20}} \put(57,23){\vector(1,0){20}} \put(57,14){\vector(1,0){20}} \put(80,10){$({\bf v}^1,{\bf v}^3,{\bf v}^5)$} \put(80,19){$({\bf v}^1,{\bf v}^4,{\bf v}^5)$} \put(80,28){$({\bf v}^1,{\bf v}^2,{\bf v}^3)$} \put(80,37){$({\bf v}^1,{\bf v}^2,{\bf v}^4)$} \put(160,30){{\scriptsize $\left[ \begin{array}{c} {\bf b}^1 \\ {\bf b}^2 \\ {\bf b}^3 \\ {\bf b}^4 \\ \end{array} \right]$}} \put(185,45){\vector(1,0){20}} \put(185,36){\vector(1,0){20}} \put(185,27){\vector(1,0){20}} \put(185,18){\vector(1,0){20}} \put(210,41){$({\overline {\bf v}},{\bf v}^2,{\bf v}^4)$} \put(210,32){$({\overline {\bf v}},{\bf v}^2,{\bf v}^3)$} \put(210,23){$({\overline {\bf v}},{\bf v}^4,{\bf v}^5)$} \put(210,14){$({\overline {\bf v}},{\bf v}^3,{\bf v}^5)$} \end{picture} } \end{picture} \caption{Refinement of a normal cone of a nonsimple polytope by cutting the corner at~$(1,1,1)$. The points ${\bf b}^i$ lie on the intersection of the hyperplane cut with the edges that contain $(1,1,1)$. } \label{figconeref} \end{center} \end{figure} 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~\cite[page~55]{zie95}, edges are in one-to-one correspondence with vertices of the polytope cut out by the separating hyperplane. Recursive application proves the following. \smallskip \begin{proposition} \label{propcostconeref} (Cost of Step 3 in Algorithm~\ref{algtorsymnewt}.) To triangulate an $n$-dimensional cone at most $n-2$ separating hyperplanes are needed. \end{proposition} 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 $\Sigma_P$, decomposing only those normal cones when needed in the scaling procedure. \medskip Just as projective Newton method, we switch to an appropriate affine chart of dimension~$n$. \begin{algorithm} \label{algtornumnewt} (Toric Newton -- number crunching) {\rm \newline \begin{tabular}{lcr} Input: $\phi = \phi_A(F)$; & \ \ \ \ \ \ \ \ \ \ \ & {\em system in homogeneous coordinates} \\ ~~~~~~~~~ $({\bf x},{\bf z})$; & & {\em approximate root with $|x_i| = 1$ and $|z_j| \leq 1$} \\ ~~~~~~~~~ $\Sigma_P$. & & {\em simplicial normal fan} \\ Output: $({\bf x},{\bf z})$. & & {\em refined approximation with $|x_i| = 1$ and $|z_j| \leq 1$} \\ \end{tabular} \begin{tabular}{lcr} $F := \left. \phi \right|_{Fix~{\bf z}}$; & & {\em evaluation map} \\ ${\bf x} := {\bf x} - D_{\bf x}^{-1} F({\bf x}) F({\bf x})$; & \ \ \ \ \ \ \ \ \ \ \ \ & {\em affine Newton method} \\ ${\bf v} := -\log|{\bf x}|$; & & {\em direction of correction} \\ $J = (j_1,j_2,\ldots,j_n)$ & & {\em find $\sigma \in \Sigma_P : {\bf v} \in \sigma$} \\ ~~~ := Localize\_Position$(\Sigma_P,{\bf v})$; & & {\em $\sigma = {\rm span}({\bf v}^{j_1},{\bf v}^{j_2},\ldots,{\bf v}^{j_n})$} \\ ${\displaystyle {\bf c} = (c_1,c_2,\ldots,c_n) : {\bf v} = \sum_{i=1}^n c_i {\bf v}^{j_i}}$; & & {\em decompose ${\bf v}$} \\ for all $i \in \{1,2,\ldots,n \}$ do & & \\ ~~ $x_i := \frac{x_i}{|x_i|}$; $z_{j_i} := (0.1)^{c_i} z_{j_i}$; & & {\em scale solution and coefficients} \\ end for. & & \end{tabular} } \end{algorithm} In analyzing the cost of localizing the position of the direction vector ${\bf v}$ we assume that the fan $\Sigma_P$ has a recursive simplicial structure, as illustrated in Figure~\ref{figconeref}. \begin{proposition} (Cost of {\rm Localize\_Position} in Algorithm~\ref{algtornumnewt}.) For ${\bf v} \in \rr^n$ and $\Sigma_P$, $P = {\rm conv}(A)$, finding $\sigma \in \Sigma_P$ so that ${\bf v} \in \sigma$ takes at least $\#A$ and at most $\#A (n-1)$ inner products. \end{proposition} {\em Proof.} If $P$ is simple, then $\sigma$ is the normal cone at the point ${\bf a} \in A$ for which $\langle {\bf a} , {\bf v} \rangle = \min_{{\bf b} \in A} \langle {\bf b} , {\bf v} \rangle$. Computing this point requires $\#A$ inner products. If $P$ is not simple, there are at most $n-2$ refinement levels, according to Proposition~\ref{propcostconeref}. At every level, there are at most $\#A$ points to consider.~\qed \medskip For $r$-tuples of supports~$\cal A$, we construct $\Sigma_{\cal P} = (\Sigma_{P_1},\Sigma_{P_2}, \ldots, \Sigma_{P_r})$. The cost of the geometric constructions is scaled by a factor~$r$. From our experiences with mixed-volume computations~\cite{vgc96} we find linear upscaling in the cost rather remarkable. \section{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. \smallskip Affine path following methods deal well with paths diverging outside $(\cc^*)^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~\cite{hv98} 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. \smallskip The dynamic of the scaling procedure of Proposition~\ref{propscaling} provides the proof of the numerical convergence properties of Toric Newton. \smallskip \begin{theorem} \label{theoseparate} The solution paths defined by a homotopy in generic position compactified torically remain separated, except possibly at the end when converging to multiple solutions. \end{theorem} {\em Proof}. At the start of the paths the solutions $({\bf x},{\bf z})$ are all mutually distinct tuples. We will show that they stay separated, as the continuation parameter $t$ progresses. By Proposition~\ref{propscaling}, the components of ${\bf x}$ remain on the unit circle and the components of ${\bf z}$ may possibly go from 1 to~0. \smallskip Either a path stays inside $(\cc^*)^n$ or diverges, as $t \rightarrow 1$. Only in the latter case, according to Theorem~\ref{theoberntwo}, the solution paths have the same direction. So if the path stays inside $(\cc^*)^n$, its direction ${\bf v}$ in Proposition~\ref{propscaling} is different from the other paths yielding distinct values for ${\bf z}$. Thus all paths that stay inside $(\cc^*)^n$ remain separated, unless they converge to a multiple solution. \smallskip Still according to Theorem~\ref{theoberntwo}, only paths that diverge to solutions of the same $\partial_{\bf v} F({\bf x}) = {\bf 0}$ can have the same direction~${\bf v}$. However, this direction becomes only apparent at the end of the paths, when the monomials of $\partial_{\bf v} F({\bf x}) = {\bf 0}$ are dominating the rest of the system. In that case, even when two paths would converge to the same ${\bf z}$-vector, they would converge to the same ${\bf x}$-vector only if that vector would be a multiple solution of~$\partial_{\bf v} F({\bf x}) = {\bf 0}$.~\qed \medskip The following is an immediate reformulation of Theorem~\ref{theoseparate}. \begin{corollary} 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. \end{corollary} 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~\ref{propscaling}. To deal with extremal values of solutions, we usually work with a logarithmic scale, as the magnitude of the solutions is most relevant. \section{Computational Experiences} Here we report on the geometric computations showing the practical feasibility of the method. \subsection{Implementation Issues} The experiments are conducted with the polynomial homotopy continuation environment~PHC~\cite{ver97} created, developed and maintained by the author. On the occasion of reported calculations with SAGBI homotopies~\cite{ver98}, a new version of the package was recently released. \smallskip The construction of the triangulation is performed with the dynamic lifting algorithm~\cite{vgc96}. 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 $(\theta,{\bf z}) \in [0,2 \pi)^n \times [0,1]^N$, where $x_i = e^{i \theta_i}$. To minimize numerical round off, the affine Newton's method must be executed in~$[0,2 \pi)^n \times [0,1]^N$. This stage of implementation in~PHC is still under construction. \smallskip 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. \subsection{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~\ref{tabgeosamples} summarizes the characteristics and timings. The computations were done on a 166 MHz Pentium II processor with 64 Mb internal memory running Linux. {\small \begin{table}[ht] \begin{center} \begin{tabular}{|c||r|r|r||r|r|r|r||r|r|} \hline & \multicolumn{7}{c||}{15 points in $\{ 0,1,2,3 \}^n$ for $n=3,4,5,6,7$} & \multicolumn{2}{c|}{cpu min, sec, millisec} \\ \cline{2-10} n & $\# \triangle$ & vol & $\# V$ & max $\omega$ & ${\scriptsize \deg(\phi_A,z_i)}$ & $\# \sigma$ & $\# \Sigma_P$ & construct & decompose \\ \hline \hline 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 \\ \hline \end{tabular} \caption{Characteristics and timings for random samples. The first part lists the number of cells~$\# \triangle$, the volume of the polytope and the number of facet normals~$\# V$. In the second part we compare the maximal height (max~$\omega$) to construct~$\triangle$ to the highest degree of a $z$-variable in the homogenization. The sixth column displays the maximum number of generating facet normals~$\# \sigma$ of a cone and the total number of cones in the refinement of the fan~$\Sigma_P$. The last part shows timings for constructing this refinement and for decomposing 1000 random vectors.} \label{tabgeosamples} \end{center} \end{table} } From Table~\ref{tabgeosamples} we observe the growth in complexity as the dimension increases. The maximal lifting value~$\omega$ is generated by the dynamic lifting algorithm~\cite{vgc96}. Although the integer lifting values generated by the implementation of the dynamic lifting algorithm can be improved by balanced floating-point liftings developed in~\cite{glvw97}, it is noteworthy that the degrees of the $z$-variables are of a much lesser magnitude. \smallskip Newton polytopes are in general not simple, as experienced from the~$\# \sigma$. 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. \subsection{A Numerical Example} The following system is inspired on an example in~\cite{kp96}. \begin{displaymath} \left\{ \begin{array}{c} x_1 - 10 = 0 \\ x_2 - x_1^2 = 0 \\ x_3 - x_2^2 = 0 \\ x_4 - x_3^2 = 0 \\ \end{array} \right. \quad \Leftrightarrow \quad \left\{ \begin{array}{c} \gamma_{11} (x_1 - 10) + \gamma_{12} (x_2 - x_1^2) + \gamma_{13} (x_3 - x_2^2) + \gamma_{14} (x_4 - x_3^2) = 0 \\ \gamma_{21} (x_1 - 10) + \gamma_{22} (x_2 - x_1^2) + \gamma_{23} (x_3 - x_2^2) + \gamma_{24} (x_4 - x_3^2) = 0 \\ \gamma_{31} (x_1 - 10) + \gamma_{32} (x_2 - x_1^2) + \gamma_{33} (x_3 - x_2^2) + \gamma_{34} (x_4 - x_3^2) = 0 \\ \gamma_{41} (x_1 - 10) + \gamma_{42} (x_2 - x_1^2) + \gamma_{43} (x_3 - x_2^2) + \gamma_{44} (x_4 - x_3^2) = 0 \\ \end{array} \right. \end{displaymath} where the $\gamma$-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. \section{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. \smallskip 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. \bigskip \noindent {\bf Acknowledgments}. This work could not have been started without Birk Huber with whom I first learned in~\cite{hv98} 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. {\small \begin{thebibliography}{10} \bibitem{ber75} D.N. Bernshte{\v{\i}}n. \newblock The number of roots of a system of equations. \newblock {\em Functional Anal. Appl.}, 9(3):183--185, 1975. \newblock Translated from {\em Funktsional. Anal. i Prilozhen.}, 9(3):1--4,1975. \bibitem{bcss97} L.~Blum, F.~Cucker, M.~Shub, and S.~Smale. \newblock {\em Complexity and Real Computation}. \newblock Springer--Verlag, New York, 1997. \bibitem{cd97} E.~Cattani and A.~Dickenstein. \newblock A global view of residues in the torus. \newblock {\em Journal of Pure and Applied Algebra}, 117--118:119--144, 1997. \newblock Proceedings of MEGA'96. \bibitem{cds96} E.~Cattani, A.~Dickenstein, and B.~Sturmfels. \newblock Computing multidimensional residues. \newblock In L.~Gonzalez-Vega and T.~Recio, editors, {\em Algorithms in Algebraic Geometry and Applications}, volume 143 of {\em Progress in Mathematics}, pages 135--164. Birkh\"auser, Basel, 1996. \newblock Proceedings of MEGA'94. \bibitem{cls98} D.~Cox, J.~Little, and D.~O'Shea. \newblock {\em Using Algebraic Geometry}, volume 185 of {\em Graduate Texts in Mathematics}. \newblock Springer--Verlag, New York, 1998. \bibitem{cox95} D.A. Cox. \newblock The homogeneous coordinate ring of a toric variety. \newblock {\em J. Algebraic Geometry}, 4(1):17--50, 1995. \bibitem{dan78} V.~Danilov. \newblock The geometry of toric varieties. \newblock {\em Russ. Math. Surveys}, 33:97--154, 1978. \bibitem{ded97} J.P. Dedieu. \newblock Condition number analysis for sparse polynomial systems. \newblock In F.~Cucker and M.~Shub, editors, {\em Foundations of Computational Mathematics. Selected Papers of a Conference, Held at IMPA in Rio de Janeiro, January 1997}, pages 75--101. Springer--Verlag, 1997. \bibitem{ds97} J.P. Dedieu and M.~Shub. \newblock Multihomogeneous {N}ewton methods. \newblock {IBM} {T}echnical {R}eport, IBM Research Division, 1997. \bibitem{ds98} J.P. Dedieu and M.~Shub. \newblock Newton and predictor-corrector methods for over-determined systems of equations. \newblock {IBM} {T}echnical {R}eport, IBM Research Division, 1998. \bibitem{ful93} W.~Fulton. \newblock {\em Introduction to Toric Varieties}, volume 131 of {\em Annals of Mathematics Studies}. \newblock Princeton University Press, Princeton, New Jersey, 1993. \bibitem{glvw97} T.~Gao, T.Y. Li, J~Verschelde, and M.~Wu. \newblock Balancing the lifting values to improve the numerical stability of polyhedral homotopy continuation methods. \newblock Manuscript, 1997. \bibitem{gkz94} I.M. {Gel'fand}, M.M. Kapranov, and A.V. Zelevinsky. \newblock {\em {D}iscriminants, {R}esultants and {M}ultidimensional {D}eterminants}. \newblock Birkh\"auser, Basel, 1994. \bibitem{hs95} B.~Huber and B.~Sturmfels. \newblock A polyhedral method for solving sparse polynomial systems. \newblock {\em Math. Comp.}, 64(212):1541--1555, 1995. \bibitem{hv98} B.~Huber and J.~Verschelde. \newblock Polyhedral end games for polynomial continuation. \newblock {\em Numerical Algorithms}, 18(1):91--108, 1998. \bibitem{kp96} T.~Krick and L.M. Pardo. \newblock A computational method for diophantine approximation. \newblock In L.~Gonzalez-Vega and T.~Recio, editors, {\em Algorithms in Algebraic Geometry and Applications}, volume 143 of {\em Progress in Mathematics}, pages 193--253. Birkh\"auser, Basel, 1996. \newblock Proceedings of MEGA'94. \bibitem{li97} T.Y. Li. \newblock Numerical solution of multivariate polynomial systems by homotopy continuation methods. \newblock {\em Acta Numerica}, 6:399--436, 1997. \bibitem{mm87} K.~Meintjes and A.P. Morgan. \newblock A methodology for solving chemical equilibrium systems. \newblock {\em Appl. Math. Comput.}, 22(4):333--361, 1987. \bibitem{mor87} A.~Morgan. \newblock {\em Solving polynomial systems using continuation for engineering and scientific problems}. \newblock Prentice-Hall, Englewood Cliffs, N.J., 1987. \bibitem{oda88} T.~Oda. \newblock {\em Convex {B}odies and {A}lgebraic {G}eometry: an {I}ntroduction to the {T}heory of {T}oric {V}arieties}. \newblock Springer--Verlag, Berlin Heidelberg New York, 1988. \bibitem{roj98} J.M. Rojas. \newblock Solving degenerate sparse polynomial systems faster. \newblock Manuscript available at http://www.cityu.edu.hk/ma/staff/rojas/. \bibitem{rs98} M.F. Roy and A.~Szpirglas. \newblock Bezoutiens et r{\'e}sidus, affines, projectifs et dans le tore. \newblock Manuscript, 1998. \bibitem{ver97} J.~Verschelde. \newblock {PHCPACK}: A general-purpose solver for polynomial systems by homotopy continuation. \newblock Paper and software available at the author's web-pages. \bibitem{ver98} J.~Verschelde. \newblock Numerical evidence for a conjecture in real algebraic geometry. \newblock Preprint \#1998-064, MSRI, 1998. \newblock Paper and software available at the author's web-pages. \bibitem{vgc96} J.~Verschelde, K.~Gatermann, and R.~Cools. \newblock Mixed-volume computation by dynamic lifting applied to polynomial system solving. \newblock {\em Discrete Comput. Geom.}, 16(1):69--112, 1996. \bibitem{vvc94} J.~Verschelde, P.~Verlinden, and R.~Cools. \newblock Homotopies exploiting {N}ewton polytopes for solving sparse polynomial systems. \newblock {\em SIAM J. Numer. Anal.}, 31(3):915--930, 1994. \bibitem{wbm87} L.T. Watson, S.C. Billups, and A.P. Morgan. \newblock Algorithm 652: {HOMPACK}: a suite of codes for globally convergent homotopy algorithms. \newblock {\em ACM Trans. Math. Softw.}, 13(3):281--310, 1987. \bibitem{wsmmw97} L.T. Watson, M.~Sosonkina, R.C. Melville, A.P. Morgan, and H.F. Walker. \newblock {HOMPACK90}: A suite of {F}ortran 90 codes for globally convergent homotopy algorithms. \newblock {\em ACM Trans. Math. Softw.}, 23(4):514--549, 1997. \newblock Available at http://www.cs.vt.edu/{\~{}}ltw/. \bibitem{zie95} G.M. Ziegler. \newblock {\em Lectures on {P}olytopes}, volume 152 of {\em Graduate Texts in Mathematics}. \newblock Springer--Verlag, New York, 1995. \end{thebibliography} } \end{document}