Mixed Precision Path Tracking for Polynomial Homotopy Continuation

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.


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) : C n × C → 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 H(x(t), t) = 0 for all t ∈ [0, 1) and x(0) = x 0 .
(1.1) 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 ∆t > 0, a number of corrector iterations N ≥ 1 allowed per step, a step adjustment factor λ ∈ (0, 1), a step expansion integer E ≥ 1, and a minimum step size t min . Additionally, there is a tracking tolerance ε > 0. This means that for a given t an approximate solution x ≈ x(t) has to satisfy a normwise absolute error x − x(t) ∞ ≤ ε.
Given an approximate solution x ≈ x(t) the prediction method provides an initial guess x(∆t) ≈ x(t + ∆t). Then, Newton's method iteratively improves the approximation x(∆t). If the required tracking tolerance ε is achieved with a most N iterations, then the solution is updated and t = t + ∆t. If there are E successes in a row, then the step size is expanded to ∆t = λ −1 ∆t. If on the other hand the tolerance is not achieved with at most N iterations, then the step size is reduced to ∆t = λ∆t. If ∆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 x(∆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 ε for some t * ∈ (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 ε smaller than the minimal pairwise distance of any two paths. However, knowing the optimal choice of ε a-priori is impossible. Thus, one has to use either a pessimistic value for ε or resort to trial and error. But choosing ε 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 ε. 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 ε.
This article introduces a new robust path tracking algorithm, Algorithm 4.3, which does not require the choice of a path tracking tolerance ε 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 x(∆t) should only be accepted if the Newton iterates for j = 1, 2, . . . and some fixed constant a ∈ (0, 1 2 ]. If the initial guess satisfies (1.3), then it is an approximate zero. This notion was introduced by Smale [Sma86] for a = 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, . . . , 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 ∆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.

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 ∆t, a predictor that produces an initial guess x of x(t + ∆t), and a corrector that takes x and returns either an approximation of x(t + ∆t) or rejects 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 ⊆ C n an open set. Let F : D ⊆ C n → C n be an analytic function and J F : C n → C n its Jacobian. Consider the Newton iteration starting at the initial guess x (0) ∈ 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) ∈ D and A ∈ GL n (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].
For some x (0) ∈ D, assume that J F (x (0) ) is invertible and that for all x, y ∈ D

(2.3)
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 ω is required. The necessity of local information motivated Smale to develop his α-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) ∈ C n is an approximate zero of F if the Newton iterate x (j) is defined for j = 1, 2, . . . and satisfies If x (0) is an approximate zero, then the true zero x * ∈ C n of F to which the iterates are converging is the associated zero of x (0) .
Smale's α-theorem gives a sufficient condition for x (0) to be an approximate zero. The theorem uses where D k F is the tensor of order-k derivatives of F and the tensor J −1 F D k F is understood as a multilinear map A : (C n ) k → C n with norm A := max v =1 A(v, . . . , v) . Theorem 2.3 (Smale's α-theorem [Sma86]). There is a naturally defined number α 0 approximately equal to 0.130707 such that if α(F, x (0) ) := β(F, x (0) ) γ(F, x (0) ) < α 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 for a parameter 0 < a < 1. Then, the contraction factors and the error bounds are satisfied. In particular, Proof. From the error estimate (2.2) follows The Newton-Kantorovich theorem and Smale's α-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 ω needs to be obtained and for Smale's α-theorem γ needs to be computed. The path tracking algorithm developed in this paper is based on the Newton-Kantorovich theorem since a (rough) estimate of ω can be computed with almost no additional cost during the Newton iteration.
A computational estimate [ω] of ω is This can be seen as follows. Using the error estimate (2.2) (2.8) The computational estimate (2.7) is now obtained by upper bounding (2.8) with

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 ω, 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 where u is the unit roundoff. In standard double precision arithmetic u = 2 −53 ≈ 2.2 · 10 −16 . In [Tis01], Tisseur analyzed the limit accuracy of Newton's method in floatingpoint arithmetic. Let x * ∈ C n be a zero of F with J F (x * ) non-singular, and let x (0) ∈ C n be an approximate zero of F with associated zero x * . In floating point arithmetic, we have • 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 ∆x (j) , • ε j is the error made when adding the correction to x (j) . Assume that F (x (j) ) is computed in the possibly extended precisionū ≤ u before rounding back to working precision u and assume that there exists a function ψ depending on F , x (j) , u andū such that Similarly, assume that the error E j satisfies for some function φ 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].
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 In the following the value µ(x * , 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ū), 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 then the normwise relative accuracy of x (j) in floating point arithmetic is only of order µ(x * , u,ū). Assume that for a given j (2.10) is satisfied. Then a computational estimate [µ] of µ(x * , u,ū) can be obtained by computing ∆x (j) / x (j+1) .
The following lemma shows that using extended precision improves the limit accuracy.
If the working precision u is standard double-precision arithmetic, then computing with extended precision can be accomplished by using double-double arithmetic. A doubledouble 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 ∈ (0, 1 2 ] the Newton iterates starting at the initial guess x (0) are required to satisfy the contraction factors If the Newton iterates are computed with precisionū = u then (2.10) implies together with Lemma 2.4 that if ωµ(x * , u, u) > a 2 k −1 h(a) x * (2.11) 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 ωµ(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ū < 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 ω and the limit accuracy µ = µ(x * , u,ū) are known. If for any j the contraction factor 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., ω ∆x (j−1) 2 An additional Newton update is computed to obtain a computational estimate of µ.
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 ω, µ and ∆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 has to be imposed. For instance, d min = max( √ u max i |x i |, u). To not introduce rounding errors, the scaling factors D should be powers of the floating-point radix β (β = 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 ∆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) , . . . do not satisfy for j = 1, 2, . . . , m and some fixed constant a ∈ (0, 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 µ and the Lipschitz constant ω 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.

Algorithm 2.7 Newton Corrector
Input: F : C n → C n , x (0) ∈ C n , a ∈ (0, 1 2 ], estimate [µ] > 0 of the limit accuracy µ, estimate [ω] > 0 of ω, evaluation precisionū ≤ u, and positive scaling factors D such that the coordinates of D −1 x (0) are of unit order. Output: Boolean indicating whether x (0) was accepted, approximationx ∈ C n of a zero x * of F , updated estimate of the (limit) accuracy µ at x * , updated estimate of ω, number of updates j and last contraction factor Θ j−2 . r ← Evaluate F at x (j) with precisionū and round result to precision u 5: then Approaching limit accuracy 13: r ← Evaluate F at x (j+1) with precisionū and round result to precision u 14: Solve J F (x (j+1) )∆x (j+1) = r 15: 16: [µ] ← D −1 ∆x (j+2) Update of the limit accuracy The algorithm requires estimates of the limit accuracy µ and the Lipschitz constant ω. 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 µ and ω are already available. If this is not the case, the following heuristic, Algorithm 2.8, to determine values for [µ] and [ω] 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 [ω] and [µ] 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.

Algorithm 2.8 Model Initialization Heuristic
Input: Candidate x (0) ∈ C n , a ∈ (0, 1), scaling factors D such that the coordinates of D −1 x (0) are of unit order. Output: Boolean indicating whether the initialization was successful, estimate [µ] of the limit accuracy µ of the associated zero of x (0) , and an estimate [ω] of ω.

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 ∆t an initial guess x (0) such that the corrector converges sufficiently fast. The choice of the step size ∆t will be addressed in the next section, but before it is essential to understand the influence of ∆t on the distance of the initial guess x (0) to the solution path.
Consider the homotopy H(x, t) : C n × C → C n and a constantt > 0. Given a solution s ∈ C n of the system H(x, 0) = 0 assume that there is a solution path x(t) : [0,t ] → C n implicitly defined by the conditions H(x(t), t) = 0 for all t ∈ [0,t ] and x(0) = s . (3.1) Also assume that H x (x(t), t) is nonsingular for all t ∈ [0,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 ∆t = t. A predictor generates a prediction path x(t) : [0,t ] → C n with x(0) = x(0) and they can be classified by the local order of the prediction error x(t) − x(t) .
Definition 3.1 (Local order of a predictor). A predictor is of local order p if there exists a τ > 0 and a constant η p ≥ 0 such that for all t ∈ [0, τ ] The constant τ is the trust region of the predictor.
Example 3.2. The Euler predictor x(t) = x(0) + tẋ(0) is of order p = 2 with τ =t since 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. a 0 + a 1 t + a 2 t 2 + · · · + a L t L 1 + b 1 t + b 2 t 2 + · · · + b M t M such that x(t) and [L/M ] x (considered as formal power series) satisfy 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) = ∞ =0 c t for the Taylor expansion of the coordinate function x j (t) at 0. 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 = 0, Hence the pole of [L/1] x j is c L /c L+1 (or it is ∞ if c L+1 = 0). Fabry's ratio theorem [Fab96] now states that if the limit lim L→∞ c L /c L+1 exists it is a singularity of x(t).
Theorem 3.5. Suppose that the coefficients of the power series ∞ =0 c t are such that the limit lim L→∞ c L /c L+1 = λ = 0 exists. Then, the series converges uniformly inside the disk {|t| < |λ|} and λ is a singular point of x j (t) = ∞ =0 c t . 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 trustregion τ 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 ( ) (0) for = 1, . . . , L + M . For this Mackens [Mac89] observed the following useful identity.
Lemma 3.6 ( [Mac89]). The local derivatives x ( ) (t) can be computed using the formula In [Mac89] this identity is used for the computation of x ( ) (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 using automatic differentiation is at most 2 2 + O( ) times the cost of evaluating H by a straight-lineprogram. The dominating factor for the accuracy of x ( ) (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 (0) = 0 for some 1 ≤ j ≤ n and 1 ≤ ≤ 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 where D > 0 are the same scaling factors as used for the Newton corrector in Section 2.

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) : C n × C → C n with a given solution s ∈ C n of the system H(x, 0) = 0 and constantt > 0. Assume there is a solution path x(t) : [0,t] → 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 ∈ [0,t ]. Denote by x(t) : [0,τ ] → C n the prediction path produced by the Padé approximant. As in Section 3 only the situation at t = 0 is considered such that ∆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 x(t) = x (0) satisfy for j = 0, 1, 2, . . . the contraction factors for a fixed parameter a ∈ (0, 1), for instance a = 0.2. Recall from Lemma 2.4 that if a ≤ 1 2 then x(t) is an approximate zero. Assuming knowledge about the Lipschitz constant ω 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 ω, η p , and τ to determine a maximal feasible step size such that Newton's method converges was pioneered by Deuflhard in [Deu79].
Moreover, assume for all t ∈ [0,t ] the affine covariant Lipschitz condition For fixed h ≤ 1 2 let t max = min(t * , τ,t) where and for all t ≤ t max let B(t) denote a ball around x(t) with radius (1 − √ 1 − 2h)/ω and assume B(t) ⊆ D(t).
Then, for all step sizes t ≤ t max the Newton iterates starting at x(t) = x (0) are welldefined, remain in B(t), converge towards x(t) and satisfy H x ( x(t), t) −1 H( x(t), t) ≤ h ω . Proof. For h = 1 2 the statement is Theorem 1.3 in [Deu79]. The more general maximal step size (4.2) and the inequality H x ( x(t), t) −1 H( x(t), t) ≤ h ω 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 ≤ 1 2 . If this is sufficiently small then the contraction factors (4.1) are satisfied. The following corollary makes this precise. ∈ (0, 1). Then, the Newton iterates starting at x(t) satisfy the contraction factors ω follows h 0 ≤ 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 ω, η p and τ is very hard. Instead the theoretical quantities are replaced by easy to obtain computational estimates [ω], [η p ] and [τ ]. Using the computational estimates and Corollary 4.2, an estimate [t max ] of the maximal feasible step size t max is given by where 0 < β τ < 1 and β ω ≥ 1 are additional safety factors, for instance β τ = 0.75 and β ω = 10. Instead of choosing β ω fixed it seems worthwhile to develop a more adaptive criterion for choosing β ω in the future. Since [η p ] and [ω] are only lower bounds and [τ ] 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 Θ k > a 2 k for some k. In this case a suitable step size correction formula is which is a clear reduction since Θ 2 −k 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 ∆t 1 computed as follows. For ∆t 1 an estimate δ 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 ∆t 1 = (β 1 δ/[η p ]) 1/p where β 1 is a safety factor to unknown region of convergence of Newton's method, for instance β 1 = 0.005. The cost of computing ∆t 1 is 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 ∆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 [ω] and [µ] 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. (t, ∆t) ← (0, ∞) 3: x ← s 4:ū ← u

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:

HC.jl
The (5.1) See Figure 5 for an illustration. For general t * the homotopy H ρ has the two solutions Note that the Jacobian H x = ∂ ∂x H is singular for t * = 1 2 +ρ √ −1 and x = (t * − 1 2 ) 2 + ρ 2 . Consider ρ = 0 here, such that x(t) : [0, 1] → C is a smooth path. The smaller |ρ|, the closer the two branch points move to the line segment [0, 1]. Figure 5 shows that as the value of ρ > 0 decreases, the two solution paths approach each other for parameter values t * ≈ 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 fourbar 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 (δ, δ) ∈ C 8 × C 8 is considered. Since the problem is formulated in isotropic coordinates, the physically meaningful configurations correspond to choices of parameters such that δ and δ are complex conjugates. Consider the 'general' situation where solutions are tracked from generic parameter values δ 1 ∈ C 8 ×C 8 to generic physically meaningful parameter values (δ 0 ,δ 0 ) with δ 0 ∈ 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 ω and the limit accuracy µ. The limit accuracy decreases sharply as soon as the algorithm switches to extended precision. After passing  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.
the problematic region, the algorithm quickly switches back to only using double precision arithmetic, as indicated by a sharp increase of µ. 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 κ(J), the componentwise relative condition number cond(J) and cond(J, x). See [Hig02, Sec. 7] for the definitions. The values of κ(J) and cond(J) obtain a maximum of 2.1×10 25 resp. 1.6×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 (κ(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 cond(J, x) which is at most 5.8 × 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. Figure 2. Behavior of the path tracking algorithm along a single challenging path.

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 ∈ C 30 to generic real parameters r ∈ 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.  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.
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. Runtime of HC.jl on a single core and wall clock time on 44 cores of phc -u for the katsura problem. 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.

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.