A Special Homotopy Continuation Method For A Class of Polynomial Systems

TL;DR

This paper introduces a novel homotopy continuation method combining polyhedral and linear product techniques to efficiently compute solutions of specific polynomial systems, with applications in real algebraic geometry.

Contribution

It proposes a new homotopy method with a root bound between total degree and mixed volume, and implements it as the LPH program for improved efficiency.

Findings

The method computes all isolated solutions efficiently.

The root number bound is easily calculated and tighter than traditional bounds.

The algorithm can find witness points on real variety components.

Abstract

A special homotopy continuation method, as a combination of the polyhedral homotopy and the linear product homotopy, is proposed for computing all the isolated solutions to a special class of polynomial systems. The root number bound of this method is between the total degree bound and the mixed volume bound and can be easily computed. The new algorithm has been implemented as a program called LPH using C++. Our experiments show its efficiency compared to the polyhedral or other homotopies on such systems. As an application, the algorithm can be used to find witness points on each connected component of a real variety.

Tables

Table 1. Table 1: Dense Examples

$n$	$k$	T1	T2	RAT
2	1	0.125s	0.094s	1.32
3	1	0.125s	0.109s	1.14
3	2	0.125s	0.109s	1.14
4	1	0.125s	0.109s	1.14
4	2	0.156s	0.156s	1.00
4	3	0.202s	0.265s	0.76
5	1	0.125s	0.109s	1.14
5	2	0.187s	0.202s	0.93
5	3	0.390s	0.655s	0.59
5	4	0.687s	1.280s	0.54
6	1	0.140s	0.109s	1.28
6	2	0.281s	0.328s	0.86
6	3	0.76s	1.68s	0.45
6	4	1.61s	4.36s	0.37
6	5	1.90s	6.59s	0.29
7	1	0.14s	0.10s	1.29
7	2	0.344s	0.56s	0.61
7	3	1.3s	3.82s	0.34
7	4	3.7s	14.1s	0.27
7	5	6.318s	27.9s	0.23
7	6	6.006s	27.6s	0.21
8	1	0.15s	0.15s	1.00
8	2	0.54s	0.73s	0.74
8	3	2.29s	7.2s	0.31
8	4	8.018s	34.1s	0.23
8	5	17.6s	91.2s	0.19
8	6	23.7s	153s	0.15
8	7	19.4s	128s	0.15
9	1	0.2s	0.18s	1.08
9	2	0.7s	1.2s	0.58
9	3	4.1s	13.1s	0.31
9	4	16.1s	1m19s	0.20
9	5	46.5s	4m29s	0.17
9	6	1m20s	9m38s	0.138
9	7	1m30s	11m10s	0.135
9	8	59.9s	7m52s	0.126
10	1	0.23s	0.18s	1.25
10	2	0.98s	1.9s	0.51
10	3	5.8s	24.8s	0.24
10	4	31.1s	2m53s	0.18
10	5	1m46s	11m30s	0.15
10	6	3m41s	29m15s	0.13
10	7	5m10s	48m32s	0.107
10	8	4m57s	48m27s	0.102
10	9	3m2s	30m7.553s	0.1

Table 2. Table 2: Sparse Examples

Ex	$n$	$k$	$d$	term	#1	#2	#	T1	T2	RAT
C2	5	4	4-5	7-32	2*1767	692	383	16.6s	19.4s	0.85
M3	9	5	2	2-7	2*2240	368	32	69s	6s	11
G2	5	2	4	8-9	2*1080	17	15	8.4s	0.2s	31.8
H1	8	6	1-3	2-4	2*400	15	15	9.8s	0.18s	52
H2	8	5	2-4	3-5	2*12320	148	80	5m49s	1.2s	267

Equations

$$F(x,\lambda):=\left[{\begin{array}[]{*{20}{c}}f\\ {{\lambda_{0}}\nabla\Phi{{(x)}^{T}}+{\lambda_{1}}\nabla f_{1}^{T}+...+{\lambda_{k}}\nabla f_{k}^{T}}\end{array}}\right]$$

$$F(x,\lambda )=\{\begin{array}{c}f\text{hfill}\\ J\cdot \lambda -\beta \text{hfill}\end{array}$$

$$\left({\begin{array}[]{*{20}{c}}{n-1}\\ {n-k}\end{array}}\right){d^{n-k}}D$$

$$F=\left\{f,\sum _{i=1}^k{\lambda }_i\nabla f_i-\beta \right\}=0.$$

$$f^{(1)}=\left\{f,x\cdot \beta +c\right\}.$$

$$\left({\begin{array}[]{*{20}{c}}{n-1}\\ {n-k}\end{array}}\right){(d-1)^{n-k}}D$$

$$\sum\limits_{j=0}^{n-k}{\left({\begin{array}[]{*{20}{c}}{n-1-j}\\ {n-k-j}\end{array}}\right){{(d-1)}^{n-k-j}}D}.$$

$$MV(F)\leqslant\left({\begin{array}[]{*{20}{c}}{n-1}\\ {k-1}\end{array}}\right){(d-1)^{n-k}}D\leqslant\left({\begin{array}[]{*{20}{c}}{n-1}\\ {k-1}\end{array}}\right){(d-1)^{n-k}}\prod\limits_{i=1}^{k}{{d_{i}}}\leqslant{d^{n}}\prod\limits_{i=1}^{k}{{d_{i}}}.$$

$${V_{\bbbr}}=\left\{{\begin{array}[]{*{20}{c}}{(-1.44299,-1.32941),}&{(-0.781143,1.28371)}\end{array}}\right\}$$

$${W_{\bbbr}}({f^{(1)}})=\left\{{\begin{array}[]{*{20}{c}}{(2.4052801,1.815026),}&{(-1.992641,5.531208)}.\end{array}}\right\}$$

$$H_2(x,\lambda ,t)=G\cdot (1-t)+F\cdot {\gamma }_2\cdot t.$$

Full text

1 Introduction

In many applications in science, engineering, and economics, solving systems of polynomial equations has been a subject of great importance. The homotopy continuation method was developed in 1970s [1][2] and has been greatly expanded and developed by many reseachers (see for example [3][4][5][6][7]). Nowadays, homotopy continuation method has become one of the most reliable and efficient classes of numerical methods for finding the isolated solutions to a polynomial system and the so-called numerical algebraic geometry based on homotopy continuation method has been a blossoming area. There are many famous software packages implementing different homotopy methods, including Bertini[8], Hom4PS-2.0[9], HOMPACK[10], PHCpack[11], etc.

Classical homotopy methods compute solutions in complex spaces, while in applications, it is quite common that only real solutions have physical meaning. Computing real roots of an algebraic system is a difficult and fundamental problem in real algebraic geometry. In the field of symbolic computation, there are some famous algorithms dealing with this problem. The cylindrical algebraic decomposition algorithm [12] is the first complete algorithm which has been implemented and used successfully to solve many real problems. However, in the worst case, its complexity is of doubly exponential in the number of variables. Based on the ideas of Seidenberg [13] and others, some algorithms for computing at least one point on each connected component of an real algebraic set were proposed through developing the formulation of critical points and the notion of polar varieties, see [14][15][16] and references therein. The idea behind is studying an objective function (or map) that reaches at least one local extremum on each connected component of a real algebraic set. For example, the function of square of the Euclidean distance to a randomly chosen point was used in [17][18]. On the other hand, some homotopy based algorithms for real solving have been proposed in [19][20][21][22][23][24]. For example, in [23], a numerical homotopy method to find the extremum of Euclidean distance to a point as the objective function was presented. More recently, the Euclidean distance to a plane was proposed as a linear objective function in [25].

In this paper, we follow the work of [25][26] to extend complex homotopy methods to finding witness points on the irreducible components of real varieties. To obtain such witness points, we first need to solve a special class of polynomial systems. Combining the polyhedral homotopy and the linear product homotopy, we give a special homotopy method for solving the system of that type. The root number bound of this method is not only easy to compute but also much smaller than the total degree bound and close to the BKK bound [27] when the polynomials defining the algebraic set is not very sparse. This key observation enables us to design an efficient homotopy procedure to obtain critical points numerically. The ideas and algorithms we proposed in this article avoid a great number of divergent paths to track compared with the total degree homotopy and save the great time cost for mixed volume computation compared with the polyhedral homotopy. The new algorithm has been implemented as a program called LPH using C++. Our experiments show its efficiency compared to the polyhedral or other homotopies on such systems.

The rest of this paper is organized as follows. Section 2 describes some preliminary concepts and results. Section 3 introduces a special type of polynomial systems we are considering. The new homotopy for these polynomial systems is also presented. It naturally leads to an algorithm which is described in Section 4. Based on this algorithm, in Section 5, we present a method to find real witness points of positive dimensional varieties, together with an illustrative example. The experimental performance of the software package LPH, which is an implementation of the method in C++, is given in Section 6.

2 Preliminary

2.1 Algebraic Sets and Genericity

For a polynomial system $f:{\bbbc^{n}}\to{\bbbc^{k}}$, let $V(f)=\{x\in{\bbbc^{n}}|f(x)=0\}$ and ${V_{\bbbr}}(f)=V(f)\cap{\bbbr^{n}}=\{x\in{\bbbr^{n}}|f(x)=0\}$ be the set of complex solutions and the set of real solutions of $f(x)=0$, respectively. A set $X\subset{\bbbc^{n}}$ is called an algebraic set if $X=V(g)$, for some polynomial system $g$.

An algebraic set $X$ is irreducible if there does not exist a decomposition ${X_{1}}\cup{X_{2}}=X$ with ${X_{1}},{X_{2}}\neq X$ of $X$ as a union of two strict algebraic subsets. An algebraic set is reducible, if there exist such decomposition. For example, the algebraic set $V(xy)\subset{\bbbc^{2}}$ is consisting of the two coordinate axes, and is obviously the union of $V(x)$ and $V(y)$, hence reducible.

For an irreducible algebraic set $X$, the subset of smooth (or manifold) points ${X_{reg}}$ is dense, open and path connected (up to the Zariski topology) in $X$. The dimension of an irreducible algebraic set $X$ is the dimension of ${X_{reg}}$ as a complex manifold.

Let $\mathcal{J}_{f}(x)$ denote the $n\times k$ Jacobian matrix of $f$ evaluated at $x$. By the Implicit Function Theorem, for an irreducible algebraic set $X$ defined by a reduced system $f$, $x\in{X_{reg}}\Leftrightarrow\mathrm{rank}(\mathcal{J}_{f}(x))=n-\dim X$. When $n=k$, the system $f$ is said to be a square system. In this case, a point $x\in V(f)$ is nonsingular if $\det(\mathcal{J}_{f}(x))\neq 0$, and singular otherwise.

On irreducible algebraic set, we can define the notion of genericity, adapted from[6].

Definition 1

Let $X$ be an irreducible algebraic set. Property P holds generically on $X$, if the set of points in $X$ that do not satisfy property P are contained in a proper algebraic subset $Y$ of $X$. The points in $Y$ are called nongeneric points, and their complements $X\backslash Y$ are called generic points.

Remark 1

From the definition, one sees that the notion of generic is only meaningful in the context of property P in question.

Every algebraic set $X$ has a (uniquely up to reordering) expression $X={X_{1}}\cup\ldots\cup{X_{r}}$ with ${X_{i}}$ irreducible and ${X_{i}}\not\subset{X_{j}}$ for $i\neq j$. And ${X_{i}}$ are the irreducible components of $X$. The dimension of an algebraic set is defined to be the maximum dimension of its irreducible components. An algebraic set is said to be pure-dimensional if each of its components has the same dimension.

2.2 Trackable Paths

In homotopy continuation methods, the notion of path tracking is fundamental, the following definition of trackable solution path is adapted from [28].

Definition 2

Let $H(x,t):{\bbbc^{n}}\times\bbbc\to{\bbbc^{n}}$ be polynomial in $x$ and complex analytic in $t$, and let ${x^{*}}$ be nonsingular isolated solution of $H(x,0)=0$, we say ${x^{*}}$ is trackable for $t\in[0,1)$ from 0 to 1 using $H(x,t)$ if there is a smooth map ${\xi_{{x^{*}}}}:[0,1)\to{\bbbc^{n}}$ such that ${\xi_{{x^{*}}}}(0)={x^{*}}$, and for $t\in[0,1)$, ${\xi_{{x^{*}}}}(t)$ is a nonsingular isolated solution of $H(x,t)=0$. The solution path started at ${x^{*}}$is said to be convergent if $\mathop{\lim}\limits_{t\to 1}{\xi_{{x^{*}}}}(t)\in{\bbbc^{n}}$, and the limit is called the endpoint of the path.

2.3 Witness Set and Degree of an Algebraic Set

Let $X\subset{\bbbc^{n}}$ be a pure $i$-dimensional algebraic set, given a generic co-dimension $i$ affine linear subspace $L\subset{\bbbc^{n}}$, then $W=L\cap X$ consists of a well-defined number $d$ of points lying in ${X_{reg}}$. The number $d$ is called the degree of $X$ and denoted by $\deg(X)$. We refer to $W$ as a set of witness points of $X$, and call $L$ the associated $(n-i)$-slicing plane, or slicing plane for short [6].

It will be convenient to use the notations adapted from ([6],Chapter 8), when we prove the theorems in Section 3.

Let $\left\langle{{e_{1}},\ldots,{e_{n}}}\right\rangle$ be the $n$ dimensional vector space having basis elements ${e_{1}},\ldots,{e_{n}}$ with complex coefficients. That is, a point in this space may be written as $\sum\limits_{i=1}^{n}{{c_{i}}{e_{i}}}$, with ${c_{i}}\in\bbbc$ for $i=1,\ldots,n$. Note that we have not specified anything about the basis elements, it could be individual variables, monomials, or polynomials. 2. 2.

Let $\left\{{{p_{1}},\ldots,{p_{n}}}\right\}\otimes\left\{{{q_{1}},\ldots,{q_{m}}}\right\}$ be the product of two sets, that is, the set $\left\{{{p_{i}}\cdot{q_{j}}|i=1,\ldots,n;j=1,\ldots,m}\right\}$. In Section 3 we take this product as the image inside the ring of polynomials; that is, $x\otimes y=xy$ is just the product of two polynomials. 3. 3.

Define $P\times Q=\{pq|p\in P,q\in Q\}$. Accordingly, we have $\left\langle P\right\rangle\times\left\langle Q\right\rangle\subset\left\langle{P\otimes Q}\right\rangle$. 4. 4.

For repeated products, we use the shorthand notations ${P^{(2)}}=P\otimes P$,${\left\langle P\right\rangle^{(2)}}=\left\langle P\right\rangle\times\left\langle P\right\rangle$, and similar for three or more products. 5. 5.

For a square polynomial system $P$, we denote by $MV(P)$ the mixed volume of the system $P$.

2.4 Critical Points

Let $X\subset{\bbbc^{n}}$ be an algebraic set defined by a reduced polynomial system $f=\{{f_{1}},\ldots,{f_{k}}\}$, and objective function $\Phi$ is polynomial function restricted to $X$.

Definition 3

A point $x\in X$ is a critical point of $\Phi$ if and only if $x\in{X_{reg}}$ and $\mathrm{rank}[\nabla\Phi{(x)^{T}},\mathcal{J}_{f}(x)]=\mathrm{rank}[\nabla\Phi{(x)^{T}},\nabla f_{1}^{T},...,\nabla f_{k}^{T}]\leqslant k$, where $\nabla\Phi(x)$ is the gradient vector of $\Phi$ evaluated at $x$.

Let $Y$ denote the zero dimensional critical sets of $\Phi$. One way to compute the critical points is to introduce auxiliary unknowns and consider a zero dimensional variety $\hat{Y}$ and then project $\hat{Y}$ onto $Y$. We use Lagrange Multipliers to define a squared system as follows

[TABLE]

Note that if ${x^{*}}\in X$ is a critical point of $\Phi$, then there exist ${\lambda^{*}}\in{\bbbp^{k}}$, such that $F({x^{*}},{\lambda^{*}})=0$ by the Fritz John condition[29]. In the affine patch where ${\lambda_{0}}=1$, the system $F$ becomes a square system, and its solution $({x^{*}},{\lambda^{*}})$ projects to critical point ${x^{*}}$. We will use system (1) in Section 5 with an objective function $\Phi$ defined by a linear function, and consider the affine patch where ${\lambda_{0}}=1$, to find at least one point on each component of ${V_{\bbbr}}(f)$.

3 Main Idea

In this section, we give a description of our idea. First we introduce a family of polynomial equations that we will be considering.

We consider the following class of polynomial systems:

[TABLE]

where

$f=\left\{{{f_{1}},\ldots,{f_{k}}}\right\}$ are polynomials in $\bbbc\left[{{x_{1}},\ldots,{x_{n}}}\right]$, and $V({f_{1}},\ldots,{f_{k}})$ is a pure $n-k$ dimension algebraic set in ${\bbbc^{n}}$. 2. 2.

J=\left({\begin{array}[]{*{20}{c}}{{g_{11}}}&\cdots&{{g_{1k}}}\\ \vdots&\ddots&\vdots\\ {{g_{n1}}}&\cdots&{{g_{nk}}}\end{array}}\right) and ${g_{ij}}(1\leqslant i\leqslant n,1\leqslant j\leqslant k)$ are polynomials in $\bbbc\left[{{x_{1}},\ldots,{x_{n}}}\right]$ with $\mathop{\max}\limits_{i,j}\deg({g_{ij}})=d$. 3. 3.

$\beta={({\beta_{1}},\ldots,{\beta_{n}})^{\rm T}}$ is a nonzero constant vector in ${\bbbc^{n}}$, $\lambda=({\lambda_{1}},\ldots,{\lambda_{k}})^{\rm T}$ are unknowns, and $n>k\geqslant 1$.

Remark 2

Note that, for any invertible $n\times n$ matrix $A$, $F(x,\lambda)=\left\{{f,J\cdot\lambda-\beta}\right\}$ and $F^{\prime}(x,\lambda)=\left\{{f,A\cdot(J\cdot\lambda-\beta)}\right\}$ have the same solutions. It’s easy to know that there exists an invertible matrix $A$ such that $A\cdot\beta={(0,\ldots,0,1)^{\rm T}}$. So without loss of generality, we may assume that $\beta={(0,\ldots,0,1)^{\rm T}}$. Then, ${J\cdot\lambda-\beta}$ has $n-1$ equations in $\left\langle{{{\{{x_{1}},\ldots,{x_{n}},1\}}^{d}}\otimes\{{\lambda_{1}},\ldots,{\lambda_{k}}\}}\right\rangle$ and one equation in $\left\langle{{{\{{x_{1}},\ldots,{x_{n}},1\}}^{d}}\otimes\{{\lambda_{1}},\ldots,{\lambda_{k}},1\}}\right\rangle$.

Theorem 3.1

Let $F(x,\lambda)=\left\{{f,J\cdot\lambda-\beta}\right\}$ be given as in (2), $\beta={(0,\ldots,0,1)^{\rm T}}$, and $G=\left\{{f,g}\right\}$ where $g=\left\{{g_{1}},\ldots,{g_{n}}\right\}$. ${g_{i}}={l_{i1}}\cdots{l_{id}}{h_{i}}\in{\left\langle{{x_{1}},\ldots,{x_{n}},1}\right\rangle^{d}}\times\left\langle{{\lambda_{1}},\ldots,{\lambda_{k}}}\right\rangle$ for $i=1,\ldots,n-1$ ; where ${l_{ij}}$ are linear functions in $\bbbc[{x_{1}},\ldots,{x_{n}}]$, ${h_{i}}$ with randomly chosen coefficient and ${h_{i}}$ are homogeneous linear functions in $\bbbc[{\lambda_{1}},\ldots,{\lambda_{k}}]$, $i=1,\ldots,n-1$, $j=1,\ldots,d$ and ${g_{n}}=\sum\limits_{i=1}^{k}{{\lambda_{i}}}{g_{ni}}-1$. $H:{\bbbc^{n}}\times{\bbbc^{k}}\times\bbbc\to{\bbbc^{n+k}}$ be the homotopy defined by $H(x,\lambda,t)=G\cdot(1-t)+F\cdot\gamma\cdot t$ where $\gamma$ is a randomly chosen complex number for Gamma Trick (see [6] Chapter 7 for details). Then, generically the following items hold,

The set $S\subseteq{\bbbc^{n+k}}$ of roots of $H(x,\lambda,0)=G(x,\lambda)$ is finite and each is a nonsingular solution of $H(x,\lambda,0)$. 2. 2.

The number of points in $S$ is equal to the maximum number of isolated solutions of $H(x,\lambda,0)$ as coefficients of ${l_{ij}}$, ${h_{i}}$, ($i=1,\ldots,n-1$, $j=1,\ldots,d$) and $\gamma$ vary over $\bbbc$. 3. 3.

The solution paths defined by $H$ starting, with $t=0$, at the points in $S$ are trackable.

Proof

As for item 1, since $f$ has $k$ equations only in $x$, and $V({f_{1}},\ldots,{f_{k}})$ is a pure $n-k$ dimension algebraic set in ${\bbbc^{n}}$. To solve system $G$, it needs only $n-k$ linear functions $L$ in $g$ from different ${g_{i}}$ with $i\in\left\{{1,\ldots,n-1}\right\}$ to determine $x$. $\left\{{f,L}\right\}$ is a $n\times n$ square system, $V(f,L)$ is a finite witness set for algebraic set $V({f_{1}},\ldots,{f_{k}})$, and each of the points is a nonsingular solution of $V(f,L)$ (see [6] Chapter 13 for details). And, we finally determine $\lambda$ by solving a square linear equations. As for item 2, and item 3, it’s a trivial deduction of Coefficient-Parameter Continuation [30].∎

Remark 3

From the proof of Theorem 3.1, the number of points of the finite set $V(f,L)$ is the degree of $V(f)$, and is independent of the choice of $L$. Thus, based on the number of different choices of $L$, and item 2, we can give a root count bound of system $F(x,\lambda)$ as in the following theorem, which is similar to the bound in [26].

Theorem 3.2

For a system $F(x,\lambda)=\left\{{f,J\cdot\lambda-\beta}\right\}$ as in (2). The number of complex root of this system is bounded by

[TABLE]

where $D$ is the degree of $V(f)$.

Due to Theorem 3.1, its proof and the remarks, we can design an efficient procedure to numerically find the isolated solutions of system $F(x,\lambda)=\left\{{f,J\cdot\lambda-\beta}\right\}$ in the form of (2). First, we solve a square systems $\left\{{f,L}\right\}$, where $L$ are $n-k$ randomly generated linear functions . Then for each group of $n-k$ linear functions $L^{\prime}$ chosen in $g$ from different ${g_{i}}$ with $i\in\left\{{1,\ldots,n-1}\right\}$, we construct linear homotopy from $\left\{{f,L}\right\}$ to $\left\{{f,L^{\prime}}\right\}$, starting from points of $V(f,L)$, and solve the square linear equation of $\lambda$ respectively. Let $S$ be the set consist of all the pairs of $x$ and $\lambda$, i.e. $(x,\lambda)$. Finally construct linear homotopy $H(x,\lambda,t)=G\cdot(1-t)+F\cdot\gamma\cdot t$ starting from points in $S$, thus the endpoints of the convergent paths of homotopy $H(x,\lambda,t)$ are isolated solutions of system $F(x,\lambda)=\left\{{f,J\cdot\lambda-\beta}\right\}$. We put specific description of this procedure in the next section.

4 Algorithm

From Theorem 3.1, its proof and Remarks 2 & 3, we propose an approach for computing isolated solutions of system $F(x,\lambda)$ as described in the end of last section. For consideration of the sparsity, we use the polyhedral homotopy method for solutions of the square system $\left\{{f,L}\right\}$. Actually, we use polyhedral homotopy method only once. Now we describe our algorithms.

Remark 4

\#C=\left({\begin{array}[]{*{20}{c}}{n-1}\\ {n-k}\end{array}}\right) , and in Step 5, $I=({\alpha_{1}},\ldots,{\alpha_{n-1}})$ has and only has $n-k$ entries ${\alpha_{i}}=1$. When ${\alpha_{i}}=1$, we choose linear equation in ${g_{i}}$, and there are $d$ candidates $\{{l_{i1}},\ldots,{l_{id}}\}$ to choose. It adds up to be \left({\begin{array}[]{*{20}{c}}{n-1}\\ {n-k}\end{array}}\right){d^{n-k}} different $\left\{{f,L^{\prime}}\right\}$. Each $\left\{{f,L^{\prime}}\right\}$ has the same number $D=\deg(V(f))$ of isolated roots as $\left\{{f,L}\right\}$, so homotopy in Step 13 will have no path divergent. Thus we have \left({\begin{array}[]{*{20}{c}}{n-1}\\ {n-k}\end{array}}\right){d^{n-k}}D points in $\Omega$, which is the root bound we mention in Theorem 3.2. It would happen that some of the homotopy paths divergent in Step 16, the method of end games for homotopy should be used [31][32][33][34].

5 Real Critical Set

In this section, we will combine the LPH algorithm in Section 4 and methods in [25] to compute a real witness set which has at least one point on each irreducible component of a real algebraic set, and give an illustrative example.

5.1 Critical Points on a Real Algebraic Set

We make the following assumptions (adapted from [25]). Let $f:{\bbbc^{n}}\to{\bbbc^{k}}$ be a polynomial system, and $f=({f_{1}},\ldots,{f_{k}})$ in $\bbbr[{x_{1}},\ldots,{x_{n}}]$ satisfying the so-called Full Rank Assumption:

${V_{\bbbr}}({f_{1}},\ldots,{f_{i}})$ has dimension $n-i$ for $i=1,\ldots,k$; 2. 2.

the ideal $I({f_{1}},\ldots,{f_{i}})$ is radical for $i=1,\ldots,k$.

Under these assumptions, $(\nabla f_{1}^{T},\ldots,\nabla f_{i}^{T})$ has rank $i$ for a generic point $p\in V({f_{1}},\ldots,{f_{i}})$ for $i=1,\ldots,k$.

The main problem we consider is finding at least one real witness point on each real dimensional components of ${V_{\bbbr}}(f)$. For this purpose, we choose $\Phi$ in Definition 3 to be a linear function with $\Phi=x\cdot\beta+c$, where $\beta$ is a random vector in ${\bbbr^{n}}$, and $c$ is a random real number. Then system (1) becomes

[TABLE]

It may happen that there is no critical points of $\Phi$ in some connected component of ${V_{\bbbr}}({f_{1}},\ldots,{f_{k}})$. In that case, we add $\Phi$ to $f$ and construct a system with $k+1$ equations

[TABLE]

Then recursively, we choose another linear function $\Phi_{1}$, compute the critical points of $\Phi_{1}$ with respect to $V({f^{(1)}})$; and so on.

We give a concrete definition of the set of real witness points ${W_{\bbbr}}(f)$ we are going to compute (see [25]).

Definition 4

Let $f:{\bbbc^{n}}\to{\bbbc^{k}}$ be a polynomial system, $k\leqslant n$, and $f=({f_{1}},\ldots,{f_{k}})$ in $\bbbr[{x_{1}},\ldots,{x_{n}}]$ satisfying Full Rank Assumption. $F$ and ${f^{(1)}}$ defined as (4) and (5). We define ${W_{\bbbr}}(f)$ as follows:

${W_{\bbbr}}(f)={V_{\bbbr}}(f)$ if $n=k$; 2. 2.

${W_{\bbbr}}(f)={V_{\bbbr}}(F)\cup{W_{\bbbr}}({f^{(1)}})$ if $k<n$.

It is obvious from the definition that we can recursively solve the square system (4), and apply plane distance critical points formulation of ${f^{(1)}}$ to finally get the set of witness points ${W_{\bbbr}}(f)$ which contains finitely many real points on ${V_{\bbbr}}(f)$, and there is at least one point on each connected component of ${V_{\bbbr}}(f)$. Since the formulation introduces auxiliary unknowns, it increases the size of the system and leads to computational difficulties. For example, when $n=15$ and $k=10$, the size of system (4) becomes $25$, which is challenging for general homotopy software. Combining the LPH algorithm, Theorem 3.1 and Theorem 3.2, we have the following algorithm and an upper bound of number of points in ${W_{\bbbr}}(f)$, as in [26].

Remark 5

Algorithm 2 is essentially a recursive calling of Algorithm 1.

Theorem 5.1

([26] Theorem 2.1) For a system $f=({f_{1}},\ldots,{f_{k}})$ with $n$ variables and degrees ${d_{i}}=\deg({f_{i}})$ for $i=1,\ldots,k$. The number of complex root of system (4) is bounded by

[TABLE]

*where $d=\max\{{d_{1}},\ldots,{d_{k}}\}>1$ and $n>k>0$, $D$ is the degree of the pure $n-k$ dimensional component of $V=V(f)$.

Moreover, the total number of points in ${W_{\bbbr}}(f)$ is bounded by*

[TABLE]

Obviously we have the following inequalities:

[TABLE]

If $f$ is dense, the equalities hold. And if $f$ is sparse, they vary considerably most of the time. For example, let $f=\{-62xy+97y-4xyz-4,80x-44xy+71{y^{2}}-17{y^{3}}+2\}$ with $d=3,n=3,k=2$. We have $MV(F)=11$, \left({\begin{array}[]{*{20}{c}}{n-1}\\ {k-1}\end{array}}\right){(d-1)^{n-k}}D=28, \left({\begin{array}[]{*{20}{c}}{n-1}\\ {k-1}\end{array}}\right){(d-1)^{n-k}}\prod\limits_{i=1}^{k}{{d_{i}}}=36, and ${d^{n}}\prod\limits_{i=1}^{k}{{d_{i}}}=243$.

5.2 Illustrative Example

In this subsection, we present an illustrative example for Algorithm 2.

Example 1

Consider the hypersurface defined by $f=({y^{2}}-{x^{3}}-ax-b)\cdot({(x-y+e)^{3}}+x+y)$, $e=6,a=-4,b=-1$. Clearly, ${V_{\bbbr}}(f)$ is the combination of a cubic ellipse $({y^{2}}-{x^{3}}-ax-b)$, and a cubic curve ${(x-y+e)^{3}}+x+y$, as plotted in Fig. 1. We show how to compute ${W_{\bbbr}}(f)$ by Algorithm 2.

• –

For computing ${V_{\bbbr}}=\texttt{LPH}(f)$, we randomly choose a line $l$ in ${\bbbc^{2}}$ and solve $L=\left\{{f,l}\right\}$ by polyhedral homotopy, which follows $D=6$ paths. Then to compute $\Omega$ by linear homotopy, we follow \left({\begin{array}[]{*{20}{c}}{2-1}\\ {2-1}\end{array}}\right){(6-1)^{2-1}}6=30 convergent paths, and for ${V_{\bbbr}}$ by linear homotopy, we follow 30 paths, of which 6 are convergent and 19 divergent. Then

[TABLE]

• –

For computing ${W_{\bbbr}}(f)$, we solve ${f^{(1)}}=\left\{{f,x\cdot\beta+c}\right\}$ by polyhedral homotopy, with $x\cdot\beta+c=0.874645x+1.0351y-3.9825$ and

[TABLE]

So ${W_{\bbbr}}(f)={W_{\bbbr}}({f^{(1)}})\cup{V_{\bbbr}}$, which has at least one point in each connected component of ${V_{\bbbr}}(f)$ as in Fig. 1.

6 Experiment Performance

As shown in Section 5, to compute the set ${W_{\bbbr}}(f)$, the key and most time consuming steps are solving the system $F=\left\{{f,\sum\limits_{i=1}^{k}{{\lambda_{i}}\nabla{f_{i}}}-\beta}\right\}$ in Algorithm 1. In this section, given $f=\left\{{{f_{1}},\ldots,{f_{k}}}\right\}$, we solve the square system $F=\left\{{f,\sum\limits_{i=1}^{k}{{\lambda_{i}}\nabla{f_{i}}}-\beta}\right\}$. We compare our program LPH which implements Algorithm 1 to Hom4PS-2.0 (available at http://www.math.msu.edu/ li). All the examples were computed on a PC with Intel Core i5 processor (2.5GHz CPU, 4 Cores and 6 GB RAM) in the Windows environment. We mention that LPH is a program written in C++, available at http://arcnl.org/PDF/LHP.zip, and an interface of Maple is provided on this site.

6.1 Dense Examples

In Table 1, we provide the timings of LPH and Hom4ps-2.0 for solving systems $F=\left\{{f,\sum\limits_{i=1}^{k}{{\lambda_{i}}\nabla{f_{i}}}-\beta}\right\}$, where $f=({f_{1}},\ldots,{f_{k}})$ consists of dense polynomials of degree 2, $n=2,...,14$ and $1\leqslant k\leqslant n-1$. T1 ,T2 are the the timings for LPH and Hom4ps-2.0, respectively, and RAT is the ratio of T1 to T2. ‘overflow’ means running out of memory. When T2=overflow, we set RAT=$\varepsilon$.

It may be observed that LPH is much faster than Hom4ps-2.0 when $k>1$. Note also that LPH is a little bit slower than Hom4ps-2.0 when $k=1$. The main reason is obvious. That is, the root number bound of LPH, i.e. \left({\begin{array}[]{*{20}{c}}{n-1}\\ {n-k}\end{array}}\right){(d-1)^{n-k}}D, is close to the mixed volume $MV(F)$ when $F$ is dense but the computation of $MV(F)$ is very time-consuming.

6.2 Sparse Examples

In Table 2, we provide the timings of LPH and Hom4ps-2.0 on sparse examples: Czapor Geddes2, Morgenstern AS(3or), Gerdt2, Hairer1, and Hawes2 which are available at : http://www-sop.inria.fr/saga/POL/. #1 and #2 is the number of curves followed by LPH and Hom4ps-2.0, respectively. # is the number of roots of the Jacobian systems constructed from the examples. ‘$d$’ means the minimal and maximal degree of the example. “term” means the minimal and maximal number of terms of the example. T1 and T2 are the timings of LPH and Hom4ps-2.0, respectively. RAT means the ratio of T1 to T2.

Note that LPH is much slower than Hom4ps on these sparse examples. The main reason is that LPH pays the overhead cost for the $\Omega$ and homotopy

[TABLE]

Moreover, LPH executes 2*\left({\begin{array}[]{*{20}{c}}{n-1}\\ {k-1}\end{array}}\right){(d-1)^{n-k}}D times of curve following, while Hom4ps does only $MV(F)$ times of curve following. \left({\begin{array}[]{*{20}{c}}{n-1}\\ {k-1}\end{array}}\right){(d-1)^{n-k}}D is not tight for these sparse examples and much greater than $MV(F)$.

6.3 RAT/Density

In Fig. 2, we present the changes of ratio of T1 to T2 as terms increase. We randomly generate $f=({f_{1}},\ldots,{f_{k}})$ with different $n,k$ and degrees, and increase the number of terms from 2 to dense.

It can be observed that, when the polynomials are not very sparse, e.g. the number of terms are more than $10\%$ of \left({\begin{array}[]{*{20}{c}}{n+d}\\ {d}\end{array}}\right), LPH is faster than Hom4ps-2.0. Actually, when the polynomials are not very sparse, the root number bound \left({\begin{array}[]{*{20}{c}}{n-1}\\ {k-1}\end{array}}\right){(d-1)^{n-k}}D is close to $MV(F)$.

7 Acknowledgement

We gratefully acknowledge the very helpful suggestions of Hoon Hong on this paper with emphasize on Section 6. We also thank Changbo Chen for his helpful comments.