Mixed Precision Path Tracking for Polynomial Homotopy Continuation

TL;DR

This paper introduces a novel predictor-corrector algorithm with adaptive step size and mixed precision arithmetic for polynomial homotopy continuation, improving robustness and efficiency in challenging numerical scenarios.

Contribution

It presents a new predictor-corrector method with a Newton corrector that rejects poor guesses and employs mixed precision to enhance numerical stability.

Findings

Demonstrates improved robustness in challenging cases

Shows efficiency gains through adaptive step size

Validates approach with multiple numerical examples

Abstract

This article develops a new predictor-corrector algorithm for numerical path tracking in the context of polynomial homotopy continuation. In the corrector step it uses a newly developed Newton corrector algorithm which rejects an initial guess if it is not an approximate zero. The algorithm also uses an adaptive step size control which builds on a local understanding of the region of convergence of Newton's method and the distance to the closest singularity. To handle numerically challenging situations the algorithm uses mixed precision arithmetic. The efficiency and robustness are demonstrated in several numerical examples.

Tables

Table 1. Table 1. Results of the experiment of Subsection 5.1 for ρ = 10 − k 𝜌 superscript 10 𝑘 \rho=10^{-k} and k = 1 , … , 7 𝑘 1 … 7 k=1,\ldots,7 . A ‘✗’ indicates that path jumping happened.

	1	2	3	4	5	6	7
Bertini DP	✓	✓	✓	✗	✗	✗	✗
Bertini AP	✓	✓	✓	✓	✓	✗	✗
phc -u	✓	✓	✓	✓	✓	✓	✓
HC.jl 1.4	✓	✗	✗	✗	✗	✗	✗
HC.jl	✓	✓	✓	✓	✓	✓	✓

Table 2. Table 2. Results for 10 runs of Alt’s problem using the same generic start solutions to a generic physically meaningful configuration. The tol column refers to the assigned path tracking tolerance.

		runtime (seconds)					# solutions
	tol	mean	median	min	max		median	min	max
HC.jl	-	6.38	6.41	5.08	7.12		1442	1442	1442
Bertini AP	default	1931.63	1832.59	1120.43	2695.13		1436	1433	1441
Bertini AP	1e-8	10950.09	11218.97	8371.76	12115.66		1441	1437	1442

Table 3. Table 3. Results for 50 runs of Steiner’s problem from generic complex parameter values to generic real parameter values. The tol column refers to the assigned path tracking tolerance.

		runtime (seconds)					# solutions
	tol	mean	median	min	max		median	min	max
HC.jl		2.47	2.49	1.64	3.11		3264	3264	3264
Bertini DP	default	34.42	34.57	23.85	45.34		3189	3094	3216
Bertini AP	default	130.53	126.61	78.70	251.01		3261	3256	3264
Bertini AP	1e-8	688.50	691.59	368.36	1125.72		3264	3261	3264

Table 4. Table 4. Runtime of HC.jl on a single core and wall clock time on 44 cores of phc -u for the katsura problem.

$n$	#sols	#real	#imag	HC.jl time (seconds)		phc -u time (seconds)
12	4,096	582	3,514	2.85	3s	7.925E+01	1m 19s
13	8,192	900	7,292	7.22	8s	2.081E+02	3m 28s
14	16,384	1,606	14,778	26.96	19s	5.065E+02	8m 27s
15	32,768	2,542	30,226	42.87	45s	1.456E+03	24m 16s
16	65,536	4,440	61,096	113.21	1m 59s	4.156E+03	1h 09m 16s
17	131,072	7,116	123,956	231.49	4m 48s	1.001E+04	2h 46m 50s
18	262,144	12,458	249,686	855.90	14m 16s	2.308E+04	6h 24m 15s
19	524,288	20,210	504,078	1330.31	22m 10s	5.696E+04	15h 49m 20s
20	1,048,576	35,206	1,013,370	3772.08	1h 02m 52s	1.317E+05	36h 34m 11s

Equations

$$H(x(t),t)=0\text{ for all }t\in [0,1)\text{ and }x(0)=x_0.$$

$$H(x(t_k),t_k)=0,k=0,1,2,\dots$$

$$\Vert x^{(j+1)}-x^{(j)}\Vert \le a^{2^j-1}\Vert x^{(1)}-x^{(0)}\Vert$$

$$\begin{array}[]{rl}J_{F}(x^{(j)})\Delta x^{(j)}&=F(x^{(j)})\\ x^{(j+1)}&=x^{(j)}-\Delta x^{(j)}\\ \end{array},\quad j=0,1,2,\ldots$$

$$\Vert J_F(x^{(0)}{)}^{-1}(J_F(x)-J_F(y))\Vert \le \omega \Vert x-y\Vert ,$$

$$\Vert \Delta x^{(0)}\Vert =\Vert J_F(x^{(0)}{)}^{-1}F(x^{(0)})\Vert \le \beta$$

$$\Vert \Delta x^{(j)}\Vert \le \frac{1}{2}\omega \frac{\sqrt{1-2h_j}}{\sqrt{1-2h_0}}\Vert \Delta x^{(j-1)}{\Vert }^2,j=1,2,3,\dots$$

$$\Vert x^{(j)}-x^{\ast }\Vert \le \frac{2\Vert \Delta x^{(j)}\Vert }{1+\sqrt{1-2\omega \frac{\sqrt{1-2h_j}}{\sqrt{1-2h_0}}\Vert \Delta x^{(j)}\Vert }},j=0,1,2,\dots$$

$$\Vert \Delta x^{(j)}\Vert \le {\left(\frac{1}{2}\right)}^{2^j-1}\Vert \Delta x^{(0)}\Vert .$$

$$\gamma (F,x)$$

$$h_0=\omega \Vert \Delta x^{(0)}\Vert \le 2(\sqrt{4a^4+a^2}-2a^2)=:h(a)$$

$${\Theta }_j:=\frac{\Vert \Delta x^{(j+1)}\Vert }{\Vert \Delta x^{(j)}\Vert }$$

$$\Vert \Delta x^{(j)}\Vert$$

$${\Theta }_j$$

$$=\frac{1}{2}\omega \Vert \Delta x^{(0)}\Vert \sqrt{\frac{1}{1-2h_0}}\prod _{\ell =0}^{j-1}{\Theta }_j=\frac{h_0}{2\sqrt{1-2h_0}}\prod _{\ell =0}^{j-1}{\Theta }_{\ell }.$$

$${\Theta }_j\le a\prod _{\ell =0}^{j-1}{\Theta }_{\ell }\le a{a}^{{\sum }_{\ell =0}^{j-1}{2}^{\ell }}=a^{2^j}.$$

$$\frac{\Vert \Delta x^{(j)}\Vert }{\Vert \Delta x^{(0)}\Vert }=\prod _{\ell =0}^{j-1}{\Theta }_{\ell }\le a^{2^j-1}.$$

$$[\omega ]=2\frac{\Vert \Delta x^{(1)}\Vert }{\Vert \Delta x^{(0)}{\Vert }^2}.$$

$$\Vert \Delta x^{(1)}\Vert \le \frac{1}{2}\omega \sqrt{\frac{1-2h_1}{1-2h_0}}\Vert \Delta x^{(0)}{\Vert }^2$$

$$\Vert \Delta x^{(1)}\Vert \le \frac{1}{2}\omega \frac{1}{1-\omega \Vert \Delta x^{(0)}\Vert }\Vert \Delta x^{(0)}{\Vert }^2$$

$$\frac{2\Vert \Delta x^{(1)}\Vert }{\Vert \Delta x^{(0)}{\Vert }^2+2\Vert \Delta x^{(0)}\Vert \Vert \Delta x^{(1)}\Vert }\le \omega .$$

$$\frac{2\Vert \Delta x^{(1)}\Vert }{\Vert \Delta x^{(0)}{\Vert }^2+2\Vert \Delta x^{(0)}\Vert \Vert \Delta x^{(1)}\Vert }\le 2\frac{\Vert \Delta x^{(1)}\Vert }{\Vert \Delta x^{(0)}{\Vert }^2}=[\omega ].$$

$$fl(x\mathrm{op}y)=(x\mathrm{op}y)(1+\delta ),|\delta |\le u,\mathrm{op}=+,-,\ast ,/$$

$$x^{(j+1)}=x^{(j)}-(J_F(x^{(j)})+E_j{)}^{-1}(F(x^{(j)})+e_j)+{\epsilon }_j$$

$$\Vert e_j\Vert \le u\Vert F(x^{(j)})\Vert +\psi (F,x^{(j)},u,\bar{u}).$$

$$\Vert E_j\Vert \le u\varphi (F,x^{(j)},n,u)$$

$$u\Vert J_F(x^{(j)}{)}^{-1}\Vert \varphi (F,x^{(j)},n,u)\le \frac{1}{8}.$$

$$\frac{\Vert x^{(j+1)}-x^{\ast }\Vert }{\Vert x^{\ast }\Vert }\approx \frac{\Vert J_F(x^{\ast }{)}^{-1}\Vert }{\Vert x^{\ast }\Vert }\psi (F,x^{\ast },u,\bar{u})+u=:\mu (x^{\ast },u,\bar{u}).$$

$$\Vert \Delta x^{(j)}\Vert \le \frac{\omega \Vert \Delta x^{(j-1)}{\Vert }^2}{2\sqrt{1-2h_0}}\le \mu (x^{\ast },u,\bar{u})\Vert x^{\ast }\Vert$$

$${\Theta }_j=\frac{\Vert \Delta x^{(j+1)}\Vert }{\Vert \Delta x^{(j)}\Vert }\le a^{2^j}j=0,1,2,\dots .$$

$$\omega \mu (x^{\ast },u,u)>a^{2^k-1}h(a)\Vert x^{\ast }\Vert$$

$${\Theta }_j=\frac{\Vert \Delta x^{(j+1)}\Vert }{\Vert \Delta x^{(j)}\Vert }\le a^{2^j}$$

Full text

1. Introduction

Systems of polynomial equations arise in many applications across the sciences including computer vision [HS97, SSN05], chemistry [Mor87], kinematics [WS11] and biology [NTI16]. A numerical method for finding all isolated solutions of a system $F$ of $n$ polynomials in $n$ unknowns is homotopy continuation [SW05]. This method constructs a homotopy $H(x,t):\mathbb{C}^{n}\times\mathbb{C}\rightarrow\mathbb{C}^{n}$ such that $H(x,t)$ is a polynomial system for all $t$, $H(x,1)=F(x)$ and $H(x,0)$ is a system whose isolated solutions are known. There is a well-developed theory [SW05] on how to construct such homotopies to guarantee, with probability one, that every isolated solution of $F(x)=0$ is the endpoint of at least one smooth solution path $x(t)$. A solution path $x(t)$ is implicitly defined by the conditions

[TABLE]

In order to trace a path $x(t)$, the problem (1.1) is treated as a sequence of problems

[TABLE]

with an (a-priori unknown) subdivision $0=t_{0}<t_{1}<\ldots<t_{M}=1$ of the interval $[0,1]$. Each of the problems (1.2) is then solved by a correction method, usually Newton’s method, under the assumption that a prediction method, e.g., Euler’s method, provides a good starting point. The choice of step size $\Delta t_{k}=t_{k+1}-t_{k}$ is usually given by an adaptive step size control. The step size must be chosen appropriately: if the step size is too large, the prediction can be outside the zone of convergence of the corrector, while a too small step size means progress is slow. There have been many efforts to design such adaptive step size controls [SC87, KX94, GS04].

In the context of polynomial homotopy continuation methods, two phenomena need particular attention. Polynomial systems often have singular solutions, and thus, the paths leading to these solutions are necessarily ill-conditioned at the end. While endgame methods [MSW90, MSW92] exist to compute singular solutions, these still require to track the solution path sufficiently close to the singularity. Usually, homotopies guarantee, with probability one, that no path passes through a singularity before reaching its endpoint. However, there is a non-negligible chance that a near-singular condition is encountered during the tracking.

Also, if two different solution paths are near to each other, then this can cause path jumping. That is, the solution that is tracked ‘jumps’ from one path to another. The typical reason is that starting from a point on the tracked path the prediction method returns a point which, according to the correction method, is a numerical approximation of a point that is on a different path than the tracked one. A possible result of path jumping is that not all isolated solutions of a polynomial system are computed.

Therefore, path tracking algorithms are required to reduce the risk of path jumping and they need to be able to handle ill-conditioned situations during the tracking. Existing software packages, e.g., Bertini [BHSW], use a version of the following path tracking algorithm. The algorithm has the following parameters: An initial step size $\Delta t>0$, a number of corrector iterations $N\geq 1$ allowed per step, a step adjustment factor $\lambda\in(0,1)$, a step expansion integer $E\geq 1$, and a minimum step size $t_{\min}$. Additionally, there is a tracking tolerance $\varepsilon>0$. This means that for a given $t$ an approximate solution $x\approx x(t)$ has to satisfy a normwise absolute error $\|x-x(t)\|_{\infty}\leq\varepsilon$.

Given an approximate solution $x\approx x(t)$ the prediction method provides an initial guess $\widehat{x}(\Delta t)\approx x(t+\Delta t)$. Then, Newton’s method iteratively improves the approximation $\widehat{x}(\Delta t)$. If the required tracking tolerance $\varepsilon$ is achieved with a most $N$ iterations, then the solution is updated and $t=t+\Delta t$. If there are $E$ successes in a row, then the step size is expanded to $\Delta t=\lambda^{-1}\Delta t$. If on the other hand the tolerance is not achieved with at most $N$ iterations, then the step size is reduced to $\Delta t=\lambda\Delta t$. If $\Delta t<t_{\min}$, the algorithm terminates with a failure. Otherwise the procedure is repeated until $t=1$ is reached.

The key to avoid path jumping is to allow only a small number of Newton iterations, typically only $N=2$ or $N=3$. In practice, this is often sufficient for the initial guess $\widehat{x}(\Delta t)$ to stay within a small enough region surrounding the path such that no path jumping occurs. However, if two paths are closer than the required tracking tolerance $\varepsilon$ for some $t^{*}\in(0,1)$, then this algorithm tends to fail for these paths. This is shown in the computational experiments in Section 5 for two different examples. Therefore, it is necessary to choose the path tracking tolerance $\varepsilon$ smaller than the minimal pairwise distance of any two paths. However, knowing the optimal choice of $\varepsilon$ a-priori is impossible. Thus, one has to use either a pessimistic value for $\varepsilon$ or resort to trial and error. But choosing $\varepsilon$ small does not only slow down the tracking of all paths, it also can result in new tracking failures. The reason for this is that Newton’s method in floating point arithmetic cannot always produce solutions whose relative normwise error is smaller than $\varepsilon$. This was shown by Tisseur in [Tis01] and is explained in detail in Section 2.2.

To avoid tracking failures due to insufficient precision Bates, Hauenstein, Sommese, and Wampler [BHSW08] developed an adaptive precision version of the above described path tracking algorithm. During the tracking, the algorithm dynamically changes the working precision such that Newton’s method can theoretically always produce solutions accurate enough for the desired tracking tolerance. This eliminates the problem of insufficient precision in exchange of a possibly high computational cost. But it also still leaves open the problem of picking a suitable tolerance $\varepsilon$.

This article introduces a new robust path tracking algorithm, Algorithm 4.3, which does not require the choice of a path tracking tolerance $\varepsilon$ or a maximal number $N$ of corrector iterations allowed per step. The key idea is to use a more intrinsic measure for accepting an initial guess in the Newton corrector: An initial guess $\widehat{x}(\Delta t)$ should only be accepted if the Newton iterates $\widehat{x}(\Delta t)=x^{(0)},x^{(1)},x^{(2)},\ldots$ satisfy

[TABLE]

for $j=1,2,\ldots$ and some fixed constant $a\in(0,\frac{1}{2}]$. If the initial guess satisfies (1.3), then it is an approximate zero. This notion was introduced by Smale [Sma86] for $a=\frac{1}{2}$ and plays an important role in the complexity analysis of polynomial homotopy continuation methods [BC11, Lai17]. Based on this idea, this article develops a new Newton corrector algorithm, Algorithm 2.7. The algorithm rejects an initial guess if (1.3) is not satisfied for some $j=1,\ldots,m$ where $m$ is dynamically chosen as the maximal number of iterations for which (1.3) can be satisfied in fixed precision floating point arithmetic.

The proposed path tracking algorithm combines the new Newton corrector algorithm with an adaptive step size control that chooses $\Delta t$ based on local geometric informations. The step size control extends an adaptive step size control developed by Deuflhard [Deu79] and combines it with the insight of [TVBV19] to use Padé approximants as prediction methods. In particular, the algorithm builds a local understanding of the region of convergence of Newton’s method and the distance to the closest singularity. This makes the developed path tracking algorithm robust against path jumping. To handle numerically challenging situations the algorithm uses mixed precision arithmetic. That is, while the bulk of the computations is performed in double precision some computations are performed, if necessary, in extended precision. This article is accompanied by a prototype implementation of the algorithm and an implementation will also be available in version 2.0 of HomotopyContinuation.jl [BT18].

This article is organized as follows. Section 2.1 reviews a Kantorovich style convergence theory of Newton’s method and Section 2.2 develops a new Newton corrector algorithm, Algorithm 2.7, based on requirement (1.3). In Section 3 the use of Padé approximants as prediction method is developed. In Section 4 the results from the previous sections are used to develop an adaptive step size control. Finally, the new path tracking algorithm, Algorithm 4.3, is stated. The algorithm’s effectiveness and robustness is shown through several numerical experiments in Section 5.

2. Newton’s Method: Theory and Computational Aspects

The path tracking algorithm for a path $x(t)$ consists of three main components: An adaptive step size routine that provides a step size $\Delta t$, a predictor that produces an initial guess $\widehat{x}$ of $x(t+\Delta t)$, and a corrector that takes $\widehat{x}$ and returns either an approximation of $x(t+\Delta t)$ or rejects $\widehat{x}$. This section focuses on Newton’s method as a corrector. The goal is to understand the size of the region of convergence of Newton’s method as well as the behavior of Newton’s method in floating point arithmetic and to translate this into a Newton corrector algorithm.

2.1. Convergence results

Let $D\subseteq\mathbb{C}^{n}$ an open set. Let $F:D\subseteq\mathbb{C}^{n}\rightarrow\mathbb{C}^{n}$ be an analytic function and $J_{F}:\mathbb{C}^{n}\rightarrow\mathbb{C}^{n}$ its Jacobian. Consider the Newton iteration

[TABLE]

starting at the initial guess $x^{(0)}\in D$. In 1948 Kantorovich [Kan48] already showed sufficient conditions for the convergence of Newton’s method and the existence of solutions. He also showed the uniqueness region of solutions and provided error estimates. A particular property of Newton’s method is that the iterates (2.1) are invariant under general linear transformations of $F$. That is, given a start value $x^{(0)}\in D$ and $A\in{\mathrm{GL}}_{n}(\mathbb{C})$ the Newton iterates of $AF(x)$ and $F(x)$ coincide. This property is referred to as affine covariance [Deu11]. In the following an affine covariant version of a Kantorovich style convergence theorem for Newton’s method is given. The statement is due to Deuflhard and Heindl [DH79] with error bounds from Yamamoto [Yam85].

Theorem 2.1** (Newton-Kantorovich [DH79, Yam85]).**

Let $F:D\subseteq\mathbb{C}^{n}\rightarrow\mathbb{C}^{n}$ be analytic. For some $x^{(0)}\in D$, assume that $J_{F}(x^{(0)})$ is invertible and that for all $x,y\in D$

[TABLE]

and $h_{0}:=\omega\beta\leq\frac{1}{2}$. Let $r^{*}=(1-\sqrt{1-2h_{0}})/\omega$ and $\bar{S}(x^{(0)},r^{*})=\{x\,|\,\|x-x^{(0)}\|\leq r^{*}\}\subseteq D$. Then:

1. (1)

The iterates (2.1) are well-defined, remain in $\bar{S}(x^{(0)},r^{*})$ and converge to a solution $x^{*}$ of $x^{*}$ of $F(x)=0$. 2. (2)

The solution is unique in $S(x^{(0)},r^{**})\cap D$ where $r^{**}=(1+\sqrt{1-2h_{0}})/\omega$.

Furthermore, assume $h<\frac{1}{2}$ and define the recursive sequence $h_{j}=\frac{h_{j-1}^{2}}{2(1-h_{j-1})^{2}}$. Then also the following error estimates hold.

[TABLE]

[TABLE]

A drawback of the Newton-Kantorovich theorem is that it is not possible to obtain sufficient conditions for the convergence of Newton’s method by only using data from the initial guess $x^{(0)}$. Instead, local information about the Lipschitz constant $\omega$ is required. The necessity of local information motivated Smale to develop his $\alpha$-theory [Sma86], which only requires data from the initial guess $x^{(0)}$ to compute sufficient conditions for the convergence of Newton’s method. This point of view has valuable features for the theory of computation. In particular, it is the building block for the complexity analysis of polynomial homotopy continuation methods [BP09, BC11, BL13, Lai17], certified path tracking algorithms [BL13] and the a posteriori certification of zeros [HS12].

In order to talk about initial guesses which are super-convergent starting with the first iteration, Smale [Sma86] introduced the concept of an approximate zero.

Definition 2.2** (Approximate zero).**

The point $x^{(0)}\in\mathbb{C}^{n}$ is an approximate zero of $F$ if the Newton iterate $x^{(j)}$ is defined for $j=1,2,\ldots$ and satisfies

[TABLE]

If $x^{(0)}$ is an approximate zero, then the true zero $x^{*}\in\mathbb{C}^{n}$ of $F$ to which the iterates are converging is the associated zero of $x^{(0)}$. ∎

Smale’s $\alpha$-theorem gives a sufficient condition for $x^{(0)}$ to be an approximate zero. The theorem uses

[TABLE]

where $D^{k}F$ is the tensor of order-$k$ derivatives of $F$ and the tensor $J_{F}^{-1}D^{k}F$ is understood as a multilinear map $A:(\mathbb{C}^{n})^{k}\to\mathbb{C}^{n}$ with norm $\|A\|:=\max_{\|v\|=1}\|A(v,\ldots,v)\|$.

Theorem 2.3** (Smale’s $\alpha$-theorem [Sma86]).**

There is a naturally defined number $\alpha_{0}$ approximately equal to 0.130707 such that if $\alpha(F,x^{(0)}):=\beta(F,x^{(0)})\,\gamma(F,x^{(0)})<\alpha_{0}$, then $x^{(0)}$ is an approximate zero of $F$.

It is also possible to give sufficient conditions for $x^{(0)}$ to be an approximate zero under the assumptions of the Newton-Kantorovich Theorem 2.1.

Lemma 2.4**.**

Using notation from Theorem 2.1 assume

[TABLE]

for a parameter $0<a<1$. Then, the contraction factors

[TABLE]

are satisfied. In particular, $x^{(0)}$ is an approximate zero if $h_{0}\leq\sqrt{2}-1$.

Proof.

From the error estimate (2.2) follows

[TABLE]

From $h_{0}\leq 2(\sqrt{4a^{4}+a^{2}}-2a^{2})<\frac{1}{2}$ follows $\frac{h_{0}}{2\sqrt{{1-2h_{0}}}}\leq a$ and therefore

[TABLE]

Statement (2.5) follows from (2.4) by observing

[TABLE]

∎

The Newton-Kantorovich theorem and Smale’s $\alpha$-theorem both give sufficient conditions for an initial guess to be an approximate zero. For the Newton-Kantorovich theorem, a (local) estimate of the Lipschitz constant $\omega$ needs to be obtained and for Smale’s $\alpha$-theorem $\gamma$ needs to be computed. The path tracking algorithm developed in this paper is based on the Newton-Kantorovich theorem since a (rough) estimate of $\omega$ can be computed with almost no additional cost during the Newton iteration.

A computational estimate $[\omega]$ of $\omega$ is

[TABLE]

This can be seen as follows. Using the error estimate (2.2)

[TABLE]

together with the observation $\frac{1-2h_{1}}{1-2h_{0}}={(1-h_{0})^{-2}}$ follows

[TABLE]

which is equivalent to

[TABLE]

The computational estimate (2.7) is now obtained by upper bounding (2.8) with

[TABLE]

2.2. Computational Aspects and Floating Point Arithmetic

After establishing the theoretical foundations of Newton’s method as well as a method to obtain a computational estimate of the Lipschitz constant $\omega$, these results are now used to guide the development of a Newton corrector algorithm. For this, the behavior of Newton’s method in floating-point arithmetic has to be taken into account.

Limit Accuracy

The following assumes the standard model of floating point arithmetic [Hig02, section 2.3]

[TABLE]

where $u$ is the unit roundoff. In standard double precision arithmetic $u=2^{-53}\approx 2.2\cdot 10^{-16}$. In [Tis01], Tisseur analyzed the limit accuracy of Newton’s method in floating-point arithmetic. Let $x^{*}\in\mathbb{C}^{n}$ be a zero of $F$ with $J_{F}(x^{*})$ non-singular, and let $x^{(0)}\in\mathbb{C}^{n}$ be an approximate zero of $F$ with associated zero $x^{*}$. In floating point arithmetic, we have

[TABLE]

where

• •

$e_{j}$ is the error made when computing the residual $F(x^{(j)})$,

• •

$E_{j}$ is the error incurred in forming $J_{F}(x^{(j)})$ and solving the linear system for $\Delta x^{(j)}$,

• •

$\varepsilon_{j}$ is the error made when adding the correction to $x^{(j)}$.

Assume that $F(x^{(j)})$ is computed in the possibly extended precision $\bar{u}\leq u$ before rounding back to working precision $u$ and assume that there exists a function $\psi$ depending on $F$, $x^{(j)}$, $u$ and $\bar{u}$ such that

[TABLE]

Similarly, assume that the error $E_{j}$ satisfies

[TABLE]

for some function $\phi$ that reflects both the instability of the linear solver and the error made when forming $J_{F}(x^{(j)})$. Then the following statement holds [Tis01, Corollary 2.3].

Theorem 2.5** ([Tis01]).**

Let $x^{(0)}$ be an approximate zero with associated zero $x^{*}$, $x^{*}\neq 0$, assume that $J_{F}(x^{*})$ is non-singular, satisfies $u\,\kappa(J_{F}(x^{*}))\leq\frac{1}{8}\,$ and assume that for all $j$

[TABLE]

Then, Newton’s method in floating point arithmetic generates a sequence of iterates $x^{(j+1)}$ whose normwise relative error decreases until the first $j$ for which

[TABLE]

In the following the value $\mu(x^{*},u,\bar{u})$ is referred to as the limit accuracy. Theorem 2.5 shows that the limit accuracy is influenced by three factors: the working precision $u$, the accuracy of the evaluation of the residual (in possibly extended precision $\bar{u}$), and the conditioning of the Jacobian. The essential consequence of this is that Newton’s method cannot always produce solutions whose normwise relative error is on the order of the working precision $u$. From the error estimate (2.2) follows that if for a given $j$

[TABLE]

then the normwise relative accuracy of $x^{(j)}$ in floating point arithmetic is only of order $\mu(x^{*},u,\bar{u})$. Assume that for a given $j$ (2.10) is satisfied. Then a computational estimate $[\mu]$ of $\mu(x^{*},u,\bar{u})$ can be obtained by computing $\|\Delta x^{(j)}\|/\|x^{(j+1)}\|$.

The following lemma shows that using extended precision improves the limit accuracy.

Lemma 2.6**.**

For extended precision $\bar{u}\leq u$ it holds $\mu(x^{*},u,\bar{u})\approx\mu(x^{*},u,u)\frac{\bar{u}}{u}+u$.

Proof.

From Section 4.3.2 in [BHSW08] follows that for a system of polynomials given as a straight line program $\psi(F,x^{(j)},u,\bar{u})$ in Theorem 2.5 is a linear function in $\bar{u}$, i.e., $\psi(F,x^{(j)},u,\bar{u})=\frac{\bar{u}}{u}\psi(F,x^{(j)},u,u)$. The statement then follows from (2.9). ∎

If the working precision $u$ is standard double-precision arithmetic, then computing with extended precision can be accomplished by using double-double arithmetic. A double-double number is represented as an unevaluated sum of a leading double and a trailing double, resulting in a unit roundoff of $2^{-106}=u^{2}$. Bailey [Bai95] pioneered double-double arithmetic, and implementations are nowadays available for a wide variety of programming languages and architectures.

Assume that for a fixed parameter $a\in(0,\frac{1}{2}]$ the Newton iterates starting at the initial guess $x^{(0)}$ are required to satisfy the contraction factors

[TABLE]

If the Newton iterates are computed with precision $\bar{u}=u$ then (2.10) implies together with Lemma 2.4 that if

[TABLE]

then there does not need to exist an initial guess $x^{(0)}$ such that the first $k$ contraction factors are satisfied. Given a fixed $k$, for instance, $k=2$, this gives a criterion when to use extended precision. Similarly, if the Newton iteration is performed with extended precision, then it is possible to use only working precision again if $\omega\mu(x^{*},u,u)<a^{2^{k}-1}h(a)$.

The working precision $u$ is insufficient if the combination of the error in the evaluation of the Jacobian and the instability in the linear system solver become too large. In this case, a multi-precision path tracking algorithm as [BHSW08] is necessary. However, as demonstrated in Section 5, using only double precision arithmetic for the linear system solver is sufficient for most applications. Nevertheless, even if the precision $u$ is sufficient to achieve the limit accuracy, the analysis of Tisseur also shows that the convergence speed of Newton’s method can decrease due to a too unstable linear system solver. In this case, the theoretical convergence speed may not be achieved which in turn can lead to not satisfying the required contraction factors. To circumvent this, the Newton updates are improved using mixed precision iterative refinement [Hig97] if $\bar{u}<u$. This stabilizes the linear system solver sufficiently to achieve the theoretical convergence speed.

Stopping Criteria

Criteria for stopping the Newton iteration are now derived. Assume that $\omega$ and the limit accuracy $\mu=\mu(x^{*},u,\bar{u})$ are known. If for any $j$ the contraction factor

[TABLE]

is not satisfied then the iteration is stopped and the initial guess is rejected. The iteration is stopped successfully at step $j$ if the next update would be smaller than the limit accuracy, i.e.,

[TABLE]

An additional Newton update is computed to obtain a computational estimate of $\mu$.

Scaling

Before the full Newton corrector algorithm is stated, a final point is addressed. So far, a simple rescaling of variables can change the behavior of the algorithm since $\omega$, $\mu$ and $\|\Delta x^{(j)}\|$ are not invariant under rescaling of variables. Additionally, if $x^{*}=0$ the accuracy needs to be measured with an absolute, and not a relative, normwise error. A rescaling of variables is formally the change of variables

[TABLE]

With $x=(x_{1},\ldots,x_{n})$ and $|x_{i}|\neq 0$ the choice $d_{i}\approx|x_{i}|$ results in new coordinates $y_{i}$ of unit order. To deal with the case $|x_{i}|=0$ as well as with possible overflows in floating point arithmetic an absolute threshold value $d_{\min}>0$ of the form

[TABLE]

has to be imposed. For instance, $d_{\min}=\max(\sqrt{u}\max_{i}|x_{i}|,u)$. To not introduce rounding errors, the scaling factors $D$ should be powers of the floating-point radix $\beta$ ($\beta=2$ in the case of IEEE-754 floating point standard arithmetic). Instead of performing the change of variables explicitly in Newton’s method the size of the Newton updates can also be measured with the weighted error $\|D^{-1}\Delta x^{(j)}\|\,.$ Using the scaling factors $D$ allows the algorithm to perform independent of the initial provided variable scaling (assuming that the initial scaling is not too extreme).

The Algorithm

Finally, a new Newton corrector algorithm, Algorithm 2.7, is stated. The algorithm builds on the results developed in this section. The idea of the algorithm is to reject an initial guess $x^{(0)}$ if the Newton iterates $x^{(0)},x^{(1)},x^{(2)},\ldots$ do not satisfy

[TABLE]

for $j=1,2,\ldots,m$ and some fixed constant $a\in(0,\frac{1}{2}]$. Here, $m$ is decided dynamically based on equation (2.13). The rejection of an initial guess is performed using the slightly stricter criterion (2.12). The algorithm needs as input estimates of the limit accuracy $\mu$ and the Lipschitz constant $\omega$ and also returns updated estimates of these quantities. During the path tracking, estimates are available by using the returned estimates from the previous steps. What to do at the beginning of the tracking, if these estimates are not available, is addressed after the algorithm.

The algorithm requires estimates of the limit accuracy $\mu$ and the Lipschitz constant $\omega$. During the path tracking, the computational estimates for both of them are available by using the computed estimates of the Newton corrector from the previous step. However, this leaves open what to do for the first step. There are two possibilities. If the start solution is the solution of a previous tracking, then computational estimates of $\mu$ and $\omega$ are already available. If this is not the case, the following heuristic, Algorithm 2.8, to determine values for $[\mu]$ and $[\omega]$ proved to be helpful. The idea is to add a small perturbation to the provided start solution and to perform two Newton steps. If the perturbation is sufficiently small, then the perturbed solution still converges to the provided start solution, and an estimate of $[\omega]$ and $[\mu]$ can be obtained. As an added benefit, this provides a test to point out invalid start solutions, e.g., due to user error. For simplicity it is assumed that it is sufficient to compute the residual with precision $u$.

3. Predictors and Padé Approximants

After carefully studying the Newton corrector in the previous section, the attention now shifts to the predictor. Recall that the role of the predictor is to produce for a given step size $\Delta t$ an initial guess $x^{(0)}$ such that the corrector converges sufficiently fast. The choice of the step size $\Delta t$ will be addressed in the next section, but before it is essential to understand the influence of $\Delta t$ on the distance of the initial guess $x^{(0)}$ to the solution path.

Consider the homotopy $H(x,t):\mathbb{C}^{n}\times\mathbb{C}\rightarrow\mathbb{C}^{n}$ and a constant $\bar{t}>0$. Given a solution $s\in\mathbb{C}^{n}$ of the system $H(x,0)=0$ assume that there is a solution path $x(t):[0,\bar{t}\,]\rightarrow\mathbb{C}^{n}$ implicitly defined by the conditions

[TABLE]

Also assume that $H_{x}(x(t),t)$ is nonsingular for all $t\in[0,\bar{t}\,]$. Then $x(t)$ can be extended to a holomorphic function with $H(x(t^{*}),t^{*})=0$ for all $t^{*}$ in some nonempty open neighborhood of 0.

Without loss of generality, in the following, only the situation at $t=0$ is considered, thus $\Delta t=t$. A predictor generates a prediction path $\widehat{x}(t):[0,\bar{t}\,]\rightarrow\mathbb{C}^{n}$ with $\widehat{x}(0)=x(0)$ and they can be classified by the local order of the prediction error $\|\widehat{x}(t)-x(t)\|$.

Definition 3.1** (Local order of a predictor).**

A predictor is of local order $p$ if there exists a $\tau>0$ and a constant $\eta_{p}\geq 0$ such that for all $t\in[0,\tau]$

[TABLE]

The constant $\tau$ is the trust region of the predictor. ∎

Example 3.2**.**

The Euler predictor $\widehat{x}(t)=x(0)+t\dot{x}(0)$ is of order $p=2$ with $\tau=\bar{t}$ since

[TABLE]

∎

There are many different families of predictors known in the literature, with the most famous ones probably being (embedded) Runge-Kutta methods. In [BHS11] it is shown that for polynomial homotopy continuation higher-order Runge-Kutta methods are substantially more efficient than the Euler predictor. However, in the following another particular class of predictors is considered: Padé approximants. See [BGM96] for an exhaustive treatment of Padé approximants.

Definition 3.3** (Padé approximant).**

Let $x(t)=\sum_{\ell=0}^{\infty}c_{\ell}t^{\ell}$ be a convergent power series. The type $(L,M)$ Padé approximant is the rational function

[TABLE]

such that $x(t)$ and $[L/M]_{x}$ (considered as formal power series) satisfy

[TABLE]

∎

Remark 3.4**.**

A type $(L,M)$ Padé approximant is a predictor of order $L+M+1$.

[SC87] demonstrates the effectiveness of a $(2,1)$ Padé approximant as a predictor in path tracking algorithms. However, Padé approximants are not only efficient predictors. They have the additional property to indicate the distance to the most nearby singularity. This property was recently highlighted by Telen, Van Barel, and Verschelde in [TVBV19] where it is used to develop a path tracking algorithm that is robust against path jumping.

Since $x(t)$ is holomorphic in a nonempty open neighborhood of 0, there is a coordinatewise expansion of $x(t)$ as a convergent power series around 0. Write $x_{j}(t)=\sum_{\ell=0}^{\infty}c_{\ell}t^{\ell}$ for the Taylor expansion of the coordinate function $x_{j}(t)$ at [math]. For sufficiently large $L$ the pole of the Padé approximant $[L/1]_{x_{j}}$ indicates the distance to the nearest singularity (also if it is a branch point). This is seen as follows. A computation shows that if $c_{L}\neq 0$,

[TABLE]

Hence the pole of $[L/1]_{x_{j}}$ is $c_{L}/c_{L+1}$ (or it is $\infty$ if $c_{L+1}=0$). Fabry’s ratio theorem [Fab96] now states that if the limit $\lim_{L\rightarrow\infty}c_{L}/c_{L+1}$ exists it is a singularity of $x(t)$.

Theorem 3.5**.**

Suppose that the coefficients of the power series $\sum_{\ell=0}^{\infty}c_{\ell}t^{\ell}$ are such that the limit $\lim_{L\rightarrow\infty}c_{L}/c_{L+1}=\lambda\neq 0$ exists. Then, the series converges uniformly inside the disk $\{|t|<|\lambda|\}$ and $\lambda$ is a singular point of $x_{j}(t)=\sum_{\ell=0}^{\infty}c_{\ell}t^{\ell}$.

For a fixed $L$ the modulus $|c_{L}/c_{L+1}|$ therefore can be assumed to be an approximation of the distance to the nearest singularity of $x(t)$ and a computational estimate of the trust-region $\tau$ of the Padé approximant. For a more extensive treatment of Padé approximants in the context of homotopy continuation, see Section 3 of [TVBV19].

For the computation of a Padé approximant of type $(L,M)$ it is necessary to compute the local derivatives $x^{(\ell)}(0)$ for $\ell=1,\ldots,L+M$. For this Mackens [Mac89] observed the following useful identity.

Lemma 3.6** ([Mac89]).**

The local derivatives $x^{(\ell)}(t)$ can be computed using the formula

[TABLE]

where

[TABLE]

In [Mac89] this identity is used for the computation of $x^{(\ell)}(0)$ by numerical differentiation. A downside of numerical differentiation is that it can suffer from catastrophic cancellation resulting in useless results. Instead of using numerical differentiation the expression (3.3) can be computed efficiently and accurately by using automatic differentiation [GW08, Chapter 13]. In particular, the cost of computing $R_{\ell}$ using automatic differentiation is at most $2\ell^{2}+\mathcal{O}(\ell)$ times the cost of evaluating $H$ by a straight-line-program. The dominating factor for the accuracy of $x^{(\ell)}(t)$ is the forward error of the linear system solving. To ensure that computed derivatives are sufficiently accurate the forward error of the linear system solving should be monitored and if necessary be reduced by using mixed precision iterative refinement [Hig97].

A robust Padé approximant implementation also needs to handle the edge cases that $x_{j}^{(\ell)}(0)=0$ for some $1\leq j\leq n$ and $1\leq\ell\leq L+M$. In [GGT13] a robust algorithm is proposed for computing Padé approximants. The provided implementation uses this algorithm for the computation of the Padé approximants.

It is also possible to obtain an estimate of the local approximation error of a Padé approximant. By comparing for each coordinate function $x_{j}(t)$ the coefficient of $t^{L+M+1}$ in

[TABLE]

it follows

[TABLE]

Considering the Taylor expansion of $[L/M]_{x_{j}}$ at [math] it follows

[TABLE]

Therefore for a Padé approximant a computational estimate $[\eta_{L+M+1}]$ of $\eta_{L+M+1}$ is

[TABLE]

where $D>0$ are the same scaling factors as used for the Newton corrector in Section 2.

4. Step Size Control and Path tracking algorithm

After studying Newton’s method as a correction method in Section 2 and Padé approximants as a prediction method in Section 3, the results are now combined to derive an adaptive step size control. Afterward, the path tracking algorithm is stated.

Consider the homotopy $H(x,t):\mathbb{C}^{n}\times\mathbb{C}\rightarrow\mathbb{C}^{n}$ with a given solution $s\in\mathbb{C}^{n}$ of the system $H(x,0)=0$ and constant $\bar{t}>0$. Assume there is a solution path $x(t):[0,\bar{t}]\rightarrow\mathbb{C}^{n}$ implicitly defined by the conditions $x(0)=s$ and (3.1), and assume that $H_{x}(x(t),t)$ is non-singular for all $t\in[0,\bar{t}\,]$. Denote by $\widehat{x}(t):[0,\bar{\tau}\,]\rightarrow\mathbb{C}^{n}$ the prediction path produced by the Padé approximant. As in Section 3 only the situation at $t=0$ is considered such that $\Delta t=t$.

The goal of the step size routine is to provide a step size $t$ such that the Newton iterates $x^{(j)}$ starting at the initial guess $\widehat{x}(t)=x^{(0)}$ satisfy for $j=0,1,2,\ldots$ the contraction factors

[TABLE]

for a fixed parameter $a\in(0,1)$, for instance $a=0.2$. Recall from Lemma 2.4 that if $a\leq\frac{1}{2}$ then $\widehat{x}(t)$ is an approximate zero. Assuming knowledge about the Lipschitz constant $\omega$ in a neighborhood of the path $x(t)$ and the theoretical quantities introduced in Section 3 it is possible to give a maximal theoretical feasible step size $t_{\max}$ such that this is the case. The approach to use the theoretical quantities $\omega$, $\eta_{p}$, and $\tau$ to determine a maximal feasible step size such that Newton’s method converges was pioneered by Deuflhard in [Deu79].

Theorem 4.1** ([DH79]).**

Let $D\subseteq\mathbb{C}^{n}$ such that for all $x\in D$ and $t\in[0,\bar{t}\,]$, $\bar{t}>0$, the Jacobian $H_{x}(x,t)$ is non-singular. Assume that for each $t\in[0,\bar{t}\,]$ there exists a convex subset $D(t)\subseteq D$ with $x(t)\in D(t)$ where $x(t)$ denotes the unique solution path in $D\times[0,\bar{t}\,]$. Let $\widehat{x}(t):[0,\bar{t}\,]\rightarrow D$ denote a prediction path of order $p$ with trust-region $\tau$, i.e., with

[TABLE]

Moreover, assume for all $t\in[0,\bar{t}\,]$ the affine covariant Lipschitz condition

[TABLE]

For fixed $h\leq\frac{1}{2}$ let $t_{\max}=\min(t^{*},\tau,\bar{t})$ where

[TABLE]

and for all $t\leq t_{\max}$ let $B(t)$ denote a ball around $\widehat{x}(t)$ with radius $(1-\sqrt{1-2h})/\omega$ and assume $B(t)\subseteq D(t)$.

Then, for all step sizes $t\leq t_{\max}$ the Newton iterates starting at $\widehat{x}(t)=x^{(0)}$ are well-defined, remain in $B(t)$, converge towards ${x}(t)$ and satisfy $\|H_{x}(\widehat{x}(t),t)^{-1}H(\widehat{x}(t),t)\|\leq\frac{h}{\omega}$.

Proof.

For $h=\frac{1}{2}$ the statement is Theorem 1.3 in [Deu79]. The more general maximal step size (4.2) and the inequality $\|H_{x}(\widehat{x}(t),t)^{-1}H(\widehat{x}(t),t)\|\leq\frac{h}{\omega}$ follows from equations (1.14a) and (1.14b) in the proof of Theorem 1.3 in [Deu79]. ∎

In Theorem 4.1 there is a choice of the parameter $h\leq\frac{1}{2}$. If this is sufficiently small then the contraction factors (4.1) are satisfied. The following corollary makes this precise.

Corollary 4.2**.**

In Theorem 4.1 choose $h\leq h(a)=2(\sqrt{4a^{4}+a^{2}}-2a^{2})$ for $a\in(0,1)$. Then, the Newton iterates starting at $\widehat{x}(t)$ satisfy the contraction factors

[TABLE]

Proof.

Since $\|H_{x}(\widehat{x}(t),t)^{-1}H(\widehat{x}(t),t)\|\leq\frac{h}{\omega}$ follows $h_{0}\leq h$ and using Lemma 2.4 follows the statement. ∎

If the step size is chosen according to Theorem 4.1, then path jumping cannot happen. However, to obtain the theoretical quantities $\omega$, $\eta_{p}$ and $\tau$ is very hard. Instead the theoretical quantities are replaced by easy to obtain computational estimates $[\omega]$, $[\eta_{p}]$ and $[\tau]$. Using the computational estimates and Corollary 4.2, an estimate $[t_{\max}]$ of the maximal feasible step size $t_{\max}$ is given by

[TABLE]

where

[TABLE]

where $0<\beta_{\tau}<1$ and $\beta_{\omega}\geq 1$ are additional safety factors, for instance $\beta_{\tau}=0.75$ and $\beta_{\omega}=10$. Instead of choosing $\beta_{\omega}$ fixed it seems worthwhile to develop a more adaptive criterion for choosing $\beta_{\omega}$ in the future.

Since $[\eta_{p}]$ and $[\omega]$ are only lower bounds and $[\tau]$ is only an upper bound for the theoretical quantities it is possible that the step size $t=[t_{\max}]$ is larger than $t_{\max}$. Then it can happen that the Newton corrector algorithm 2.7 rejects the initial guess since $\Theta_{k}>a^{2^{k}}$ for some $k$. In this case a suitable step size correction formula is

[TABLE]

which is a clear reduction since $\Theta_{k}^{2^{-k}}>a$.

Before the path tracking algorithm is stated, the similarities and differences to the step size control developed in [TVBV19] are stated. In [TVBV19] the authors develop an adaptive step size control similar to (4.3) with the difference that their algorithm uses instead of $[t^{*}]$ the step size candidate $\Delta t_{1}$ computed as follows. For $\Delta t_{1}$ an estimate $\delta$ of the distance to the nearest path is computed based on a second-order Taylor expansion around $x(0)$. This involves the computation of the Hessian of $H$ and multiple singular value decompositions. Then $\Delta t_{1}=\left(\beta_{1}\delta/[\eta_{p}]\right)^{1/p}$ where $\beta_{1}$ is a safety factor to unknown region of convergence of Newton’s method, for instance $\beta_{1}=0.005$. The cost of computing $\Delta t_{1}$ is $\mathcal{O}(n^{4})$, whereas computing $[t^{*}]$ has almost no overhead. The results in Section 5 demonstrate that the approach developed in this article is much more efficient than using $\Delta t_{1}$ without substantially increasing the risk of path jumping.

Finally, the path tracking algorithm, Algorithm 4.3, is stated. It is assumed that the computational estimates $[\omega]$ and $[\mu]$ for the start solution $s$ are available. These are either available as a result of a previous path tracking or by using Algorithm 2.8.

5. Computational Experiments

In this section, numerical experiments are shown to illustrate the effectiveness of the proposed path tracking algorithm. A prototype implementation of the path tracking algorithm together with all the data necessary to run the experiments is available at

https://doi.org/10.5281/zenodo.3667414  .

The path tracking algorithm will also be implemented in version 2.0 of the Julia package HomotopyContinuation.jl [BT18].

In the experiments, the implementation is compared with the state of the art. For the different solvers the following notations are used in the experiments:

[TABLE]

All solvers are used intentionally with the default settings unless otherwise mentioned since for a non-expert user it is very hard to understand which parameters need to be changed. Note that Bertini uses by default a path tracking tolerance of $10^{-5}$ before and $10^{-6}$ during the endgame. The experiments are performed on a 24 GB RAM machine with an Intel Core i5-7500 CPU working at 3.40 GHz. All solvers use only one core for all the experiments unless stated otherwise. The experiments are designed such that the tracked solution paths are all smooth. Therefore endgame algorithms are not necessary. In all experiments the implementation of the proposed algorithm uses a $(2,1)$ Padé approximant with parameters $a=0.2$, $\beta_{\omega}=10$ and $\beta_{\tau}=0.75$.

5.1. A Family of Hyperbolas

This is an experiment from [TVBV19] which constructs a situation which causes danger for path jumping. Consider a family of hyperbolas parametrized by a real parameter $\rho$ and defined by

[TABLE]

See Figure 1 for an illustration. For general $t^{*}$ the homotopy $H_{\rho}$ has the two solutions

[TABLE]

Note that the Jacobian $H_{x}=\frac{\partial}{\partial x}H$ is singular for $t^{*}=\frac{1}{2}+\rho\sqrt{-1}$ and $x=\sqrt{(t^{*}-\frac{1}{2})^{2}+\rho^{2}}$. Consider $\rho\neq 0$ here, such that $x(t):[0,1]\rightarrow\mathbb{C}$ is a smooth path. The smaller $|\rho|$, the closer the two branch points move to the line segment $[0,1]$. Figure 1 shows that as the value of $\rho>0$ decreases, the two solution paths approach each other for parameter values $t^{*}\approx 0.5$ which causes danger for path jumping. The experiments confirm this with the results shown in Table 1. Note that HC.jl 1.4 and Bertini AP/DP are prone to path jumping. The new path tracking algorithm, on the other hand, handles this situation correctly.

5.2. Alt’s problem

Alt’s problem, formulated in 1923, is to count the number of four-bar linkages whose coupler curve interpolates nine general points in the plane. In 1992, Morgan, Sommese, and Wampler [WMS92] provided a numerical proof to Alt’s problem that there are generically 1442 non-degenerate four-bar linkages. Due to Roberts cognates and a two-fold symmetry, the resulting polynomial system generically has 8652 regular solutions. Here the formulation of the polynomial system as an affine polynomial system in 24 variables and the 16 parameters $(\delta,\widehat{\delta})\in\mathbb{C}^{8}\times\mathbb{C}^{8}$ is considered. Since the problem is formulated in isotropic coordinates, the physically meaningful configurations correspond to choices of parameters such that $\delta$ and $\widehat{\delta}$ are complex conjugates. Consider the ‘general’ situation where solutions are tracked from generic parameter values $\delta_{1}\in\mathbb{C}^{8}\times\mathbb{C}^{8}$ to generic physically meaningful parameter values $(\delta_{0},\bar{\delta_{0}})$ with $\delta_{0}\in\mathbb{C}^{8}$. The results in Table 2 show that even this general situation results in numerically challenging paths which Bertini AP cannot handle with its default settings. After decreasing the path tracking tolerance to $10^{-8}$, in half of the cases all solutions are found. The proposed algorithm, on the other hand, reliably computes all 1442 solutions without any path jumping in a fraction of the time.

Figure 2 illustrates the behavior of the proposed algorithm for one particular numerically challenging path. This path closely passes at around $t=0.93$ a singularity. As expected, this increases the Lipschitz constant $\omega$ and the limit accuracy $\mu$. The limit accuracy decreases sharply as soon as the algorithm switches to extended precision. After passing the problematic region, the algorithm quickly switches back to only using double precision arithmetic, as indicated by a sharp increase of $\mu$. The lower part of Figure 2 depicts different values associated with the Jacobian $J=H_{x}(x(t),t)$ along the path: the condition number $\kappa(J)$, the componentwise relative condition number $\mathrm{cond}(J)$ and $\mathrm{cond}(J,x)$. See [Hig02, Sec. 7] for the definitions. The values of $\kappa(J)$ and $\mathrm{cond}(J)$ obtain a maximum of $2.1\times 10^{25}$ resp. $1.6\times 10^{19}$. From these values, it seems hopeless to track the path in double precision arithmetic due to the well-known rule of thumb that one expects to lose around $\log_{10}(\kappa(J))$ digits of accuracy in the linear system solving. However, the componentwise relative forward error of the computed Newton updates is only governed by the much tamer $\mathrm{cond}(J,x)$ which is at most $5.8\times 10^{10}$. This explains why it is still possible to track the path by only using mixed precision iterative refinement. Using the proposed path tracking algorithm the path needs $253$ steps in total with only two rejected steps and a total runtime of 13 milliseconds. Trying to compute the path with Bertini AP results in a path failure after around 90 seconds and over 5000 steps due to supposedly insufficient precision of at most 1024 bits. After changing the required tracking tolerance to $10^{-12}$ Bertini AP successfully tracks the path in $120$ seconds.

5.3. Steiner’s Conic Problem

A classic problem in enumerative geometry is Steiner’s conic problem. It asks: How many plane conics are tangent to five given conics in general position? See [EH16] for historic remarks and a modern intersection theory treatment of this classic problem. Here the formulation of Steiner’s conic problem from [BST20] is used where the resulting polynomial system has 15 variables and 30 parameters. To test the path tracking algorithm consider the case of a parameter homotopy [MS89] from generic complex parameters $c\in\mathbb{C}^{30}$ to generic real parameters $r\in\mathbb{R}^{30}$. The results for 50 parameter homotopies are shown in Table 3. As for Alt’s problem, the proposed path tracking algorithm handles all instances without any failure or path jumping. However, even these generic instances pose problems for other path tracking algorithms with Bertini AP losing solutions almost always solutions using the default settings. After the path tracking tolerance is manually decreased to $10^{-8}$ in almost all instances all solutions are found.

5.4. The Katsura Family

The katsura family of systems is named after the problem posed by Katsura [Kat94], see [Kat90] for a description of its relevance to applications. The katsura-$n$ problem consists of $n$ quadratic equations and one linear equation. The number of solutions equals $2^{n}$, the product of the degrees of all polynomials in the system.

Table 4 summarizes the characteristics and run times on katsura-$n$, for $n$ ranging from 12 to 20 on a single CPU. For reference the runtime of phc -u is stated as reported in [TVBV19]. Note that these computations were performed on two 22-core 2.2 GHz processors, i.e., using 44 CPUs in total. In [TVBV19] it is also mentioned that the computation required the use of homogeneous coordinates since otherwise intermediate coordinate growth resulted in tracking failures for some paths. However, this is not necessary if a dynamic rescaling of the variables is performed as derived in Section 2.2.

6. Conclusion and Future Work

This article proposed a mixed precision path tracking algorithm for numerical path tracking in polynomial homotopy continuation. The results of the computational experiments demonstrate that the algorithm is very robust, can handle numerically challenging situations and, at the same time, is faster than existing software implementations. An implementation is available, and the algorithm will also be integrated into version 2.0 of HomotopyContinuation.jl. It is expected that the techniques in this article are also helpful in developing more efficient and robust endgame algorithms to deal with singular solutions and diverging paths.

Acknowledgement

The author wants to thank Paul Breiding, Peter Deuflhard, Michael Joswig, Simon Telen, and Marc Van Barel for helpful discussions. The author was supported by the Deutsche Forschungsgemeinschaft (German Research Foundation) Graduiertenkolleg Facets of Complexity (GRK 2434).