# Solving polynomial systems via homotopy continuation and monodromy

Solving polynomial systems via homotopy continuation and monodromy Abstract We study methods for finding the solution set of a generic system in a family of polynomial systems with parametric coefficients. We present a framework for describing monodromy-based solvers in terms of decorated graphs. Under the theoretical . that monodromy actions are generated uniformly, we show that the expected number of homotopy paths tracked by an algorithm following this framework is linear in the number of solutions. We demonstrate that our software implementation is competitive with the existing state-of-the-art methods implemented in other software packages. 1. Introduction Homotopy continuation has become a standard technique to find approximations of solutions of polynomial systems. There is an early popular text on the subject and its applications in the study by Morgan (1987). This technique is the backbone of Numerical Algebraic Geometry, the area that classically addresses the questions of complex algebraic geometry through algorithms that employ numerical approximate computation. The chapter in the study by Sommese et al. (2005, Section 8) is the earliest introduction and the book by Sommese & Wampler (2005) is the primary reference in the area. Families of polynomial systems with parametric coefficients play one of the central roles. Most homotopy continuation techniques could be viewed as going from a generic system in the family to a particular one. This process is commonly referred to as degeneration. Going in the reverse direction it may be called deformation, undegeneration or regeneration depending on the literature. Knowing the solutions of a generic system, one can use coefficient-parameter homotopy (Sommese & Wampler, 2005, Section 7) to get to the solution of a particular one. The main problem that we address here is how to solve a generic system in a family of systems, $$F_{p} = \left(f_{p}^{(1)},\ldots,\,f_{p}^{(N)}\right) = 0, \quad f_{p}^{(i)}\in{\mathbb{C}}[p][x],\ i=1,\ldots,N,$$ with finitely many parameters p and n variables x. In the main body of the paper we restrict our attention to linear parametric families of systems, defined as systems with affine linear parametric coefficients, such that for a generic p we have a nonempty finite set of solutions x to $$F_{p}(x)=0$$. This implies N ≥ n. The number of parameters is arbitrary, but we require that for a generic x there exists p with $$F_{p}(x)=0$$. These restrictions are made for the sake of simplicity. We explain what modifications are needed to apply our approach in more general settings in Section 7. Linear parametric systems form a large class that includes sparse polynomial systems. These are square (n = N) systems with a fixed monomial support for each equation and a distinct parameter for the coefficient of each monomial. Polyhedral homotopy methods for solving sparse systems stem from the BKK (Bernstein, Khovanskii, Kouchnirenko) bound on the number of solutions (Bernstein, 1975); the early work on algorithm development was done in the studies by Verschelde et al. (1994) and Huber & Sturmfels (1995). Polyhedral homotopies provide an optimal solution to sparse systems in the sense that they are designed to follow exactly as many paths as the number of solutions of a generic system (the BKK bound). The method that we propose is clearly not optimal in the above sense. The expected number of homotopy paths followed can be larger than the number of solutions, though not significantly larger. We also use linear segment homotopies that are significantly simpler and less expensive to follow in practice. Our current implementation shows it is competitive with the state-of-the-art implementations of polyhedral homotopies in PHCpack (Verschelde, 1999) and HOM4PS2 (Lee et al., 2008) for solving sparse systems. In a setting more general than sparse, we demonstrate examples of linear parametric systems for which our implementation exceeds the capabilities of the existing sparse system solvers and blackbox solvers based on other ideas. The idea of using the monodromy action induced by the fundamental group of the regular locus of the parameter space has been successfully employed throughout Numerical Algebraic Geometry. One of the main tools in the area, numerical irreducible decomposition, can be efficiently implemented using the monodromy breakup algorithm, which first appeared in the study by Sommese et al. (2001). One parallel incarnation of the monodromy breakup algorithm is described in the study by Leykin & Verschelde (2009). In fact, the main idea in that work is close in spirit to what we propose in this article. The idea to use monodromy to find solutions drives numerical implicitization (Chen & Kileel, 2016) and appears in other works such as del Campo & Rodriguez (2017). Computing monodromy groups numerically, as in the studies by Leykin & Sottile (2009) and Hauenstein et al. (2017), requires more computation than just finding solutions. One can approach this computation with the same methodology as we propose; see 5 of Section 7. Our main contribution is a new framework to describe algorithms for solving polynomial systems using monodromy; we call it the Monodromy Solver (MS) framework. We analyze the complexity of our main algorithm both theoretically assuming a certain statistical model and experimentally on families of examples. The analysis gives us grounds to say that the expected number of paths tracked by our method is linear, with a small coefficient, as the number of solutions grows. Our method and its implementation not only provide a new general tool for solving polynomial systems, but also can solve some problems out of reach for other existing software. The structure of the paper is as follows. We give a brief overview of the MS method intermingled with some necessary preliminaries in Section 2. An algorithm following the MS framework depends on a choice of strategy with several possibilities outlined in Section 3. Statistical analysis of the method is the topic of Section 4. The implementation is discussed in Section 5 together with the side topic of certification of the solution set. The results of our experiments on selected example families highlighting various practical computational aspects are in Section 6. The reader may also want look at examples of systems in Sections 6.1 and 6.2 before reading some earlier sections. Possible generalizations of the MS technique and the future directions to explore are presented in Section 7. 2. Background and framework overview Let $$m, n\in {\mathbb{N}}$$. We consider the complex linear space of square systems $$F_{p}$$, $$p\in {\mathbb{C}}^{m}$$, where the monomial support of $$f_{p}^{(1)},\dots ,\,f_{p}^{(n)}$$ in the variables $$x=(x_{1},\dots ,x_{n})$$ is fixed and the coefficients vary. By a base space B we mean a parametrized linear variety of systems. We think of it as the image of an affine linear map $$\varphi : p\mapsto F_{p}$$ from a parameter space $${\mathbb{C}}^{m}$$ with coordinates $$p=(p_{1},\ldots ,p_{m})$$ to the space of systems. We assume the structure of our family is such that the projection $$\pi$$ from the solution variety $$V = \left\{\left(F_{p},x\right) \in B\times{\mathbb{C}}^{n} \mid F_{p}(x)=0\right\}$$ to B gives us a branched covering, i.e., the fibre $$\pi ^{-1}\left (F_{p}\right )$$ is finite of the same cardinality for a generic p. The discriminant variety D in this context is the subset of the systems in the base space with nongeneric fibres; it is also known as the branch locus of $$\pi$$. The fundamental group$$\pi _{1}(B\setminus D)$$—note that $$\pi _{1}$$ is a usual topological notation that is not related to the map $$\pi$$ above—as a set consists of loops, i.e., paths in B ∖ D starting and finishing at a fixed p ∈ B ∖ D considered up to homotopy equivalence. The definition, more details to which one can find in Section 2.1, does not depend on the point p, since B ∖ D is connected. Each loop induces a permutation of the fibre $$\pi ^{-1}\left (F_{p}\right )$$, which is referred to as a monodromy action. Our goal is to find the fibre of one generic system in our family. Our method is to find one pair $$(p_{0},x_{0})\in V$$ and use the monodromy action on the fibre $$\pi ^{-1}\left (F_{p_{0}}\right )$$ to find its points. We assume that this action is transitive, which is the case if and only if the solution variety V is irreducible. If V happens to be reducible we replace V with its unique dominant irreducible component as explained in Remark 2.2. 2.1 Monodromy We briefly review the basic facts concerning monodromy groups of branched coverings. With notation as before fix a system $$F_{p}\in B \setminus D$$ and consider a loop $$\tau$$ without branch points based at $$F_{p}$$; that is, a continuous path $$\tau : [0,1] \rightarrow B \setminus D$$ such that $$\tau (0) = \tau (1) = F_{p}.$$ Suppose we are also given a point $$x_{i}$$ in the fibre $$\pi ^{-1} \left (F_{p}\right )$$ with d points $$x_{1}, x_{2}, \ldots , x_{d}.$$ Since $$\pi$$ is a covering map the pair $$(\tau , x_{i})$$ corresponds to a unique lifting$$\widetilde{\tau _{i}},$$ a path $$\widetilde{\tau_{i}} : [0,1] \rightarrow V$$ such that $$\widetilde{\tau _{i}} (0) = x_{i}$$ and $$\widetilde{\tau _{i}} (1) = x_{j}$$ for some 1 ≤ j ≤ d. Note that the reversal of $$\tau$$ and $$x_{j}$$ lifts to a reversal of $$\widetilde{\tau _{i}}$$. Thus, the loop $$\tau$$ induces a permutation of the set $$\pi ^{-1} \left (F_{p}\right ).$$ We have a group homomorphism $$\varphi : \pi_{1} \left(B\setminus D, F_{p}\right) \rightarrow{\mathrm{S}}_{d}$$ whose domain is the usual fundamental group of B ∖ D based at $$F_{p}$$. The image of $$\varphi$$ is the monodromy group associated to $$\pi ^{-1} \left (F_{p}\right ).$$ The monodromy group acts on the fibre $$\pi ^{-1} \left (F_{p}\right )$$ by permuting the solutions of $$F_{p}$$. Remark 2.1 A reader familiar with the notion of a monodromy loop in the discussion of Sommese & Wampler (2005, Section 15.4) may think of this keyword referring to a representative of an element of the fundamental group together with its liftings to the solution variety and the induced action on the fibre. For the purposes of this article we need to be clear about the ingredients bundled in this term. We have not used any algebraic properties so far. The construction of the monodromy group above holds for an arbitrary covering with finitely many sheets. The monodromy group is a transitive subgroup of $${\mathrm{S}}_{d}$$ whenever the total space is connected. In our setting since we are working over $${\mathbb{C}}$$, this occurs precisely when the solution variety is irreducible. Remark 2.2 For a linear family we can show that there is at most one irreducible component of the solution variety V for which the restriction of the projection $$\left (F_{p},x\right )\mapsto x$$ is dominant (that is, its image is dense). We call such component the dominant component. Indeed, let U be the locus of points $$\left (F_{p},x\right )\in \pi ^{-1}(B\setminus D)$$ such that the restriction of the x-projection map is locally surjective, and the solution to the linear system of equations $$F_{p}(x)=0$$ in p has the generic dimension. Being locally surjective could be interpreted either in the sense of Zariski topology or as inducing surjection on the tangent spaces. Then either U is empty or $$\overline{U}$$ is the dominant component we need, since it is a vector bundle over an irreducible variety, and is hence irreducible. In the rest of the paper when we say solution variety, we mean the dominant component of the solution variety. In particular, for sparse systems restricting the attention to the dominant component translates into looking for solutions only in the torus $$\left ( {\mathbb{C}}^{\ast }\right )^{n}$$. 2.2 Homotopy continuation Given two points $$F_{p_{1}}$$ and $$F_{p_{2}}$$ in the base space B, we may form the family of systems $$H(t) = (1-t)F_{p_{1}} + t F_{p_{2}}\,,\quad t\in[0,1],$$ known as the linear segment homotopy between the two systems. If $$p_{1}$$ and $$p_{2}$$ are sufficiently generic for each t ∈ [0, 1], we have H(t) outside the real codimension 2 set D. Consequently, each system H(t) has a finite and equal number of solutions. This homotopy is a path in B; a lifting of this path in the solution variety V is called a homotopy path. The homotopy paths of H(t) establish a one-to-one correspondence between the fibres $$\pi ^{-1}\left (F_{p_{1}}\right )$$ and $$\pi ^{-1}\left (F_{p_{2}}\right )$$. Remark 2.3 Note that $$\gamma F_{p}$$ for $$\gamma \in {\mathbb{C}}\setminus \{0\}$$ has the same solutions as $$F_{p}$$. Let us scale both ends of the homotopy by taking a homotopy between $$\gamma _{1}F_{p_{1}}$$ and $$\gamma _{2}F_{p_{2}}$$ for generic $$\gamma _{1}$$ and $$\gamma _{2}$$. If the coefficients of $$F_{p}$$ are homogeneous in p then $$H^{\prime}(t) = (1-t)\gamma_{1}F_{p_{1}} + t \gamma_{2}F_{p_{2}} = F_{(1-t)\gamma_{1}p_{1}+t\gamma_{2}p_{2}}\,, \quad t\in[0,1],$$ is a homotopy matching solutions $$\pi ^{-1}\left (F_{p_{1}}\right )$$ and $$\pi ^{-1}\left (F_{p_{2}}\right )$$, where the matching is potentially different from that given by H(t). Similarly, for an affine linear family, $$F_{p} = F^{\prime}_{p} + C$$, where $$F^{\prime}_{p}$$ is homogeneous in p and C is a constant system, we have $$H^{\prime}(t) = (1-t)\gamma_{1}F_{p_{1}} + t \gamma_{2}F_{p_{2}} = F^{\prime}_{(1-t)\gamma_{1}p_{1}+t\gamma_{2}p_{2}}+\left((1-t)\gamma_{1}+t\gamma_{2}\right)C.$$ We ignore the fact that H′(t) may go outside B for t ∈ (0, 1), since its rescaling, \begin{align*} H^{\prime\prime}(t) &= \frac{1}{(1-t)\gamma_{1}+t\gamma_{2}}H^{\prime}(t)\\ &=F^{\prime}_{\frac{(1-t)\gamma_{1}p_{1}+t\gamma_{2}p_{2}}{(1-t)\gamma_{1}+t\gamma_{2}}}+C = F_{\frac{(1-t)\gamma_{1}p_{1}+t\gamma_{2}p_{2}}{(1-t)\gamma_{1}+t\gamma_{2}}}\,, \quad t\in[0,1], \end{align*} does not leave B and clearly has the same homotopy paths. Note that H″(t) is well defined as $$(1-t)\gamma _{1}+t\gamma _{2}\neq 0$$ for all t ∈ [0, 1] for generic $$\gamma _{1}$$ and $$\gamma _{2}$$. One may use methods of numerical homotopy continuation, described, for instance in the book by Sommese & Wampler (2005, Section 2.3), to track the solutions as t changes from 0 to 1. In some situations the path in B may pass close to the branch locus D and numerical issues must be considered. Remark 2.4 If the family $$F_{p}$$ is nonlinear in the parameters p one has to take the parameter linear segment homotopy in the parameter space, i.e., $$H(t) = F_{(1-t)p_{1}+tp_{2}}$$, t ∈ [0, 1]. This does not change the overall construction; however, the freedom to replace the systems $$F_{p_{1}}$$ and $$F_{p_{2}}$$ at the ends of the homotopy with their scalar multiples as in Remark 2.3 is lost. 2.3 Graph of homotopies: main ideas Some readers may find it helpful to use the examples of Section 2.4 for graphical intuition as we introduce notation and definitions below. To organize the discovery of new solutions we represent the set of homotopies by a finite undirected graph G. Let E(G) and V(G) denote the edge and vertex set of G, respectively. Any vertex v in V(G) is associated to a point $$F_{p}$$ in the base space. An edge e in E(G) connecting $$v_{1}$$ and $$v_{2}$$ in V(G) is decorated with two complex numbers, $$\gamma _{1}$$ and $$\gamma _{2}$$, and represents the linear homotopy connecting $$\gamma _{1}F_{p_{1}}$$ and $$\gamma _{2} F_{p_{2}}$$ along a line segment (Remark 2.3). We assume that both $$p_{i}$$ and $$\gamma _{i}$$ are chosen so that the segments do not intersect the branch locus. Choosing these at random (see Section 5.1 for a possible choice of distribution) satisfies the assumption, since the exceptional set of choices where such intersections happen is contained in a real Zariski closed set, see Sommese & Wampler (2005, Lemma 7.1.3). We allow multiple edges between two distinct vertices but no loops, since the latter induce trivial homotopies. For a graph G to be potentially useful in a monodromy computation it must contain a cycle. Some of the general ideas behind the structure of a graph G are listed below. For each vertex $$v_{i}$$ we maintain a subset of known points $$Q_{i}\subset \pi ^{-1}\left (F_{p_{i}}\right )$$. For each edge e between $$v_{i}$$ and $$v_{j}$$ we record the two complex numbers $$\gamma _{1}$$ and $$\gamma _{2}$$, and we store the known partial correspondences $$C_{e}\!\subset\! \pi ^{-\!1}\!\left (F_{p_{i}}\right )\!\times \pi ^{-\!1}\!\left (F_{p_{j}}\right )$$ between known points $$Q_{i}$$ and $$Q_{j}$$. At each iteration we pick an edge and direction, track the corresponding homotopy starting with yet unmatched points, and update known points and correspondences between them. We may obtain the initial ‘knowledge’ as a seed pair$$(p_{0},x_{0})$$ by picking $$x_{0}\in {\mathbb{C}}^{n}$$ at random and choosing $$p_{0}$$ to be a generic solution of the linear system $$F_{p}(x_{0})=0$$. We list basic operations that result in transition between one state of our algorithm captured by G, $$Q_{i}$$ for $$v_{i}\in V(G)$$ and $$C_{e}$$ for e ∈ E(G) to another. For an edge $$e = v_{i}\overset{(\gamma _{1},\gamma _{2})}{\overleftrightarrow{\hphantom{(\gamma _{1},\gamma _{2})}}} v_{j}$$ consider the homotopy $$H^{(e)} = (1-t)\gamma_{1}F_{p_{i}} + t \gamma_{2}F_{p_{j}}$$ where $$\left (\gamma _{1},\gamma _{2}\right )\in {\mathbb{C}}^{2}$$ is the label of e. Take start points $$S_{i}$$ to be a subset of the set of known points $$Q_{i}$$ that do not have an established correspondence with points in $$Q_{j}$$. Track $$S_{i}$$ along $$H^{(e)}$$ for t ∈ [0, 1] to get $$S_{j} \subset \pi ^{-1}\left (F_{p_{j}}\right )$$. Extend the known points for $$v_{j}$$, that is, $$Q_{j} := Q_{j} \cup S_{j}$$ and record the newly established correspondences. Add a new vertex corresponding to $$F_{p}$$ for a generic p ∈ B ∖ D. Add a new edge $$e = v_{i}\overset{(\gamma _{1},\gamma _{2})}{\overleftrightarrow{\hphantom{(\gamma _{1},\gamma _{2})}}}v_{j}$$ between two existing vertices decorated with generic $$\gamma _{1},\gamma _{2}\!\in\! {\mathbb{C}}$$. At this point a reader who is ready to see an algorithm based on these ideas may skip to Algorithm 3.1. 2.4 Graph of homotopies: examples We demonstrate the idea of graphs of homotopies, the core idea of the MS framework, by giving two examples. Example 2.5 Figure 1 shows a graph G with two vertices and three edges embedded in the base space B with paths partially lifted to the solution variety, which is a covering space with three sheets. The two fibres $$\{x_{1},x_{2},x_{3}\}$$ and $$\{y_{1},y_{2},y_{3}\}$$ are connected by three partial correspondences induced by the liftings of three egde paths. Note that several aspects in this illustration are fictional. There is only one branch point in the actual complex base space B that we would like the reader to imagine. The visible self-intersections of the solution variety V are an artifact of drawing the picture in the real space. Also, in practice we use homotopy paths as simple as possible, however, here the paths are more involved for the purpose of distinguishing them in print. An algorithm that we envision may hypothetically take the following steps: (1) seed the first fibre with $$x_{1}$$; (2) use a lifting of edge $$e_{a}$$ to get $$y_{1}$$ from $$x_{1}$$; (3) use a lifting of edge $$e_{b}$$ to get $$x_{2}$$ from $$y_{1}$$; (4) use a lifting of edge $$e_{c}$$ to get $$y_{2}$$ from $$x_{1}$$; (5) use a lifting of edge $$e_{a}$$ to get $$x_{3}$$ from $$y_{2}$$. Note that it is not necessary to complete the correspondences (a), (b) and (c). Doing so would require tracking nine continuation paths, while the hypothetical run above uses only four paths to find a fibre. Example 2.6 Figure 2 illustrates two partial correspondences associated to two edges $$e_{a}$$ and $$e_{b}$$, both connecting two vertices $$v_{1}$$ and $$v_{2}$$ in V(G). Each vertex $$v_{i}$$ stores the array of known points $$Q_{i}$$, which are depicted in solid. Both correspondences in the picture are subsets of a perfect matching, a one-to-one correspondence established by a homotopy associated to the edge. Note that taking the set of start points $$S_{1}=\{x_{3}\}$$ and following the homotopy $$H^{(e_{a})}$$ from left to right is guaranteed to discover a new point in the second fibre. On the other hand, it is impossible to obtain new knowledge by tracking $$H^{(e_{a})}$$ from right to left. Homotopy $$H^{(e_{b})}$$ has a potential to discover new points if tracked in either direction. We can choose $$S_{1}=\{x_{1},x_{3}\}$$ as the start points for one direction and $$S_{2}=\{y_{3}\}$$ for the other. In this scenario, following the homotopy from left to right is guaranteed to produce at least one new point, while going the other way may either deliver a new point or just augment the correspondences between the already known points. If the correspondences in (a) and (b) are completed to one-to-one correspondences of the fibres taking the homotopy induced by the edge $$e_{a}$$ from left to right followed by the homotopy induced by edge $$e_{b}$$ from right to left would produce a permutation. However, the group generated by this permutation has to stabilize $$\{x_{2}\}$$, therefore, it would not act transitively on the fibre of $$v_{1}$$. One could also imagine a completion such that the given edges would not be sufficient to discover $$x_{5}$$ and $$y_{4}$$. Fig. 1. View largeDownload slide Selected liftings of three edges connecting the fibres of two vertices and induced correspondences. Fig. 1. View largeDownload slide Selected liftings of three edges connecting the fibres of two vertices and induced correspondences. Fig. 2. View largeDownload slide Two partial correspondences induced by edges $$e_{a}$$ and $$e_{b}$$ for the fibres of the covering map of degree d = 5 in Example 2.6. Fig. 2. View largeDownload slide Two partial correspondences induced by edges $$e_{a}$$ and $$e_{b}$$ for the fibres of the covering map of degree d = 5 in Example 2.6. Fig. 3. View largeDownload slide Graphs for the flower(4,2) strategy and completeGraph(5,1). Fig. 3. View largeDownload slide Graphs for the flower(4,2) strategy and completeGraph(5,1). In our algorithm we record and use correspondences; however, they are viewed as a secondary kind of knowledge. In particular, in Section 3.2.4 we develop heuristics driven by edge potential functions which look to maximize the number of newly discovered solutions, in other words, to extend the primary knowledge in some greedy way. 3. Algorithms and strategies The operations listed in Section 2.3 give a great deal of freedom in the discovery of solutions. However, not all strategies for applying these operations are equally efficient. We distinguish between static strategies, where the graph is fixed throughout the discovery process (only basic operation 1 of Section 2.3 is used) and dynamic strategies, where vertices and edges may be added (operations 2 and 3). 3.1 A naive dynamic strategy To visualize this strategy in our framework jump ahead and to the flower graph in Fig. 3. Start with the seed solution at the vertex $$v_{0}$$ and proceed creating loops as petals in this graph: e.g., use basic operations 2 and 3 to create $$v_{1}$$ and two edges between $$v_{0}$$ and $$v_{1}$$, track the known solutions at $$v_{0}$$ along the new petal to potentially find new solutions at $$v_{0}$$, then ‘forget’ the petal and create an entirely new one in the next iteration. This strategy populates the fibre $$\pi ^{-1}\left (F_{p_{0}}\right )$$, but how fast? Assume the permutation induced by a petal permutation on $$\pi ^{-1}\left (F_{p_{1}}\right )$$ is uniformly distributed. Then for the first petal the probability of finding a new solution equals (d − 1)/d where $$d=\left |\pi ^{-1}\left (F_{p_{1}}\right )\right |$$, which is large. However, for the other petals the probability of arriving at anything new at the end of one tracked path decreases as the known solution set grows. Finding the expected number of iterations (petals) to discover the entire fibre is equivalent to solving the coupon collector’s problem. The number of iterations is d ℓ(d) where $$\ell (d):={1\over 1}+{1\over 2}+\cdots +{1\over d}$$. The values of ℓ(d) can be regarded as lower and upper sums for two integrals of the function $$x\mapsto x^{-1}$$, leading to the bounds $$\ln (d+1)\leq \ell (d)\leq \ln (d)+1$$. Simultaneously tracking all known points along a petal gives a better complexity, since different paths cannot lead to the same solution. We remark that the existing implementations of numerical irreducible decomposition in Bertini (Bates et al., 2013), PHCpack (Verschelde, 1999) and NumericalAlgebraicGeometry for Macaulay2 (Leykin, 2011) that use monodromy are driven by a version of the naive dynamic strategy. 3.2 Static graph strategies It turns out to be an advantage to reuse the edges of the graph. In a static strategy the graph is fixed and we discover solutions according to the following algorithm. The algorithm can be specialized in several ways. We may choose the graph G, specify a stopping criterion stop, choose a strategy for picking the edge e = (j, k). We address the first choice in Section 3.2.1 by listing several graph layouts that can be used. Stopping criteria are discussed in Section 3.2.2 and Section 3.2.3, while strategies for selecting an edge are discussed in Section 3.2.4. Remark 3.2 We notice that if the stopping criterion is never satisfied the number of paths being tracked by Algorithm 3.1 is at most d|E(G)|, where d is the number of solutions of a generic system. 3.2.1 Two static graph layouts We present two graph layouts to be used for the static strategy (Fig. 3). flower(s,t) The graph consists of a central node$$v_{0}$$ and s additional vertices (number of petals), each connected to $$v_{0}$$ by t edges. completeGraph(s,t) The graph has s vertices. Every pair of vertices is connected by t edges. 3.2.2 Stopping criterion if a solution count is known Suppose the cardinality of the fibre $$\pi ^{-1}\left (F_{p}\right )$$ for a generic value of p is known. Then a natural stopping criterion for our algorithm is to terminate when the set of known solutions $$Q_{i}$$ at any node i reaches that cardinality. In particular, for a generic sparse system with fixed monomial support, we can rely on this stopping criterion due to the BKK bound (Bernstein, 1975) that can be obtained by a mixed volume computation. 3.2.3 Stopping criterion if no solution count is known For a static strategy one natural stopping criterion is saturation of the known solution correspondences along all edges. In this case the algorithm simply can’t derive any additional information. It also makes sense to consider a heuristic stopping criterion based on stabilization. The algorithm terminates when no new points are discovered in a fixed number of iterations. This avoids saturating correspondences unnecessarily. In particular, this could be useful if a static strategy algorithm is a part of the dynamic strategy of Section 3.3. Remark 3.3 In certain cases it is possible to provide a stopping criterion using the trace test (Sommese et al., 2002; Leykin et al., 2016). This is particularly useful when there is an equation in the family $$F_{p}(x)=0$$ that describes a generic hypersurface in the parameter space, e.g., an affine linear equation with indeterminate coefficients. In full generality, one could restrict the parameter space to a generic line and, hence, restrict the solution variety to a curve. Now, thinking of $$F_{p}(x)=0$$ as a system of equations bihomogeneous in p and x, one can use the multihomogeneous trace test (Hauenstein & Rodriguez, 2015; Leykin et al., 2016). We note that the multihomogeneous trace test complexity depends on the degree of the solution variety, which may be significantly higher than the degree d of the covering map, where the latter is the measure of complexity for the main problem in our paper. For instance, the system (7) corresponding to the reaction network in Fig. 4 has four solutions, but an additional set of 11 points is necessary to execute the trace test. See example-traceCRN.m2 at (Duff et al.). 3.2.4 Edge-selection strategy We propose two methods for selecting the edge e in Algorithm 3.1. The default is to select an edge and direction at random. A more sophisticated method is to select an edge and a direction based on the potential of that selection to deliver new information: see the discussion in Example 2.6. Let $$e = v_{i} \overset{(\gamma _{1},\gamma _{2})}{\overleftrightarrow{\hphantom{(\gamma _{1},\gamma _{2})}}}v_{j}$$ be an edge considered in the direction from $$v_{i}$$ to $$v_{j}$$. potentialLowerBound equals the minimal number of new points guaranteed to be discovered by following a chosen homotopy using the maximal batch of starting points $$S_{i}$$. That is, it equals the difference between the numbers of known unmatched points $$\left (|Q_{i}|-|C_{e}|\right )-\left (\left |Q_{j}\right |-|C_{e}|\right ) = \left |Q_{i}\right |-\left |Q_{j}\right |$$ if this difference is positive and 0 otherwise. potentialE equals the expected number of new points obtained by tracking one unmatched point along e. This is the ratio $$d-|Q_{j}|\over d-|C_{e}|$$ of undiscovered points among all unmatched points if $$\left |Q_{i}\right |-|C_{e}|>0$$ and 0 otherwise. Note that potentialE assumes we know the cardinality of the fibre, while potentialLowerBound does not depend on that piece of information. There is a lot of freedom in choosing potentials in our algorithmic framework. The two above potentials are natural ‘greedy’ choices that are easy to describe and implement. It is evident from our experiments (Section 6.1.1) that they may order edges differently resulting in varying performance. 3.3 An incremental dynamic graph strategy Consider a dynamic strategy that amounts to augmenting the graph once one of the above ‘static’ criteria terminates Algorithm 3.1 for the current graph. One simple way to design a dynamic stopping criterion, we call it dynamic stabilization, is to decide how augmentation is done and fix the number of augmentation steps that the algorithm is allowed to make without increasing the solution count. A dynamic strategy, which is simple to implement, is one that starts with a small graph G and augments it if necessary. We emphasize that the criteria described in this subsection and parts of Section 3.2.3 are heuristic and there is a lot of freedom in designing such. In Section 6.2 we successfully experiment using a static stabilization criterion with some examples, for which the solution count is generally not known. 4. Statistical analysis The directed cycles in the graph G starting and ending at a vertex $$v_{1}$$ give elements of the fundamental group $$\pi _{1}(B\setminus D)$$, which correspond to the elements of the monodromy subgroup M(G) of the monodromy group $$M(\pi _{1}(B\setminus D))$$. The latter is a subgroup of $${\mathrm{S}}_{d}$$, where $$d = \left |\pi ^{-1}\left (F_{p_{1}}\right )\right |$$. For example, if $$G={\tt completeGraph}(2,j+1)$$ then the j cycles produced by edges $$e_{1}$$ and $$e_{2}$$, $$e_{2}$$ and $$e_{3},\dots,e_{j}$$ and $$e_{j+1}$$ suffice to generate M(G). The minimal number j of cycles necessary to generate M(G) in the general case is $$\beta _{1}(G)$$, the first Betti number of G as a topological space. For the purpose of simplifying statistical analysis, we assume that picking a random decorated graph G with $$j=\beta _{1}(G)$$ induces uniformly and independently distributed permutations $$\sigma _{1}, \ldots , \sigma _{j} \in {\mathrm{S}}_{d}$$, where $${\mathrm{S}}_{d}$$ is the symmetric group acting on the fibre $$\pi ^{-1}\left (F_{p_{1}}\right )$$. It would be hard in practice to achieve uniformity even when the monodromy group is a full symmetric group: see Section 5.1. 4.1 The probability of a transitive action Suppose the number of solutions d is known and stop(d) denotes the corresponding stopping criterion. Our aim is to analyze the probability of producing the full solution set via Algorithm 3.1 or, equivalently, the probability of $${\tt dynamicMonodromySolve}(G,x_{1},{\tt stop(d)},{\tt augment})$$ terminating after at most j iterations, assuming that $$\beta _{1}(G)=j$$ at the jth iteration. This equals the probability of $$\langle \sigma _{1}, \ldots , \sigma _{j}\rangle$$ acting transitively, i.e., $$\Pr \left [X_{d}\leq j\right ]$$, where $$X_{d}$$ is the random variable $$X_{d} = \displaystyle\inf \left\{ i \in \mathbb{N} \mid \langle \sigma_{1}, \ldots, \sigma_{i} \rangle \textrm{ is transitive} \right\}\!.$$ When d > 1 we have $$\Pr [X_{d} =0]=0$$, while $$\Pr [X_{d} =1 ]$$ is proportional to the number of d-cycles in the monodromy group. When the monodromy group is full symmetric, we can compute and give asymptotic estimates for the distribution of $$X_{d}.$$ The following theorem is a generalization of a result by Dixon, regarding the case j = 2. The proof we give in Section 4.2 follows the strategy of the study by Dixon (1969). Theorem 4.1 For j ≥ 2, $$\Pr \left [ X_{d} \le j \right ] = 1 - d^{1-j} + R_{j} (d),$$ where the error term $$R_{j}$$ satisfies $$\left |R_{j}(d)\right |=O\left (d^{-j}\right )$$. Remark 4.2 As a corollary one can deduce that the expected value of $$X_{d}$$ is asymptotically finite and $$\mathrm{E}\left [X_{d}\right ]\to 2$$ as $$d\to \infty$$. The numerical approximations in Table 1 show that $$\mathrm{E}\left [X_{d}\right ]\leq 2.1033$$ for all d. Moreover, the proof in Section 4.2 implies that $$\left |R_{j}(d)\right |<C\left ({d\over 2}\right )^{-j}$$ with the constant C not depending on j. Therefore, $$\Pr \left [ X_{d}> j \right ]$$ decays exponentially with j. Under the idealistic assumption that new cycles in the graph lead to independently and uniformly distributed permutations of the fibre $$Q_{1}$$ the expected Betti number needed for completion in Algorithm 3.4 is at most 2.1033. If we assume that augment increases the Betti number by one by adding at most a fixed number of edges then the expected number of tracked paths is linear in d. Remark 4.3 We point out that Babai (1989) proved Dixon’s conjecture stating that the subgroup of $${\mathrm{S}}_{d}$$ generated by two random permutations is $${\mathrm{S}}_{d}$$ or $${\mathrm{A}}_{d}$$ with probability $$1 - d^{-1} + O\left (d^{-2}\right )$$. This shows that other subgroups are rare. However, it is easy to construct families with a transitive monodromy group that is neither full symmetric nor alternating. For example, take $${x_{1}^{2}}-c_{1}={x_{2}^{2}}-c_{2}=0$$ with irreducible solution variety, and four solutions for generic choices of $$c_{1}$$ and $$c_{2}$$. Tracking two solutions with the same $$x_{1}$$ coordinate, as $$c_{1}$$ and $$c_{2}$$ vary, the moving points on the tracked paths will continue to have equal projections to $$x_{1}$$. The monodromy group is $${\mathbb{Z}}_{2}\times {\mathbb{Z}}_{2}$$. We reiterate that generators for the monodromy group are seldomly known a priori. Computing them is likely to be prohibitively expensive, and the probability distribution with which our algorithm picks elements of the monodromy group is unknown, as it is prohibitively hard to analyze. 4.2 Proof of a generalization of Dixon’s theorem Fix any integer j ≥ 2. We wish to prove Theorem 4.1 by estimating the quantity $$t_{d} = \Pr \left[ \langle \sigma_{1}, \sigma_{2}, \ldots, \sigma_{j} \rangle \textrm{ is transitive}\right]\!,$$ where $$\sigma _{1}, \ldots , \sigma _{j}$$ are independent and uniformly distributed on $${\mathrm{S}}_{d}$$. Suppose we partition the set {1, 2, … , d} in such a way that there are $$k_{i}$$ classes of size i for each 1 ≤ i ≤ d. All such partitioning schemes are indexed by the set $$K_{d} = \left\{ \vec{k} \in{\mathbb{N}}^{d} \left|\ \sum \right. i \, k_{i} = d \right\}\!.$$ The number of partitions corresponding to each $$\vec{k} \in K_{d}$$ is $${d!}/\big (\!\prod _{i=1}^{d} (i!\!)^{k_{i}} \cdot k_{i}!\!\big )$$. For each $$\vec{k} \in K_{d},$$ the partition given by the orbits of $$\langle \sigma _{1}, \ldots , \sigma _{j} \rangle$$ is $$\vec{k}$$-indexed precisely when this group acts transitively on all classes of some partition associated to $$\vec{k}.$$ The number of tuples in $${S_{i}^{j}}$$ with coordinates generating a group acting transitively on $$\{1,\dots ,i\}$$ is $$t_{i} \,\left (i!\!\right )^{j}$$. Thus, we may count the set $$\underbrace{ {\mathrm{S}}_{d} \times \cdots \times {\mathrm{S}}_{d}}_{j}$$ as \begin{align*} \left(d!\!\right)^{j} &= \displaystyle\sum_{\vec{k} \in K_{d}} \displaystyle\frac{d!}{\prod_{i=1}^{d} (i!\!)^{k_{i}} \cdot k_{i}!} \cdot \displaystyle\prod_{i=1}^{d} \big( t_{i} \left(i!\!\right)^{j} \big)^{k_{i}} \\[-2pt] &= d! \cdot \displaystyle\sum_{\vec{k} \in K_{d}}\displaystyle\prod_{i=1}^{d} \displaystyle\frac{\left( t_{i} \, (i!\!)^{j-1} \right)^{k_{i}}}{k_{i}!}. \end{align*} Let $$\widehat{F}$$ denote the generating function of the sequence $$F(d) = (d!\!)^{j-1}.$$ Note the formal identity $$\exp \left( \displaystyle\sum_{i=0}^{\infty } y_{i} x^{i} \right) = \displaystyle\sum_{d=0}^{\infty} x^{d} \, \sum_{\vec{k} \in K_{d} } \prod_{i=1}^{d} \frac{y_{i}^{k_{i}}}{k_{i}!},$$ which follows by letting g(x) denote the right-hand side as a formal power series in x, $$f(x) = \sum y_{i} x^{i},$$ and noting the equivalent form f′g = g′ with $$f(0)=y_{0}.$$ We have \begin{align*} \displaystyle\sum_{d=1}^{\infty} d\cdot (d!\!)^{j-1} \, x^{d-1} &= \frac{d}{d x} \, \widehat{F} (x) \\ &= \frac{d}{dx} \, \exp \left( \displaystyle\sum_{i=1}^{\infty } t_{i} \, (i!\!)^{j-1} x^{i} \right)\\ &= \left( \displaystyle\sum_{d=0}^{\infty} (d!\!)^{j-1} \, x^{d} \right) \cdot \displaystyle\sum_{i=1}^{\infty } i \cdot t_{i} \, (i!\!)^{j-1} x^{i-1} \\ &= \displaystyle\sum_{d^{\prime}=1}^{\infty } x^{d^{\prime}-1} \, \left( \displaystyle\sum_{i=1}^{d^{\prime}} i \cdot t_{i} \, \big( i! \cdot \left(d^{\prime}-i\right)! \big)^{j-1} \right)\!, \end{align*} where the first equation follows by formal differentiation of the power series $$\widehat{F}$$, the second from the two identities above with properly substituted values for $$y_{i}$$, the third by applying the chain rule and the definition of $$\widehat{F}$$ and the fourth by rearranging terms by index substitution $$d^{\prime}=i+d$$. Upon equating coefficients of $$x^{d-1}$$ for $$d=1,\dots$$ we obtain $$d = \displaystyle\sum_{i=1}^{d} \binom{d}{i}^{1-j} \, i\, t_{i}.$$ (4.1) Table 1 Numerical approximations of $$t_{d}$$—the probability of the j random permutations acting transitively on a fibre of size d for j = 2, 3, 4. After computing these values for larger j, a numerical approximation of $$\mathrm{E}\left [X_{d}\right ]$$ is extracted d j = 2 j = 3 j = 4 $$\mathrm{E}\left [X_{d}\right ]$$ 1 1 1 1 0 2 0.75 0.875 0.9375 2 3 0.72222222 0.89814815 0.96450617 2.10000000 4 0.73958333 0.93012153 0.98262080 2.10329381 5 0.76833333 0.95334722 0.99115752 2.08926525 10 0.88180398 0.98954768 0.99898972 2.02976996 20 0.94674288 0.99747856 0.99987487 2.00591026 30 0.96536852 0.99888488 0.99996295 2.00245160 d j = 2 j = 3 j = 4 $$\mathrm{E}\left [X_{d}\right ]$$ 1 1 1 1 0 2 0.75 0.875 0.9375 2 3 0.72222222 0.89814815 0.96450617 2.10000000 4 0.73958333 0.93012153 0.98262080 2.10329381 5 0.76833333 0.95334722 0.99115752 2.08926525 10 0.88180398 0.98954768 0.99898972 2.02976996 20 0.94674288 0.99747856 0.99987487 2.00591026 30 0.96536852 0.99888488 0.99996295 2.00245160 Table 1 Numerical approximations of $$t_{d}$$—the probability of the j random permutations acting transitively on a fibre of size d for j = 2, 3, 4. After computing these values for larger j, a numerical approximation of $$\mathrm{E}\left [X_{d}\right ]$$ is extracted d j = 2 j = 3 j = 4 $$\mathrm{E}\left [X_{d}\right ]$$ 1 1 1 1 0 2 0.75 0.875 0.9375 2 3 0.72222222 0.89814815 0.96450617 2.10000000 4 0.73958333 0.93012153 0.98262080 2.10329381 5 0.76833333 0.95334722 0.99115752 2.08926525 10 0.88180398 0.98954768 0.99898972 2.02976996 20 0.94674288 0.99747856 0.99987487 2.00591026 30 0.96536852 0.99888488 0.99996295 2.00245160 d j = 2 j = 3 j = 4 $$\mathrm{E}\left [X_{d}\right ]$$ 1 1 1 1 0 2 0.75 0.875 0.9375 2 3 0.72222222 0.89814815 0.96450617 2.10000000 4 0.73958333 0.93012153 0.98262080 2.10329381 5 0.76833333 0.95334722 0.99115752 2.08926525 10 0.88180398 0.98954768 0.99898972 2.02976996 20 0.94674288 0.99747856 0.99987487 2.00591026 30 0.96536852 0.99888488 0.99996295 2.00245160 Remark 4.4 Equation (1) gives a list of linear equations in the probabilities $$t_{1},t_{2},\dots$$ allowing us to successively determine these values by backward substitution. In Table 1 we list some solutions for j = 2, 3, 4. To complete the proof of Theorem 4.1, we introduce, as in Dixon’s proof, the auxiliary quantities \begin{equation*} r_{d} = d\, \left(1 - t_{d}\right)\ \ \textrm{and}\ \ \ c_{d} = \displaystyle\sum_{i=1}^{d-1} \binom{d}{i}^{1-j} \, i . \end{equation*} Noting that $$\binom{d}{i}^{1-j}i+\binom{d}{d-i}^{1-j}(d-i)={d\over 2}\left (\binom{d}{i}^{1-j}+\binom{d}{d-i}^{1-j}\right )$$, we have \begin{align} {c_{d}\over d} =\frac{1}{2} \cdot \displaystyle\sum_{i=1}^{d-1} \binom{d}{i}^{1-j} &= d^{1-j} + \left[ \binom{d}{2}^{1-j} + \frac{1}{2} \displaystyle\sum_{i=3}^{d-3} \binom{d}{i}^{1-j} \right] \nonumber\\ &\leq d^{1-j} + \left[ \binom{d}{2}^{1-j} + \frac{1}{2} (d-5) \binom{d}{3}^{1-j} \right]\!. \end{align} (4.2) From j ≥ 2 it follows that the bracketed expression in (2) is $$O\left (d^{-j}\right ).$$ Using (1), $$t_{i}=1-{r_{i}\over i}\leq 1-{r_{i}\over d}$$ and the definition of $$c_{d}$$ we may bound $$r_{d}$$: \begin{align*} r_{d} &= d(1-t_{d})=-dt_{d}+d=-dt_{d}+\displaystyle\sum_{i=1}^{d} \binom{d}{i}^{1-j} \, i\cdot t_{i} \\ &=\displaystyle\sum_{i=1}^{d-1} \binom{d}{i}^{1-j} \, i \cdot t_{i} \leq\displaystyle\sum_{i=1}^{d-1} \binom{d}{i}^{1-j} \, i \cdot \left(1-{r_{i}\over d}\right)\nonumber \\ &= \sum_{i=1}^{d-1} \binom{d}{i}^{1-j}i -{1\over d}\displaystyle\sum_{i=1}^{d-1} r_{i} \, \binom{d}{i}^{1-j}= c_{d} -{1\over d}\displaystyle\sum_{i=1}^{d-1} r_{i} \, \binom{d}{i}^{1-j}\leq c_{d}. \end{align*} (3) (4) To bound the error term $$R_{j}(d):=t_{d} - \left (1-d^{1-j}\right )$$, we consider first the case where its sign is positive. Expanding $$t_{d}=1-{r_{d}\over d}$$ using (3) above and $$i\, t_{i}=i-{r_{i}}$$, \begin{align*} t_{d} - \big(1-d^{1-j}\big) &= 1-\sum_{i=1}^{d-1} \binom{d}{i}^{1-j}{it_{i}\over d}-1+d^{1-j}=d^{1-j}-\sum_{i=1}^{d-1} \binom{d}{i}^{1-j}{it_{i}\over d}\\ &=d^{1-j}-\sum_{i=1}^{d-1} \binom{d}{i}^{1-j}{i\over d}+\sum_{i=1}^{d-1} \binom{d}{i}^{1-j}{r_{i}\over d}\\ &=\underbrace{d^{1-j}-d^{1-j}\left({1\over d}+{d-1\over d}\right)}_{=0} \, - \, \sum_{i=2}^{d-2} \binom{d}{i}^{1-j}{i\over d} \, + \, \sum_{i=1}^{d-1} \binom{d}{i}^{1-j}{r_{i}\over d}, \end{align*} and we may focus on the last summation to get, for d ≥ 2, \begin{align*} t_{d} - \big(1-d^{1-j}\big) &\leq d^{1-j}{r_{d-1}\over d-1}+\sum_{i=2}^{d-2} \binom{d}{i}^{1-j}{r_{i}\over d}\qquad\qquad(\text{note }r_{1} = 0)\\ &\leq d^{1-j}{c_{d-1}\over d-1}+\sum_{i=2}^{d-2} \binom{d}{2}^{1-j}{c_{i}\over d}\qquad\qquad\qquad(\text{by }(4))\\ &\leq d^{1-j}{c_{d-1}\over d-1}+\left({d^{2}\over 4}\right)^{1-j}\, \sum_{i=2}^{d-2}{c_{i}\over i}\\ &=O\big(d^{2-2j}\big)+O\big(d^{(2-2j)+(2-j)}\big) \\ &= O\big(d^{-j}\big).\qquad\qquad\qquad\qquad\qquad(\text{since }j\geq 2) \end{align*} The case where $$R_{j}(d)\le 0$$ may be handled similarly, using $$t_{d}=1-{r_{d}\over d}$$, $$r_{d}\leq c_{d}$$ and what we know about the content of the bracket in (2). \begin{align*} -t_{d} + \big(1-d^{1-j}\big)&= -\left(1-{r_{d} \over d}\right)+\big(1-d^{1-j}\big)={r_{d} \over d}-d^{1-j}\leq{c_{d} \over d}-d^{1-j}\\ &\leq \left[ \binom{d}{2}^{1-j} + \frac{1}{2} (d-5) \binom{d}{3}^{1-j} \right] = O\big(d^{-j}\big). \end{align*} 5. Implementation We implement the package MonodromySolver in Macaulay2 (Grayson & Stillman) using the functionality of the package NumericalAlgebraicGeometry (Leykin, 2011). The source code and examples used in the experiments in the next section are available at Duff et al. The main function monodromySolve realizes Algorithms 3.1 and 3.4, see the documentation for details and many options. The tracking of homotopy paths in our experiments is performed with the native routines implemented in the kernel of Macaulay2, however, NumericalAlgebraicGeometry provides an ability to outsource this core task to an alternative tracker (PHCpack or Bertini). Main auxiliary functions—createSeedPair, sparseSystemFamily, sparseMonodromySolve and solveSystemFamily—are there to streamline the user’s experience. The last two are blackbox routines that don’t assume any knowledge of the framework described in this paper. The overhead of managing the data structures is supposed to be negligible compared to the cost of tracking paths. However, since our implementation uses the interpreted language of Macaulay2 for other tasks, this overhead could be sizable (up to 10% for large examples in Section 6). Nevertheless, most of our experiments are focused on measuring the number of tracked paths as a proxy for computational complexity. Remark 5.1 This paper’s discussion focuses on linear parametric systems with a nonempty dominant component. However, the implementation works for other cases where our framework can be applied. For instance, if the system is linear in parameters, but has no dominant component, there may still be a unique ‘component of interest’ with a straightforward way to produce a seed pair. This is so, for instance, in the problem of finding the degree of the variety $${\mathop{SO}}(n)$$, which we use in Table 7. The point x is restricted to $${\mathop{SO}}(n)$$, the special orthogonal group, which is irreducible as a variety. This results in a unique ‘component of interest’ in the solution variety, the one that projects onto $${\mathop{SO}}(n)$$, see the study by Brandt et al. (2017) for details. In a yet more general case of a system that is nonlinear in parameters, it is still possible to use our software. We outline the theoretical issues one would need to consider in 4 of Section 7. 5.1 Randomization Throughout the paper we refer to random choices we make that we assume avoid various nongeneric loci. For implementation purposes we make simple choices. For instance, the vertices of the graph get distributed uniformly in a cube in the base space with the exception of the seeded vertex: createSeedPair picks $$(p_{0},x_{0})\in B\times {\mathbb{C}}^{n}$$ by choosing x uniformly in a cube, then choosing $$p_{0}$$ uniformly in a box in the subspace $$\left \{ p \mid F_{p}(x)=0 \right \}$$. A choice of probability distribution on B translates to some (discrete) distribution on the symmetric group $${\mathrm{S}}_{d}$$. However, it is simply too hard to analyze—there are virtually no studies in this direction. We make the simplest possible assumption of uniform distribution on $${\mathrm{S}}_{d}$$ in order to perform the theoretical analysis in Section 4, and shed some light on why our framework works well. There is an interesting, more involved, alternative to this assumption in the studies by Galligo & Poteaux (2011) and Galligo & Miclo (2012), which relies on the intuition in the case n = 1. 5.2 Solution count The BKK bound, computed via mixed volume, is used as a solution count in the examples of sparse systems in Sections 6.1.1 and 6.1.2. In the latter we compute mixed volume via a closed formula that involves permanents, while the former relies on general algorithms implemented in several software packages. Our current implementation uses PHCpack (Verschelde, 1999), which incorporates the routines of MixedVol further developed in Hom4PS-2 (Lee et al., 2008). Other alternatives are pss5 (Malajovich, 2017) and Gfanlib (Jensen, 2016). While some implementations are randomized the latter uses symbolic perturbations to achieve exactness. The computation of the mixed volume is not a bottleneck in our algorithm. The time spent in that preprocessing stage is negligible compared to the rest of the computation. 5.3 Certification The reader should realize that the numerical homotopy continuation we use is driven partly by heuristics. As a post-processing step we can certify (i.e., formally prove) the completeness and correctness of the solution set to a polynomial system computed with our main method. This is possible in the scenario when the parametric system is square, all solutions are regular (the Jacobian of the system is invertible), and the solution count is known. We can use Smale’s $$\alpha$$-theory (Blum et al., 1998, Section 8) to certify an approximation to a regular solution of a square system. In a Macaulay2 package NumericalCertification, we implement a numerical version of an $$\alpha$$-test after finding an approximate solution to certify that our solution is an approximate zero in a rigorous sense. One of the main functions of NumericalCertification is certifySolutions, which determines whether the given solution is an approximate zero of the given polynomial system. It also produces an upper bound on the distance from the approximation to the exact solution to which it is associated. See paper-examples/example-NashCertify.m2 at (Duff et al.), which is an example of an $$\alpha$$-test application to the solutions of a problem described in Section 6.1.2. In the implementation of certification all arithmetic and linear algebra operations are done over the field of Gaussian rationals, $${\mathbb{Q}}[ {\mathbf{i}}]/\big ( {\mathbf{i}}^{2}+1\big )$$. To use this certification method we first convert the coefficients of the system to Gaussian rationals, then perform certification numerically. See the study by Hauenstein & Sottile (2012) for a stand-alone software package alphaCertified and detailed implementation notes. 6. Experiments In this section we first report on experiments with our implementation and various examples in Sections 6.1 and 6.2. We then investigate the completion rate of Algorithm 3.1 in Section 6.3. Finally, we compare against other software in Section 6.4. Table 2 Cyclic 7 experimental results for the flower strategy (#vertices-1, edge multiplicity) (3,2) (4,2) (5,2) (3,3) (4,3) |E(G)| 6 8 10 9 12 $$\beta _{1}(G)$$ 3 4 5 6 8 |E(G)| ⋅ 924 5544 7392 9240 8316 11088 completion rate 100% 100% 100% 100% 100% Random Edge 5119 6341 7544 6100 7067 potentialLowerBound 5252 6738 8086 6242 7886 potentialE 4551 5626 6355 4698 5674 (#vertices-1, edge multiplicity) (3,2) (4,2) (5,2) (3,3) (4,3) |E(G)| 6 8 10 9 12 $$\beta _{1}(G)$$ 3 4 5 6 8 |E(G)| ⋅ 924 5544 7392 9240 8316 11088 completion rate 100% 100% 100% 100% 100% Random Edge 5119 6341 7544 6100 7067 potentialLowerBound 5252 6738 8086 6242 7886 potentialE 4551 5626 6355 4698 5674 View Large Table 2 Cyclic 7 experimental results for the flower strategy (#vertices-1, edge multiplicity) (3,2) (4,2) (5,2) (3,3) (4,3) |E(G)| 6 8 10 9 12 $$\beta _{1}(G)$$ 3 4 5 6 8 |E(G)| ⋅ 924 5544 7392 9240 8316 11088 completion rate 100% 100% 100% 100% 100% Random Edge 5119 6341 7544 6100 7067 potentialLowerBound 5252 6738 8086 6242 7886 potentialE 4551 5626 6355 4698 5674 (#vertices-1, edge multiplicity) (3,2) (4,2) (5,2) (3,3) (4,3) |E(G)| 6 8 10 9 12 $$\beta _{1}(G)$$ 3 4 5 6 8 |E(G)| ⋅ 924 5544 7392 9240 8316 11088 completion rate 100% 100% 100% 100% 100% Random Edge 5119 6341 7544 6100 7067 potentialLowerBound 5252 6738 8086 6242 7886 potentialE 4551 5626 6355 4698 5674 View Large Table 3 Cyclic 7 experimental results for the completeGraph strategy (#vertices, edge multiplicity) (2,3) (2,4) (2,5) (3,2) (4,1) |E(G)| 3 4 5 6 6 $$\beta _{1}(G)$$ 2 3 4 4 3 |E(G)| ⋅ 924 2772 3698 4620 5544 5544 completion rate 65% 80% 90% 100% 100% Random Edge 2728 3296 3947 4805 5165 potentialLowerBound 2727 3394 3821 4688 5140 potentialE 2692 2964 2957 3886 4380 (#vertices, edge multiplicity) (2,3) (2,4) (2,5) (3,2) (4,1) |E(G)| 3 4 5 6 6 $$\beta _{1}(G)$$ 2 3 4 4 3 |E(G)| ⋅ 924 2772 3698 4620 5544 5544 completion rate 65% 80% 90% 100% 100% Random Edge 2728 3296 3947 4805 5165 potentialLowerBound 2727 3394 3821 4688 5140 potentialE 2692 2964 2957 3886 4380 Table 3 Cyclic 7 experimental results for the completeGraph strategy (#vertices, edge multiplicity) (2,3) (2,4) (2,5) (3,2) (4,1) |E(G)| 3 4 5 6 6 $$\beta _{1}(G)$$ 2 3 4 4 3 |E(G)| ⋅ 924 2772 3698 4620 5544 5544 completion rate 65% 80% 90% 100% 100% Random Edge 2728 3296 3947 4805 5165 potentialLowerBound 2727 3394 3821 4688 5140 potentialE 2692 2964 2957 3886 4380 (#vertices, edge multiplicity) (2,3) (2,4) (2,5) (3,2) (4,1) |E(G)| 3 4 5 6 6 $$\beta _{1}(G)$$ 2 3 4 4 3 |E(G)| ⋅ 924 2772 3698 4620 5544 5544 completion rate 65% 80% 90% 100% 100% Random Edge 2728 3296 3947 4805 5165 potentialLowerBound 2727 3394 3821 4688 5140 potentialE 2692 2964 2957 3886 4380 6.1 Sparse polynomial systems The example families in this subsection have the property that the support of the equations is fixed, while the coefficients can vary freely, as long as they are generic. We run the static graph strategy Algorithm 3.1 on these examples. Our timings do not include the $$\alpha$$-test, which was only applied in Section 6.1.2. 6.1.1 Cyclic roots The cyclic n-roots polynomial system is $$\begin{cases} i = 1, 2, 3, 4, \ldots, n-1: \displaystyle\sum_{j=0}^{n-1} ~ \prod_{k=j}^{j+i-1} x_{k~\textrm{mod}~n}=0 \\ x_{0}x_{1}x_{2} \cdots x_{n-1} - 1 = 0. \end{cases}$$ (6.1) This system is commonly used to benchmark polynomial system solvers. We will study the modified system with randomized coefficients and seek solutions in $$( {\mathbb{C}}\setminus \{0\})^{n}$$. Therefore, the solution count can be computed as the mixed volume of the Newton polytopes of the left-hand sides, providing a natural stopping criterion discussed in Section 3.2.2. This bound is 924 for cyclic 7. Tables 2 and 3 contain averages of experimental data from running 20 trials of Algorithm 3.1 on cyclic 7. The main measurement reported is the average number of paths tracked, as the unit of work for our algorithm is tracking a single homotopy path. The experiments were performed with 10 different graph layouts and three edge-selection strategies. With respect to number of paths tracked, we see that it is an advantage to keep the Betti number high and edge number low. Remark 6.1 Computing the expected success rates (> 99%) using Remark 4.4, we conclude that the resulting permutations do not conform to the model of picking uniformly from $$S_{924}$$. The completion rate depends on the choice of strategy (compare Table 2 to Table 3). Nevertheless, both in theory (assuming uniform distribution as in Section 4) and in practice (with distribution unknown to us), the completion rate does converge to 100% rapidly as the Betti number grows. 6.1.2 Nash equilibria Semimixed multihomogeneous systems arise when one is looking for all totally mixed Nash equilibria (TMNE) in game theory. A specialization of mixed volume using matrix permanents gives a concise formula for a root count for systems arising from TMNE problems (Emiris & Vidunas, 2014). We provide an overview of how such systems are constructed based on the study by Emiris & Vidunas (2014). Suppose there are N players with m options each. For player i ∈ {1, … , N} using option j ∈ {1, … , m} we have the equation $$P_{j}^{(i)} = 0$$, where \begin{align} \begin{aligned} P_{j}^{(i)} = \sum_{\substack{k_{1},\ldots,k_{i-1},\\ k_{i+1},\ldots, k_{N}}} a^{(i)}_{k_{1},\ldots, k_{i-1},\, j, k_{i+1}, \ldots, k_{N}}p_{k_{1}}^{(1)}p_{k_{2}}^{(2)}\cdots p_{k_{i-1}}^{(i-1)}p_{k_{i+1}}^{(i+1)} \cdots p_{k_{N}}^{(N)}. \end{aligned} \end{align} (6.2) The parameters $$a^{(i)}_{k_{1},k_{2},\ldots ,k_{N}}$$ are the pay-off rates for player i when players $$1,\ldots ,i-1,i+1,\ldots , N$$ are using options $$k_{1}, \ldots ,k_{i-1},k_{i+1},\ldots ,k_{N}$$, respectively. Here the unknowns are $$p^{(i)}_{k_{j}}$$, representing the probability that player i will use option $$k_{j} \in \{1,\ldots , m\}$$. There is one constraint on the probabilities for each player i ∈ {1, …, N}, namely the condition that \begin{align} p_{1}^{(i)}+p_{2}^{(i)}+\cdots+p_{m}^{(i)}=1. \end{align} (6.3) The system (4) consists of N ⋅ m equations in N ⋅ m unknowns. Using condition (5) reduces the number of unknowns to N(m − 1). Finally, we eliminate the $$P_{j}^{(i)}$$ by constructing \begin{align} P_{1}^{(i)}=P_{2}^{(i)}, ~ P_{1}^{(i)}=P_{3}^{(i)}, ~ \ldots ~, P_{1}^{(i)}=P_{m}^{(i)},\quad \textrm{ for each } i \in \{1,\ldots, N\}. \end{align} (6.4) The final system is a square system of N(m − 1) equations in N(m − 1) unknowns. For one of our examples (paper-examples/example-Nash.m2 at (Duff et al.)), we chose the generic system of this form for N = 3 players with m = 3 options for each. The result is a system of six equations in six unknowns and 81 parameters with 10 solutions. We also use this example to demonstrate that these solutions can be certified using NumericalCertification (Section 5.3). 6.2 Chemical reaction networks A family of interesting examples arises from chemical reaction network theory. A chemical reaction network considered under the laws of mass action kinetics leads to a dynamical polynomial system, the solutions of which represent all the equilibria for the given reaction network (MacLean et al., 2015; Gross et al., 2016). These polynomial systems are not generically sparse and we cannot easily compute their root count. In our experiments we used the stabilization stopping criterion, terminating the algorithm after a fixed number of iterations that do not deliver new points; the default is 10 fruitless iterations. Figure 4 gives an example of a small chemical reaction network. Fig. 4. View largeDownload slide Chemical reaction network example. Fig. 4. View largeDownload slide Chemical reaction network example. Applying the laws of mass action kinetics to the reaction network above we obtain the polynomial system (7) consisting of the corresponding steady-state and conservation equations. Here the $$k_{i}$$s represent the reaction rates, $$x_{i}$$s represent species concentrations and the $$c_{i}$$s are parameters. \begin{align} \begin{aligned} \dot{x_{A}} & = k_{1}{x_{B}^{2}}-k_{2}x_{A}-k_{3}x_{A}x_{C}+k_{4}x_{D}+k_{5}x_{B}x_{E} \\ \dot{x_{B}} & = 2k_{1}x_{A}-2k_{2}{x_{B}^{2}}+k_{4}x_{D}-k_{5}x_{B}x_{E} \\ \dot{x_{C}} & = -k_{3}x_{A}x_{C}+k_{4}x_{D}+k_{5}x_{B}x_{E} \\ \dot{x_{D}} & = k_{3}x_{A}x_{C}-(k_{4}+k_{6})x_{D} \\ \dot{x_{E}} & = -k_{5}x_{B}x_{E}+k_{6}x_{D} \\ 0 & = 2x_{A}+x_{B}-x_{C}+x_{D}-c_{1} \\ 0 & = -2x_{A}-x_{B}+2x_{C}+x_{E}-c_{2}. \end{aligned} \end{align} (6.5) Typically, systems resulting from chemical reaction networks will be overdetermined. With the current implementation one needs to either square the system or use a homotopy tracker that supports following a homotopy in a space of overdetermined systems. Table 4 Katsura-(n − 1) for the flower strategy (#vertices-1, edge multiplicity) n BKK Bound (3,2) (4,2) (5,2) (3,3) (4,3) $$\beta _{1}=3$$ $$\beta _{1}=4$$ $$\beta _{1}=5$$ $$\beta _{1}=6$$ $$\beta _{1}=8$$ 5 12 96.4% 99.4% 99.6% 100% 99.8% 6 30 98.6% 100% 99.8% 100% 99.6% 7 54 97.6% 98.8% 99.4% 99.4% 98.4% 8 126 99.2% 99.8% 99.6% 99.8% 99.8% 9 240 98.8% 99.6% 98.4% 98.4% 98.6% 10 504 98.6% 98.8% 99.2% 99.4% 98.8% (#vertices-1, edge multiplicity) n BKK Bound (3,2) (4,2) (5,2) (3,3) (4,3) $$\beta _{1}=3$$ $$\beta _{1}=4$$ $$\beta _{1}=5$$ $$\beta _{1}=6$$ $$\beta _{1}=8$$ 5 12 96.4% 99.4% 99.6% 100% 99.8% 6 30 98.6% 100% 99.8% 100% 99.6% 7 54 97.6% 98.8% 99.4% 99.4% 98.4% 8 126 99.2% 99.8% 99.6% 99.8% 99.8% 9 240 98.8% 99.6% 98.4% 98.4% 98.6% 10 504 98.6% 98.8% 99.2% 99.4% 98.8% Table 4 Katsura-(n − 1) for the flower strategy (#vertices-1, edge multiplicity) n BKK Bound (3,2) (4,2) (5,2) (3,3) (4,3) $$\beta _{1}=3$$ $$\beta _{1}=4$$ $$\beta _{1}=5$$ $$\beta _{1}=6$$ $$\beta _{1}=8$$ 5 12 96.4% 99.4% 99.6% 100% 99.8% 6 30 98.6% 100% 99.8% 100% 99.6% 7 54 97.6% 98.8% 99.4% 99.4% 98.4% 8 126 99.2% 99.8% 99.6% 99.8% 99.8% 9 240 98.8% 99.6% 98.4% 98.4% 98.6% 10 504 98.6% 98.8% 99.2% 99.4% 98.8% (#vertices-1, edge multiplicity) n BKK Bound (3,2) (4,2) (5,2) (3,3) (4,3) $$\beta _{1}=3$$ $$\beta _{1}=4$$ $$\beta _{1}=5$$ $$\beta _{1}=6$$ $$\beta _{1}=8$$ 5 12 96.4% 99.4% 99.6% 100% 99.8% 6 30 98.6% 100% 99.8% 100% 99.6% 7 54 97.6% 98.8% 99.4% 99.4% 98.4% 8 126 99.2% 99.8% 99.6% 99.8% 99.8% 9 240 98.8% 99.6% 98.4% 98.4% 98.6% 10 504 98.6% 98.8% 99.2% 99.4% 98.8% Although we may obtain large systems they typically have very low root counts compared to the sparse case. The polynomial system (7) has four solutions. A larger example is the wnt signaling pathway from Systems Biology (Gross et al., 2016) consisting of 19 polynomial equations with nine solutions. All nine solutions are obtained in less than a second with Algorithm 3.1. Table 5 Katsura-(n − 1) for the completeGraph strategy (#vertices, edge multiplicity) n BKK Bound (2,3) (2,4) (2,5) (3,2) (4,1) $$\beta _{1}=2$$ $$\beta _{1}=3$$ $$\beta _{1}=4$$ $$\beta _{1}=4$$ $$\beta _{1}=3$$ 5 12 65.6% 88.2% 95% 99.2% 98% 6 30 77.4% 95.2% 99% 99.8% 99.6% 7 54 74.4% 96.2% 99.2% 99.6% 99.8% 8 126 81.8% 97% 99.2% 100% 99.8% 9 240 85.2% 97.6% 99.4% 99% 98.2% 10 504 89.2% 98.2% 99.2% 99.4% 99% (#vertices, edge multiplicity) n BKK Bound (2,3) (2,4) (2,5) (3,2) (4,1) $$\beta _{1}=2$$ $$\beta _{1}=3$$ $$\beta _{1}=4$$ $$\beta _{1}=4$$ $$\beta _{1}=3$$ 5 12 65.6% 88.2% 95% 99.2% 98% 6 30 77.4% 95.2% 99% 99.8% 99.6% 7 54 74.4% 96.2% 99.2% 99.6% 99.8% 8 126 81.8% 97% 99.2% 100% 99.8% 9 240 85.2% 97.6% 99.4% 99% 98.2% 10 504 89.2% 98.2% 99.2% 99.4% 99% Table 5 Katsura-(n − 1) for the completeGraph strategy (#vertices, edge multiplicity) n BKK Bound (2,3) (2,4) (2,5) (3,2) (4,1) $$\beta _{1}=2$$ $$\beta _{1}=3$$ $$\beta _{1}=4$$ $$\beta _{1}=4$$ $$\beta _{1}=3$$ 5 12 65.6% 88.2% 95% 99.2% 98% 6 30 77.4% 95.2% 99% 99.8% 99.6% 7 54 74.4% 96.2% 99.2% 99.6% 99.8% 8 126 81.8% 97% 99.2% 100% 99.8% 9 240 85.2% 97.6% 99.4% 99% 98.2% 10 504 89.2% 98.2% 99.2% 99.4% 99% (#vertices, edge multiplicity) n BKK Bound (2,3) (2,4) (2,5) (3,2) (4,1) $$\beta _{1}=2$$ $$\beta _{1}=3$$ $$\beta _{1}=4$$ $$\beta _{1}=4$$ $$\beta _{1}=3$$ 5 12 65.6% 88.2% 95% 99.2% 98% 6 30 77.4% 95.2% 99% 99.8% 99.6% 7 54 74.4% 96.2% 99.2% 99.6% 99.8% 8 126 81.8% 97% 99.2% 100% 99.8% 9 240 85.2% 97.6% 99.4% 99% 98.2% 10 504 89.2% 98.2% 99.2% 99.4% 99% Table 6 Rounded expected probability of success assuming uniform distribuition of permutations and full monodromy group d $$\beta _{1}=2$$ $$\beta _{1}=3$$ $$\beta _{1}\geq 4$$ 12 90.5% 99.3% 100.0% 30 96.5% 99.9% 100.0% 54 98.1% 100.0% 100.0% 126 99.2% 100.0% 100.0% 240 99.6% 100.0% 100.0% 504 99.8% 100.0% 100.0% d $$\beta _{1}=2$$ $$\beta _{1}=3$$ $$\beta _{1}\geq 4$$ 12 90.5% 99.3% 100.0% 30 96.5% 99.9% 100.0% 54 98.1% 100.0% 100.0% 126 99.2% 100.0% 100.0% 240 99.6% 100.0% 100.0% 504 99.8% 100.0% 100.0% Table 6 Rounded expected probability of success assuming uniform distribuition of permutations and full monodromy group d $$\beta _{1}=2$$ $$\beta _{1}=3$$ $$\beta _{1}\geq 4$$ 12 90.5% 99.3% 100.0% 30 96.5% 99.9% 100.0% 54 98.1% 100.0% 100.0% 126 99.2% 100.0% 100.0% 240 99.6% 100.0% 100.0% 504 99.8% 100.0% 100.0% d $$\beta _{1}=2$$ $$\beta _{1}=3$$ $$\beta _{1}\geq 4$$ 12 90.5% 99.3% 100.0% 30 96.5% 99.9% 100.0% 54 98.1% 100.0% 100.0% 126 99.2% 100.0% 100.0% 240 99.6% 100.0% 100.0% 504 99.8% 100.0% 100.0% Table 7 Examples with solution count smaller than BKK bound (timings in seconds) problem wnt $${\mathop{SO}}(4)$$ $${\mathop{SO}}(5)$$ $${\mathop{SO}}(6)$$ $${\mathop{SO}}(7)$$ count 9 40 384 4768 111616 MonodromySolver 0.52 4 23 528 42791 Bertini 42 81 10605 out of memory PHCpack 862 103 > one day problem wnt $${\mathop{SO}}(4)$$ $${\mathop{SO}}(5)$$ $${\mathop{SO}}(6)$$ $${\mathop{SO}}(7)$$ count 9 40 384 4768 111616 MonodromySolver 0.52 4 23 528 42791 Bertini 42 81 10605 out of memory PHCpack 862 103 > one day Table 7 Examples with solution count smaller than BKK bound (timings in seconds) problem wnt $${\mathop{SO}}(4)$$ $${\mathop{SO}}(5)$$ $${\mathop{SO}}(6)$$ $${\mathop{SO}}(7)$$ count 9 40 384 4768 111616 MonodromySolver 0.52 4 23 528 42791 Bertini 42 81 10605 out of memory PHCpack 862 103 > one day problem wnt $${\mathop{SO}}(4)$$ $${\mathop{SO}}(5)$$ $${\mathop{SO}}(6)$$ $${\mathop{SO}}(7)$$ count 9 40 384 4768 111616 MonodromySolver 0.52 4 23 528 42791 Bertini 42 81 10605 out of memory PHCpack 862 103 > one day 6.3 Completion rate We investigate the completion rate of Algorithm 3.1 for the Katsura family parametrized by n with fixed support and generically chosen coefficients. Tables 4 and 5 contain the percentage of successes from 500 runs with distinct random seeds. In Table 6 we show the computed expected values using Remark 4.4. For $$\beta _{1}\geq 3$$ the observed success rates approach the expected values of Table 6. We note that the flower strategy is again closest to the estimates. We do not expect the numbers produced in experiments to match the numbers in Table 1, since the assumptions made for that statistical analysisare quite idealistic; however, both the analysis and experiments show that the probability of success approaches 100% rapidly as the number of solutions grows and the first Betti number increases. 6.4 Timings and comparison with other solvers All timings appearing in this section are done on one thread and on the same machine. Remarks 3.2 and 4.2 show that we should expect the number of tracked paths in Algorithms 3.1 and 3.4 to be linear (with a small constant!) in the number of solutions of the system. In this section we highlight the practicality of our approach in two ways. First, the monodromy method dramatically extends our computational ability for systems where the solution count turns out to be significantly smaller than the count corresponding to a more general family, for example, BKK count for sparse systems. This means that the existing blackbox methods, whose complexity relies on a larger count, are likely to spend significantly more time in computation compared to our approach. In Table 7 we collect timings on several challenging examples mentioned in recent literature where smaller solution counts are known, thus providing us with rigorous test cases for our heuristic stopping criterion. The first system in the table is that of the wnt signaling pathway reaction network mentioned in Section 6.2. The others come from the problem of computing the degree of $${\mathop{SO}}(n)$$, the special orthogonal group, as a variety (Brandt et al., 2017). Table 8 Software timings on large examples (in seconds) problem cyclic-10 cyclic-11 noon-10 BKK bound 35940 184756 59029 completeGraph(2,3) 610 7747 failed (107820 paths) (540155 paths) (59001 solutions) completeGraph(2,4) 740 8450 935 (129910 paths) (737432 paths) (236051 paths) PHCpack 538 4256 751 HOM4PS2 62 410 120 problem cyclic-10 cyclic-11 noon-10 BKK bound 35940 184756 59029 completeGraph(2,3) 610 7747 failed (107820 paths) (540155 paths) (59001 solutions) completeGraph(2,4) 740 8450 935 (129910 paths) (737432 paths) (236051 paths) PHCpack 538 4256 751 HOM4PS2 62 410 120 View Large Table 8 Software timings on large examples (in seconds) problem cyclic-10 cyclic-11 noon-10 BKK bound 35940 184756 59029 completeGraph(2,3) 610 7747 failed (107820 paths) (540155 paths) (59001 solutions) completeGraph(2,4) 740 8450 935 (129910 paths) (737432 paths) (236051 paths) PHCpack 538 4256 751 HOM4PS2 62 410 120 problem cyclic-10 cyclic-11 noon-10 BKK bound 35940 184756 59029 completeGraph(2,3) 610 7747 failed (107820 paths) (540155 paths) (59001 solutions) completeGraph(2,4) 740 8450 935 (129910 paths) (737432 paths) (236051 paths) PHCpack 538 4256 751 HOM4PS2 62 410 120 View Large Below is a list of comments on the set-up: For our implementation we chose small graphs with $$\beta _{1}\leq 4$$ and the random edge-selection strategy. The stopping criterion is ‘stabilization’ as discussed in Section 3.2.3. While the blackbox solver of PHCpack ultimately performs polyhedral homotopy continuation, Bertini relies by default on an equation-by-equation technique dubbed regeneration (see Bates et al., 2013). The latter may be faster than the former in certain cases, which this series of examples shows. Secondly, when the solution count is given by the BKK bound our method is a viable alternative to polyhedral homotopy solvers, since the number of paths we track is linear in the number of solutions. The timings on a few large benchmark problems of our current implementation and several other software packages are in Table 8. Our goal in the rest of this section is to show that our running times are in the same ballpark as polyhedral homotopies. Below is a list of comments on the set-up: For our implementation we chose two small graphs and default (random) edge-selection strategy. For PHCpack there is a way to launch a mixed volume computation, with the option of creating a system with the same support and random coefficients together with its solutions. This is the option we are using; the blackbox computation takes a little longer. HOM4PS2 (Lee et al., 2008) is not open source, unlike all other software mentioned here. (We use HOM4PS2 stock examples for all systems and call its blackbox polyhedral homotopies solver.) HOM4PS2 may use just-in-time compilation of straight-line programs used for evaluation, which speeds up computations considerably. (PHCpack does not use this technique; neither does our software, but our preliminary experiments in Macaulay2 show a potential for a 10- to 20-fold speed up over our currently reported timings.) Remark 6.2 For large examples assuming the probabilistic model leading to Theorem 4.1 and Remark 4.2, the probability of success should be extremely close to 100% even for a random graph with $$\beta _{1} = 2$$. The run of noon-10, which is an example of neural network model from the study by Noonburg (1989), demonstrates an unlikely, but possible failure for $$\beta _{1} = 2$$ followed by success at $$\beta _{1} = 3$$. On the examples in Table 8 we also ran the blackbox solvers of Bertini and Numerical Algebraic Geometry (Leykin, 2011), which use the total-degree homotopy. Both were able to finish noon-10 with timings similar to the table, but all other problems took longer than a day. This is expected, as the BKK bound of noon-10 is only slightly sharper than the Bézout bound. Remark 6.3 In comparison with the naive dynamic strategy (Section 3.1) our framework loses slightly only in one aspect: memory consumption. For a problem with d solutions the naive approach stores up to (and typically close to) 2d points. The number of points our approach stores is up to (and typically considerably fewer than) d times the number of vertices. For instance, it is up to 4d points in all runs in Table 8. The number of tracked paths is significantly lower in our framework: for example, the naive strategy tracks about 7500 paths on average for cyclic 7. Even before looking at Table 2 it is clear that running the flower strategy in combination with the incremental dynamic strategy of Section 3.3 guarantees to dominate the naive strategy. 7. Generalizations While we propose a more general algorithmic framework, a concurrent goal of this paper is to demonstrate that significant practical advantages are already apparent when we apply a relatively simple implementation and analysis to simple problems (linearly parametrized families). The following topics thus lie outside the scope of this article, but seem deserving of further study: One advantage of the MS approach is that it can tolerate numerical failures of the underlying homotopy tracker. In fact, we already implemented a simple failure resistant mechanism, and it successfully tolerates a few failures that arise in some runs for large test examples in Section 6.4. A natural extension of this paper’s statistical analysis would be to model the algorithm’s performance in the presence of failures. Ideally, heuristics such as edge potentials should incorporate information such as the failures discussed above. It is also of interest to adapt potentials to the parallel setting discussed below. The parallelization of the MS approach is not as straightforward as that of other homotopy continuation methods. The question of when speed-ups close to linear can be achieved should be addressed. Consider the generalized set-up in which the base space B is an irreducible variety and the family is given by a rational map from P into a space of systems. To apply our general framework, a major requirement is to find an effective way to parametrize a curve between two points of P. This parametrization would conceivably depend on the nature of the problem being considered. Certain other ingredients are also likely to be problem-specific—for instance, even in the case of $$P = {\mathbb{C}}^{m},$$ the construction of the initial seed $$(p_{0},x_{0})$$ is complicated by the possibility that the systems’ coefficients are nonlinear in the parameters. Nonetheless, this is one of the strengths of the MS framework—once all required ‘oracles’ are supplied, the procedures become effective. In the classical language of enumerative geometry, the monodromy groups we consider are isomorphic to Galois groups of incidence varieties (essentially solution varieties in our terminology). For a large class of Schubert problems and other interesting incidence varieties the associated Galois group turns out to be the full symmetric group (Leykin & Sottile, 2009). A suitable modification of our dynamic strategy is one practical approach to verifying this in conjectural cases. Our paper demonstrates the strength of our method relative to other techniques such as polyhedral homotopy and regeneration. Building on our framework one could use polyhedral homotopy as a subroutine to quickly populate a partial solution set (quickly discarding any path that becomes poorly conditioned). Further advantages may be achievable by using different techniques in parallel. These and other hybrid approaches have the potential to produce even faster and more robust blackbox solvers. Funding Research of T.D., C.H., K.L., and A.L. is supported in part by National Science Foundation grants DMS-1151297 and DMS-1719968. References Babai , L. ( 1989 ) The probability of generating the symmetric group . J. Comb. Theory Ser. A , 52 , 148 – 153 . Google Scholar CrossRef Search ADS Bates , D. J. , Hauenstein , J. D. , Sommese , A. J. & Wampler , C. W. ( 2013 ) Numerically Solving Polynomial Systems With Bertini , vol. 25. Philadelphia : SIAM . Bernstein , D. N. ( 1975 ) The number of roots of a system of equations . Funk. Anal. Priložen , 9 , 1 – 4 . Blum , L. , Cucker , F. , Shub , M. & Smale , S. ( 1998 ) Complexity and Real Computation . New York : Springer . Google Scholar CrossRef Search ADS Brandt , M. , Bruce , J. , Brysiewicz , T. , Krone , R. & Robeva , E. ( 2017 ) The degree of SO($$n,\mathbb{C}$$) . Combinatorial Algebraic Geometry (G. Smith & B. Sturmfels eds). Fields Institute Communications, vol. 80 . New York City : Springer , pp. 229 – 246 . Chen , J. & Kileel , J. ( 2016 ) Numerical implicitization for Macaulay2 . arXiv preprint arXiv:1610.03034 . del Campo , A. M. & Rodriguez , J. I. ( 2017 ) Critical points via monodromy and local methods . J. Symbolic Comput. , 79 , 559 – 574 . Google Scholar CrossRef Search ADS Dixon , J. D. ( 1969 ) The probability of generating the symmetric group . Math. Z. , 110 , 199 – 205 . Google Scholar CrossRef Search ADS Duff , T. , Hill , C. , Jensen , A. , Lee , K. , Leykin , A. & Sommars , J. ( 2016 ) MonodromySolver: a Macaulay2 package for solving polynomial systems via homotopy continuation and monodromy . Available at http://people.math.gatech.edu/aleykin3/MonodromySolver. Emiris , I. Z. & Vidunas , R. ( 2014 ) Root counts of semi-mixed systems, and an application to counting nash equilibria . Proceedings of the 39th International Symposium on Symbolic and Algebraic Computation, ISSAC ’14 . New York, NY, USA : Association for Computing Machinery , pp. 154 – 161 . Galligo , A. & Miclo , L. ( 2012 ) On the cut-off phenomenon for the transitivity of randomly generated subgroups . Random Struct. Alg. , 40 , 182 – 219 . Google Scholar CrossRef Search ADS Galligo , A. & Poteaux , A. ( 2011 ) Computing monodromy via continuation methods on random Riemann surfaces . Theor. Comput. Sci. , 412 , 1492 – 1507 . Google Scholar CrossRef Search ADS Grayson , D. R. & Stillman , M. E. Macaulay2, a software system for research in algebraic geometry . Available at http://www.math.uiuc.edu/Macaulay2/. Gross , E. , Harrington , H. A. , Rosen , Z. & Sturmfels , B. ( 2016 ) Algebraic systems biology: a case study for the Wnt pathway . Bull. Math. Biol. , 78 , 21 – 51 . Google Scholar CrossRef Search ADS PubMed Hauenstein , J. D. & Rodriguez , J. I. ( 2015 ) Multiprojective witness sets and a trace test . arXiv preprint arXiv:1507.07069 . Hauenstein , J. D. , Rodriguez , J. I. & Sottile , F. ( 2017 ) Numerical computation of Galois groups . Found. Comput. Math . To appear . Hauenstein , J. D. & Sottile , F. ( 2012 ) Algorithm 921: alphaCertified: certifying solutions to polynomial systems . ACM Trans. Math. Softw. , 38 , 28 . Google Scholar CrossRef Search ADS Huber , B. & Sturmfels , B. ( 1995 ) A polyhedral method for solving sparse polynomial systems . Math. Comp. , 64 , 1541 – 1555 . Google Scholar CrossRef Search ADS Jensen , A. N. ( 2016 ) An implementation of exact mixed volume computation . Mathematical Software - ICMS 2016 - 5th International Conference, Berlin, Germany, July 11–14, 2016, Proceedings . Cham, Switzerland : Springer , pp. 198 – 205 . Lee , T.-L. , Li , T.-Y. & Tsai , C.-H. ( 2008 ) HOM4PS-2.0: a software package for solving polynomial systems by the polyhedral homotopy continuation method . Computing , 83 , 109 – 133 . Google Scholar CrossRef Search ADS Leykin , A. ( 2011 ) Numerical algebraic geometry . J. Softw. Alg. Geom. , 3 , 5 – 10 . Google Scholar CrossRef Search ADS Leykin , A. , Rodriguez , J. I. & Sottile , F. ( 2016 ) Trace test . arXiv preprint arXiv:1608.00540 . Leykin , A. & Sottile , F. ( 2009 ) Galois groups of Schubert problems via homotopy computation . Math. Comp. , 78 , 1749 – 1765 . Google Scholar CrossRef Search ADS Leykin , A. & Verschelde , J. ( 2009 ) Decomposing solution sets of polynomial systems: a new parallel monodromy breakup algorithm . Int. J. Comput. Sci. Eng. , 4 , 94 – 101 . Google Scholar CrossRef Search ADS MacLean , A. L. , Rosen , Z. , Byrne , H. M. & Harrington , H. A. ( 2015 ) Parameter-free methods distinguish Wnt pathway models and guide design of experiment . Proc. Natl. Acad. Sci. , 112 , 2652 – 2657 . Google Scholar CrossRef Search ADS Malajovich , G. ( 2017 ) Computing mixed volume and all mixed cells in quermassintegral time . Found. Comput. Math. , 17 , 1293 – 1334 . Google Scholar CrossRef Search ADS Morgan , A. ( 1987 ) Solving Polynomial Systems Using Continuation for Engineering and Scientific Problems . Englewood Cliffs, NJ : Prentice Hall Inc . Noonburg , V. W. ( 1989 ) A neural network modeled by an adaptive Lotka–Volterra system . SIAM J. Appl. Math. , 49 , 1779 – 1792 . Google Scholar CrossRef Search ADS Sommese , A. J. , Verschelde , J. & Wampler , C. W. ( 2001 ) Using monodromy to decompose solution sets of polynomial systems into irreducible components . Application of Algebraic Geometry to Coding Theory, Physics, and Computation , Dordrecht : Springer Netherlands , pp. 297 – 315 . Google Scholar CrossRef Search ADS Sommese , A. J. , Verschelde , J. & Wampler , C. W. ( 2002 ) Symmetric functions applied to decomposing solution sets of polynomial systems . SIAM J. Numer. Anal. , 40 , 2026 – 2046 . Google Scholar CrossRef Search ADS Sommese , A. J. , Verschelde , J. & Wampler , C. W. ( 2005 ) Introduction to numerical algebraic geometry . Solving Polynomial Equations: Foundations, Algorithms, and Applications (A. Dickenstein & I. Z. Emiris eds). Berlin, Heidelberg : Springer Berlin Heidelberg , pp. 301 – 337 . Sommese , A. J. & Wampler II , C. W. ( 2005 ) The Numerical Solution of Systems of Polynomials . Hackensack, NJ : World Scientific Publishing Co. Pte. Ltd . Google Scholar CrossRef Search ADS Verschelde , J. ( 1999 ) Algorithm 795: PHCpack: a general-purpose solver for polynomial systems by homotopy continuation . ACM Trans. Math. Softw. , 25 , 251 – 276 . Google Scholar CrossRef Search ADS Verschelde , J. , Verlinden , P. & Cools , R. ( 1994 ) Homotopies exploiting Newton polytopes for solving sparse polynomial systems . SIAM J. Numer. Anal. , 31 , 915 – 930 . Google Scholar CrossRef Search ADS © The Author(s) 2018. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices) For permissions, please e-mail: journals. permissions@oup.com http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png IMA Journal of Numerical Analysis Oxford University Press

# Solving polynomial systems via homotopy continuation and monodromy

, Volume Advance Article – Apr 13, 2018
26 pages

/lp/ou_press/solving-polynomial-systems-via-homotopy-continuation-and-monodromy-iSVXMp1ly9
Publisher
Oxford University Press
ISSN
0272-4979
eISSN
1464-3642
D.O.I.
10.1093/imanum/dry017
Publisher site
See Article on Publisher Site

### Abstract

Abstract We study methods for finding the solution set of a generic system in a family of polynomial systems with parametric coefficients. We present a framework for describing monodromy-based solvers in terms of decorated graphs. Under the theoretical . that monodromy actions are generated uniformly, we show that the expected number of homotopy paths tracked by an algorithm following this framework is linear in the number of solutions. We demonstrate that our software implementation is competitive with the existing state-of-the-art methods implemented in other software packages. 1. Introduction Homotopy continuation has become a standard technique to find approximations of solutions of polynomial systems. There is an early popular text on the subject and its applications in the study by Morgan (1987). This technique is the backbone of Numerical Algebraic Geometry, the area that classically addresses the questions of complex algebraic geometry through algorithms that employ numerical approximate computation. The chapter in the study by Sommese et al. (2005, Section 8) is the earliest introduction and the book by Sommese & Wampler (2005) is the primary reference in the area. Families of polynomial systems with parametric coefficients play one of the central roles. Most homotopy continuation techniques could be viewed as going from a generic system in the family to a particular one. This process is commonly referred to as degeneration. Going in the reverse direction it may be called deformation, undegeneration or regeneration depending on the literature. Knowing the solutions of a generic system, one can use coefficient-parameter homotopy (Sommese & Wampler, 2005, Section 7) to get to the solution of a particular one. The main problem that we address here is how to solve a generic system in a family of systems, $$F_{p} = \left(f_{p}^{(1)},\ldots,\,f_{p}^{(N)}\right) = 0, \quad f_{p}^{(i)}\in{\mathbb{C}}[p][x],\ i=1,\ldots,N,$$ with finitely many parameters p and n variables x. In the main body of the paper we restrict our attention to linear parametric families of systems, defined as systems with affine linear parametric coefficients, such that for a generic p we have a nonempty finite set of solutions x to $$F_{p}(x)=0$$. This implies N ≥ n. The number of parameters is arbitrary, but we require that for a generic x there exists p with $$F_{p}(x)=0$$. These restrictions are made for the sake of simplicity. We explain what modifications are needed to apply our approach in more general settings in Section 7. Linear parametric systems form a large class that includes sparse polynomial systems. These are square (n = N) systems with a fixed monomial support for each equation and a distinct parameter for the coefficient of each monomial. Polyhedral homotopy methods for solving sparse systems stem from the BKK (Bernstein, Khovanskii, Kouchnirenko) bound on the number of solutions (Bernstein, 1975); the early work on algorithm development was done in the studies by Verschelde et al. (1994) and Huber & Sturmfels (1995). Polyhedral homotopies provide an optimal solution to sparse systems in the sense that they are designed to follow exactly as many paths as the number of solutions of a generic system (the BKK bound). The method that we propose is clearly not optimal in the above sense. The expected number of homotopy paths followed can be larger than the number of solutions, though not significantly larger. We also use linear segment homotopies that are significantly simpler and less expensive to follow in practice. Our current implementation shows it is competitive with the state-of-the-art implementations of polyhedral homotopies in PHCpack (Verschelde, 1999) and HOM4PS2 (Lee et al., 2008) for solving sparse systems. In a setting more general than sparse, we demonstrate examples of linear parametric systems for which our implementation exceeds the capabilities of the existing sparse system solvers and blackbox solvers based on other ideas. The idea of using the monodromy action induced by the fundamental group of the regular locus of the parameter space has been successfully employed throughout Numerical Algebraic Geometry. One of the main tools in the area, numerical irreducible decomposition, can be efficiently implemented using the monodromy breakup algorithm, which first appeared in the study by Sommese et al. (2001). One parallel incarnation of the monodromy breakup algorithm is described in the study by Leykin & Verschelde (2009). In fact, the main idea in that work is close in spirit to what we propose in this article. The idea to use monodromy to find solutions drives numerical implicitization (Chen & Kileel, 2016) and appears in other works such as del Campo & Rodriguez (2017). Computing monodromy groups numerically, as in the studies by Leykin & Sottile (2009) and Hauenstein et al. (2017), requires more computation than just finding solutions. One can approach this computation with the same methodology as we propose; see 5 of Section 7. Our main contribution is a new framework to describe algorithms for solving polynomial systems using monodromy; we call it the Monodromy Solver (MS) framework. We analyze the complexity of our main algorithm both theoretically assuming a certain statistical model and experimentally on families of examples. The analysis gives us grounds to say that the expected number of paths tracked by our method is linear, with a small coefficient, as the number of solutions grows. Our method and its implementation not only provide a new general tool for solving polynomial systems, but also can solve some problems out of reach for other existing software. The structure of the paper is as follows. We give a brief overview of the MS method intermingled with some necessary preliminaries in Section 2. An algorithm following the MS framework depends on a choice of strategy with several possibilities outlined in Section 3. Statistical analysis of the method is the topic of Section 4. The implementation is discussed in Section 5 together with the side topic of certification of the solution set. The results of our experiments on selected example families highlighting various practical computational aspects are in Section 6. The reader may also want look at examples of systems in Sections 6.1 and 6.2 before reading some earlier sections. Possible generalizations of the MS technique and the future directions to explore are presented in Section 7. 2. Background and framework overview Let $$m, n\in {\mathbb{N}}$$. We consider the complex linear space of square systems $$F_{p}$$, $$p\in {\mathbb{C}}^{m}$$, where the monomial support of $$f_{p}^{(1)},\dots ,\,f_{p}^{(n)}$$ in the variables $$x=(x_{1},\dots ,x_{n})$$ is fixed and the coefficients vary. By a base space B we mean a parametrized linear variety of systems. We think of it as the image of an affine linear map $$\varphi : p\mapsto F_{p}$$ from a parameter space $${\mathbb{C}}^{m}$$ with coordinates $$p=(p_{1},\ldots ,p_{m})$$ to the space of systems. We assume the structure of our family is such that the projection $$\pi$$ from the solution variety $$V = \left\{\left(F_{p},x\right) \in B\times{\mathbb{C}}^{n} \mid F_{p}(x)=0\right\}$$ to B gives us a branched covering, i.e., the fibre $$\pi ^{-1}\left (F_{p}\right )$$ is finite of the same cardinality for a generic p. The discriminant variety D in this context is the subset of the systems in the base space with nongeneric fibres; it is also known as the branch locus of $$\pi$$. The fundamental group$$\pi _{1}(B\setminus D)$$—note that $$\pi _{1}$$ is a usual topological notation that is not related to the map $$\pi$$ above—as a set consists of loops, i.e., paths in B ∖ D starting and finishing at a fixed p ∈ B ∖ D considered up to homotopy equivalence. The definition, more details to which one can find in Section 2.1, does not depend on the point p, since B ∖ D is connected. Each loop induces a permutation of the fibre $$\pi ^{-1}\left (F_{p}\right )$$, which is referred to as a monodromy action. Our goal is to find the fibre of one generic system in our family. Our method is to find one pair $$(p_{0},x_{0})\in V$$ and use the monodromy action on the fibre $$\pi ^{-1}\left (F_{p_{0}}\right )$$ to find its points. We assume that this action is transitive, which is the case if and only if the solution variety V is irreducible. If V happens to be reducible we replace V with its unique dominant irreducible component as explained in Remark 2.2. 2.1 Monodromy We briefly review the basic facts concerning monodromy groups of branched coverings. With notation as before fix a system $$F_{p}\in B \setminus D$$ and consider a loop $$\tau$$ without branch points based at $$F_{p}$$; that is, a continuous path $$\tau : [0,1] \rightarrow B \setminus D$$ such that $$\tau (0) = \tau (1) = F_{p}.$$ Suppose we are also given a point $$x_{i}$$ in the fibre $$\pi ^{-1} \left (F_{p}\right )$$ with d points $$x_{1}, x_{2}, \ldots , x_{d}.$$ Since $$\pi$$ is a covering map the pair $$(\tau , x_{i})$$ corresponds to a unique lifting$$\widetilde{\tau _{i}},$$ a path $$\widetilde{\tau_{i}} : [0,1] \rightarrow V$$ such that $$\widetilde{\tau _{i}} (0) = x_{i}$$ and $$\widetilde{\tau _{i}} (1) = x_{j}$$ for some 1 ≤ j ≤ d. Note that the reversal of $$\tau$$ and $$x_{j}$$ lifts to a reversal of $$\widetilde{\tau _{i}}$$. Thus, the loop $$\tau$$ induces a permutation of the set $$\pi ^{-1} \left (F_{p}\right ).$$ We have a group homomorphism $$\varphi : \pi_{1} \left(B\setminus D, F_{p}\right) \rightarrow{\mathrm{S}}_{d}$$ whose domain is the usual fundamental group of B ∖ D based at $$F_{p}$$. The image of $$\varphi$$ is the monodromy group associated to $$\pi ^{-1} \left (F_{p}\right ).$$ The monodromy group acts on the fibre $$\pi ^{-1} \left (F_{p}\right )$$ by permuting the solutions of $$F_{p}$$. Remark 2.1 A reader familiar with the notion of a monodromy loop in the discussion of Sommese & Wampler (2005, Section 15.4) may think of this keyword referring to a representative of an element of the fundamental group together with its liftings to the solution variety and the induced action on the fibre. For the purposes of this article we need to be clear about the ingredients bundled in this term. We have not used any algebraic properties so far. The construction of the monodromy group above holds for an arbitrary covering with finitely many sheets. The monodromy group is a transitive subgroup of $${\mathrm{S}}_{d}$$ whenever the total space is connected. In our setting since we are working over $${\mathbb{C}}$$, this occurs precisely when the solution variety is irreducible. Remark 2.2 For a linear family we can show that there is at most one irreducible component of the solution variety V for which the restriction of the projection $$\left (F_{p},x\right )\mapsto x$$ is dominant (that is, its image is dense). We call such component the dominant component. Indeed, let U be the locus of points $$\left (F_{p},x\right )\in \pi ^{-1}(B\setminus D)$$ such that the restriction of the x-projection map is locally surjective, and the solution to the linear system of equations $$F_{p}(x)=0$$ in p has the generic dimension. Being locally surjective could be interpreted either in the sense of Zariski topology or as inducing surjection on the tangent spaces. Then either U is empty or $$\overline{U}$$ is the dominant component we need, since it is a vector bundle over an irreducible variety, and is hence irreducible. In the rest of the paper when we say solution variety, we mean the dominant component of the solution variety. In particular, for sparse systems restricting the attention to the dominant component translates into looking for solutions only in the torus $$\left ( {\mathbb{C}}^{\ast }\right )^{n}$$. 2.2 Homotopy continuation Given two points $$F_{p_{1}}$$ and $$F_{p_{2}}$$ in the base space B, we may form the family of systems $$H(t) = (1-t)F_{p_{1}} + t F_{p_{2}}\,,\quad t\in[0,1],$$ known as the linear segment homotopy between the two systems. If $$p_{1}$$ and $$p_{2}$$ are sufficiently generic for each t ∈ [0, 1], we have H(t) outside the real codimension 2 set D. Consequently, each system H(t) has a finite and equal number of solutions. This homotopy is a path in B; a lifting of this path in the solution variety V is called a homotopy path. The homotopy paths of H(t) establish a one-to-one correspondence between the fibres $$\pi ^{-1}\left (F_{p_{1}}\right )$$ and $$\pi ^{-1}\left (F_{p_{2}}\right )$$. Remark 2.3 Note that $$\gamma F_{p}$$ for $$\gamma \in {\mathbb{C}}\setminus \{0\}$$ has the same solutions as $$F_{p}$$. Let us scale both ends of the homotopy by taking a homotopy between $$\gamma _{1}F_{p_{1}}$$ and $$\gamma _{2}F_{p_{2}}$$ for generic $$\gamma _{1}$$ and $$\gamma _{2}$$. If the coefficients of $$F_{p}$$ are homogeneous in p then $$H^{\prime}(t) = (1-t)\gamma_{1}F_{p_{1}} + t \gamma_{2}F_{p_{2}} = F_{(1-t)\gamma_{1}p_{1}+t\gamma_{2}p_{2}}\,, \quad t\in[0,1],$$ is a homotopy matching solutions $$\pi ^{-1}\left (F_{p_{1}}\right )$$ and $$\pi ^{-1}\left (F_{p_{2}}\right )$$, where the matching is potentially different from that given by H(t). Similarly, for an affine linear family, $$F_{p} = F^{\prime}_{p} + C$$, where $$F^{\prime}_{p}$$ is homogeneous in p and C is a constant system, we have $$H^{\prime}(t) = (1-t)\gamma_{1}F_{p_{1}} + t \gamma_{2}F_{p_{2}} = F^{\prime}_{(1-t)\gamma_{1}p_{1}+t\gamma_{2}p_{2}}+\left((1-t)\gamma_{1}+t\gamma_{2}\right)C.$$ We ignore the fact that H′(t) may go outside B for t ∈ (0, 1), since its rescaling, \begin{align*} H^{\prime\prime}(t) &= \frac{1}{(1-t)\gamma_{1}+t\gamma_{2}}H^{\prime}(t)\\ &=F^{\prime}_{\frac{(1-t)\gamma_{1}p_{1}+t\gamma_{2}p_{2}}{(1-t)\gamma_{1}+t\gamma_{2}}}+C = F_{\frac{(1-t)\gamma_{1}p_{1}+t\gamma_{2}p_{2}}{(1-t)\gamma_{1}+t\gamma_{2}}}\,, \quad t\in[0,1], \end{align*} does not leave B and clearly has the same homotopy paths. Note that H″(t) is well defined as $$(1-t)\gamma _{1}+t\gamma _{2}\neq 0$$ for all t ∈ [0, 1] for generic $$\gamma _{1}$$ and $$\gamma _{2}$$. One may use methods of numerical homotopy continuation, described, for instance in the book by Sommese & Wampler (2005, Section 2.3), to track the solutions as t changes from 0 to 1. In some situations the path in B may pass close to the branch locus D and numerical issues must be considered. Remark 2.4 If the family $$F_{p}$$ is nonlinear in the parameters p one has to take the parameter linear segment homotopy in the parameter space, i.e., $$H(t) = F_{(1-t)p_{1}+tp_{2}}$$, t ∈ [0, 1]. This does not change the overall construction; however, the freedom to replace the systems $$F_{p_{1}}$$ and $$F_{p_{2}}$$ at the ends of the homotopy with their scalar multiples as in Remark 2.3 is lost. 2.3 Graph of homotopies: main ideas Some readers may find it helpful to use the examples of Section 2.4 for graphical intuition as we introduce notation and definitions below. To organize the discovery of new solutions we represent the set of homotopies by a finite undirected graph G. Let E(G) and V(G) denote the edge and vertex set of G, respectively. Any vertex v in V(G) is associated to a point $$F_{p}$$ in the base space. An edge e in E(G) connecting $$v_{1}$$ and $$v_{2}$$ in V(G) is decorated with two complex numbers, $$\gamma _{1}$$ and $$\gamma _{2}$$, and represents the linear homotopy connecting $$\gamma _{1}F_{p_{1}}$$ and $$\gamma _{2} F_{p_{2}}$$ along a line segment (Remark 2.3). We assume that both $$p_{i}$$ and $$\gamma _{i}$$ are chosen so that the segments do not intersect the branch locus. Choosing these at random (see Section 5.1 for a possible choice of distribution) satisfies the assumption, since the exceptional set of choices where such intersections happen is contained in a real Zariski closed set, see Sommese & Wampler (2005, Lemma 7.1.3). We allow multiple edges between two distinct vertices but no loops, since the latter induce trivial homotopies. For a graph G to be potentially useful in a monodromy computation it must contain a cycle. Some of the general ideas behind the structure of a graph G are listed below. For each vertex $$v_{i}$$ we maintain a subset of known points $$Q_{i}\subset \pi ^{-1}\left (F_{p_{i}}\right )$$. For each edge e between $$v_{i}$$ and $$v_{j}$$ we record the two complex numbers $$\gamma _{1}$$ and $$\gamma _{2}$$, and we store the known partial correspondences $$C_{e}\!\subset\! \pi ^{-\!1}\!\left (F_{p_{i}}\right )\!\times \pi ^{-\!1}\!\left (F_{p_{j}}\right )$$ between known points $$Q_{i}$$ and $$Q_{j}$$. At each iteration we pick an edge and direction, track the corresponding homotopy starting with yet unmatched points, and update known points and correspondences between them. We may obtain the initial ‘knowledge’ as a seed pair$$(p_{0},x_{0})$$ by picking $$x_{0}\in {\mathbb{C}}^{n}$$ at random and choosing $$p_{0}$$ to be a generic solution of the linear system $$F_{p}(x_{0})=0$$. We list basic operations that result in transition between one state of our algorithm captured by G, $$Q_{i}$$ for $$v_{i}\in V(G)$$ and $$C_{e}$$ for e ∈ E(G) to another. For an edge $$e = v_{i}\overset{(\gamma _{1},\gamma _{2})}{\overleftrightarrow{\hphantom{(\gamma _{1},\gamma _{2})}}} v_{j}$$ consider the homotopy $$H^{(e)} = (1-t)\gamma_{1}F_{p_{i}} + t \gamma_{2}F_{p_{j}}$$ where $$\left (\gamma _{1},\gamma _{2}\right )\in {\mathbb{C}}^{2}$$ is the label of e. Take start points $$S_{i}$$ to be a subset of the set of known points $$Q_{i}$$ that do not have an established correspondence with points in $$Q_{j}$$. Track $$S_{i}$$ along $$H^{(e)}$$ for t ∈ [0, 1] to get $$S_{j} \subset \pi ^{-1}\left (F_{p_{j}}\right )$$. Extend the known points for $$v_{j}$$, that is, $$Q_{j} := Q_{j} \cup S_{j}$$ and record the newly established correspondences. Add a new vertex corresponding to $$F_{p}$$ for a generic p ∈ B ∖ D. Add a new edge $$e = v_{i}\overset{(\gamma _{1},\gamma _{2})}{\overleftrightarrow{\hphantom{(\gamma _{1},\gamma _{2})}}}v_{j}$$ between two existing vertices decorated with generic $$\gamma _{1},\gamma _{2}\!\in\! {\mathbb{C}}$$. At this point a reader who is ready to see an algorithm based on these ideas may skip to Algorithm 3.1. 2.4 Graph of homotopies: examples We demonstrate the idea of graphs of homotopies, the core idea of the MS framework, by giving two examples. Example 2.5 Figure 1 shows a graph G with two vertices and three edges embedded in the base space B with paths partially lifted to the solution variety, which is a covering space with three sheets. The two fibres $$\{x_{1},x_{2},x_{3}\}$$ and $$\{y_{1},y_{2},y_{3}\}$$ are connected by three partial correspondences induced by the liftings of three egde paths. Note that several aspects in this illustration are fictional. There is only one branch point in the actual complex base space B that we would like the reader to imagine. The visible self-intersections of the solution variety V are an artifact of drawing the picture in the real space. Also, in practice we use homotopy paths as simple as possible, however, here the paths are more involved for the purpose of distinguishing them in print. An algorithm that we envision may hypothetically take the following steps: (1) seed the first fibre with $$x_{1}$$; (2) use a lifting of edge $$e_{a}$$ to get $$y_{1}$$ from $$x_{1}$$; (3) use a lifting of edge $$e_{b}$$ to get $$x_{2}$$ from $$y_{1}$$; (4) use a lifting of edge $$e_{c}$$ to get $$y_{2}$$ from $$x_{1}$$; (5) use a lifting of edge $$e_{a}$$ to get $$x_{3}$$ from $$y_{2}$$. Note that it is not necessary to complete the correspondences (a), (b) and (c). Doing so would require tracking nine continuation paths, while the hypothetical run above uses only four paths to find a fibre. Example 2.6 Figure 2 illustrates two partial correspondences associated to two edges $$e_{a}$$ and $$e_{b}$$, both connecting two vertices $$v_{1}$$ and $$v_{2}$$ in V(G). Each vertex $$v_{i}$$ stores the array of known points $$Q_{i}$$, which are depicted in solid. Both correspondences in the picture are subsets of a perfect matching, a one-to-one correspondence established by a homotopy associated to the edge. Note that taking the set of start points $$S_{1}=\{x_{3}\}$$ and following the homotopy $$H^{(e_{a})}$$ from left to right is guaranteed to discover a new point in the second fibre. On the other hand, it is impossible to obtain new knowledge by tracking $$H^{(e_{a})}$$ from right to left. Homotopy $$H^{(e_{b})}$$ has a potential to discover new points if tracked in either direction. We can choose $$S_{1}=\{x_{1},x_{3}\}$$ as the start points for one direction and $$S_{2}=\{y_{3}\}$$ for the other. In this scenario, following the homotopy from left to right is guaranteed to produce at least one new point, while going the other way may either deliver a new point or just augment the correspondences between the already known points. If the correspondences in (a) and (b) are completed to one-to-one correspondences of the fibres taking the homotopy induced by the edge $$e_{a}$$ from left to right followed by the homotopy induced by edge $$e_{b}$$ from right to left would produce a permutation. However, the group generated by this permutation has to stabilize $$\{x_{2}\}$$, therefore, it would not act transitively on the fibre of $$v_{1}$$. One could also imagine a completion such that the given edges would not be sufficient to discover $$x_{5}$$ and $$y_{4}$$. Fig. 1. View largeDownload slide Selected liftings of three edges connecting the fibres of two vertices and induced correspondences. Fig. 1. View largeDownload slide Selected liftings of three edges connecting the fibres of two vertices and induced correspondences. Fig. 2. View largeDownload slide Two partial correspondences induced by edges $$e_{a}$$ and $$e_{b}$$ for the fibres of the covering map of degree d = 5 in Example 2.6. Fig. 2. View largeDownload slide Two partial correspondences induced by edges $$e_{a}$$ and $$e_{b}$$ for the fibres of the covering map of degree d = 5 in Example 2.6. Fig. 3. View largeDownload slide Graphs for the flower(4,2) strategy and completeGraph(5,1). Fig. 3. View largeDownload slide Graphs for the flower(4,2) strategy and completeGraph(5,1). In our algorithm we record and use correspondences; however, they are viewed as a secondary kind of knowledge. In particular, in Section 3.2.4 we develop heuristics driven by edge potential functions which look to maximize the number of newly discovered solutions, in other words, to extend the primary knowledge in some greedy way. 3. Algorithms and strategies The operations listed in Section 2.3 give a great deal of freedom in the discovery of solutions. However, not all strategies for applying these operations are equally efficient. We distinguish between static strategies, where the graph is fixed throughout the discovery process (only basic operation 1 of Section 2.3 is used) and dynamic strategies, where vertices and edges may be added (operations 2 and 3). 3.1 A naive dynamic strategy To visualize this strategy in our framework jump ahead and to the flower graph in Fig. 3. Start with the seed solution at the vertex $$v_{0}$$ and proceed creating loops as petals in this graph: e.g., use basic operations 2 and 3 to create $$v_{1}$$ and two edges between $$v_{0}$$ and $$v_{1}$$, track the known solutions at $$v_{0}$$ along the new petal to potentially find new solutions at $$v_{0}$$, then ‘forget’ the petal and create an entirely new one in the next iteration. This strategy populates the fibre $$\pi ^{-1}\left (F_{p_{0}}\right )$$, but how fast? Assume the permutation induced by a petal permutation on $$\pi ^{-1}\left (F_{p_{1}}\right )$$ is uniformly distributed. Then for the first petal the probability of finding a new solution equals (d − 1)/d where $$d=\left |\pi ^{-1}\left (F_{p_{1}}\right )\right |$$, which is large. However, for the other petals the probability of arriving at anything new at the end of one tracked path decreases as the known solution set grows. Finding the expected number of iterations (petals) to discover the entire fibre is equivalent to solving the coupon collector’s problem. The number of iterations is d ℓ(d) where $$\ell (d):={1\over 1}+{1\over 2}+\cdots +{1\over d}$$. The values of ℓ(d) can be regarded as lower and upper sums for two integrals of the function $$x\mapsto x^{-1}$$, leading to the bounds $$\ln (d+1)\leq \ell (d)\leq \ln (d)+1$$. Simultaneously tracking all known points along a petal gives a better complexity, since different paths cannot lead to the same solution. We remark that the existing implementations of numerical irreducible decomposition in Bertini (Bates et al., 2013), PHCpack (Verschelde, 1999) and NumericalAlgebraicGeometry for Macaulay2 (Leykin, 2011) that use monodromy are driven by a version of the naive dynamic strategy. 3.2 Static graph strategies It turns out to be an advantage to reuse the edges of the graph. In a static strategy the graph is fixed and we discover solutions according to the following algorithm. The algorithm can be specialized in several ways. We may choose the graph G, specify a stopping criterion stop, choose a strategy for picking the edge e = (j, k). We address the first choice in Section 3.2.1 by listing several graph layouts that can be used. Stopping criteria are discussed in Section 3.2.2 and Section 3.2.3, while strategies for selecting an edge are discussed in Section 3.2.4. Remark 3.2 We notice that if the stopping criterion is never satisfied the number of paths being tracked by Algorithm 3.1 is at most d|E(G)|, where d is the number of solutions of a generic system. 3.2.1 Two static graph layouts We present two graph layouts to be used for the static strategy (Fig. 3). flower(s,t) The graph consists of a central node$$v_{0}$$ and s additional vertices (number of petals), each connected to $$v_{0}$$ by t edges. completeGraph(s,t) The graph has s vertices. Every pair of vertices is connected by t edges. 3.2.2 Stopping criterion if a solution count is known Suppose the cardinality of the fibre $$\pi ^{-1}\left (F_{p}\right )$$ for a generic value of p is known. Then a natural stopping criterion for our algorithm is to terminate when the set of known solutions $$Q_{i}$$ at any node i reaches that cardinality. In particular, for a generic sparse system with fixed monomial support, we can rely on this stopping criterion due to the BKK bound (Bernstein, 1975) that can be obtained by a mixed volume computation. 3.2.3 Stopping criterion if no solution count is known For a static strategy one natural stopping criterion is saturation of the known solution correspondences along all edges. In this case the algorithm simply can’t derive any additional information. It also makes sense to consider a heuristic stopping criterion based on stabilization. The algorithm terminates when no new points are discovered in a fixed number of iterations. This avoids saturating correspondences unnecessarily. In particular, this could be useful if a static strategy algorithm is a part of the dynamic strategy of Section 3.3. Remark 3.3 In certain cases it is possible to provide a stopping criterion using the trace test (Sommese et al., 2002; Leykin et al., 2016). This is particularly useful when there is an equation in the family $$F_{p}(x)=0$$ that describes a generic hypersurface in the parameter space, e.g., an affine linear equation with indeterminate coefficients. In full generality, one could restrict the parameter space to a generic line and, hence, restrict the solution variety to a curve. Now, thinking of $$F_{p}(x)=0$$ as a system of equations bihomogeneous in p and x, one can use the multihomogeneous trace test (Hauenstein & Rodriguez, 2015; Leykin et al., 2016). We note that the multihomogeneous trace test complexity depends on the degree of the solution variety, which may be significantly higher than the degree d of the covering map, where the latter is the measure of complexity for the main problem in our paper. For instance, the system (7) corresponding to the reaction network in Fig. 4 has four solutions, but an additional set of 11 points is necessary to execute the trace test. See example-traceCRN.m2 at (Duff et al.). 3.2.4 Edge-selection strategy We propose two methods for selecting the edge e in Algorithm 3.1. The default is to select an edge and direction at random. A more sophisticated method is to select an edge and a direction based on the potential of that selection to deliver new information: see the discussion in Example 2.6. Let $$e = v_{i} \overset{(\gamma _{1},\gamma _{2})}{\overleftrightarrow{\hphantom{(\gamma _{1},\gamma _{2})}}}v_{j}$$ be an edge considered in the direction from $$v_{i}$$ to $$v_{j}$$. potentialLowerBound equals the minimal number of new points guaranteed to be discovered by following a chosen homotopy using the maximal batch of starting points $$S_{i}$$. That is, it equals the difference between the numbers of known unmatched points $$\left (|Q_{i}|-|C_{e}|\right )-\left (\left |Q_{j}\right |-|C_{e}|\right ) = \left |Q_{i}\right |-\left |Q_{j}\right |$$ if this difference is positive and 0 otherwise. potentialE equals the expected number of new points obtained by tracking one unmatched point along e. This is the ratio $$d-|Q_{j}|\over d-|C_{e}|$$ of undiscovered points among all unmatched points if $$\left |Q_{i}\right |-|C_{e}|>0$$ and 0 otherwise. Note that potentialE assumes we know the cardinality of the fibre, while potentialLowerBound does not depend on that piece of information. There is a lot of freedom in choosing potentials in our algorithmic framework. The two above potentials are natural ‘greedy’ choices that are easy to describe and implement. It is evident from our experiments (Section 6.1.1) that they may order edges differently resulting in varying performance. 3.3 An incremental dynamic graph strategy Consider a dynamic strategy that amounts to augmenting the graph once one of the above ‘static’ criteria terminates Algorithm 3.1 for the current graph. One simple way to design a dynamic stopping criterion, we call it dynamic stabilization, is to decide how augmentation is done and fix the number of augmentation steps that the algorithm is allowed to make without increasing the solution count. A dynamic strategy, which is simple to implement, is one that starts with a small graph G and augments it if necessary. We emphasize that the criteria described in this subsection and parts of Section 3.2.3 are heuristic and there is a lot of freedom in designing such. In Section 6.2 we successfully experiment using a static stabilization criterion with some examples, for which the solution count is generally not known. 4. Statistical analysis The directed cycles in the graph G starting and ending at a vertex $$v_{1}$$ give elements of the fundamental group $$\pi _{1}(B\setminus D)$$, which correspond to the elements of the monodromy subgroup M(G) of the monodromy group $$M(\pi _{1}(B\setminus D))$$. The latter is a subgroup of $${\mathrm{S}}_{d}$$, where $$d = \left |\pi ^{-1}\left (F_{p_{1}}\right )\right |$$. For example, if $$G={\tt completeGraph}(2,j+1)$$ then the j cycles produced by edges $$e_{1}$$ and $$e_{2}$$, $$e_{2}$$ and $$e_{3},\dots,e_{j}$$ and $$e_{j+1}$$ suffice to generate M(G). The minimal number j of cycles necessary to generate M(G) in the general case is $$\beta _{1}(G)$$, the first Betti number of G as a topological space. For the purpose of simplifying statistical analysis, we assume that picking a random decorated graph G with $$j=\beta _{1}(G)$$ induces uniformly and independently distributed permutations $$\sigma _{1}, \ldots , \sigma _{j} \in {\mathrm{S}}_{d}$$, where $${\mathrm{S}}_{d}$$ is the symmetric group acting on the fibre $$\pi ^{-1}\left (F_{p_{1}}\right )$$. It would be hard in practice to achieve uniformity even when the monodromy group is a full symmetric group: see Section 5.1. 4.1 The probability of a transitive action Suppose the number of solutions d is known and stop(d) denotes the corresponding stopping criterion. Our aim is to analyze the probability of producing the full solution set via Algorithm 3.1 or, equivalently, the probability of $${\tt dynamicMonodromySolve}(G,x_{1},{\tt stop(d)},{\tt augment})$$ terminating after at most j iterations, assuming that $$\beta _{1}(G)=j$$ at the jth iteration. This equals the probability of $$\langle \sigma _{1}, \ldots , \sigma _{j}\rangle$$ acting transitively, i.e., $$\Pr \left [X_{d}\leq j\right ]$$, where $$X_{d}$$ is the random variable $$X_{d} = \displaystyle\inf \left\{ i \in \mathbb{N} \mid \langle \sigma_{1}, \ldots, \sigma_{i} \rangle \textrm{ is transitive} \right\}\!.$$ When d > 1 we have $$\Pr [X_{d} =0]=0$$, while $$\Pr [X_{d} =1 ]$$ is proportional to the number of d-cycles in the monodromy group. When the monodromy group is full symmetric, we can compute and give asymptotic estimates for the distribution of $$X_{d}.$$ The following theorem is a generalization of a result by Dixon, regarding the case j = 2. The proof we give in Section 4.2 follows the strategy of the study by Dixon (1969). Theorem 4.1 For j ≥ 2, $$\Pr \left [ X_{d} \le j \right ] = 1 - d^{1-j} + R_{j} (d),$$ where the error term $$R_{j}$$ satisfies $$\left |R_{j}(d)\right |=O\left (d^{-j}\right )$$. Remark 4.2 As a corollary one can deduce that the expected value of $$X_{d}$$ is asymptotically finite and $$\mathrm{E}\left [X_{d}\right ]\to 2$$ as $$d\to \infty$$. The numerical approximations in Table 1 show that $$\mathrm{E}\left [X_{d}\right ]\leq 2.1033$$ for all d. Moreover, the proof in Section 4.2 implies that $$\left |R_{j}(d)\right |<C\left ({d\over 2}\right )^{-j}$$ with the constant C not depending on j. Therefore, $$\Pr \left [ X_{d}> j \right ]$$ decays exponentially with j. Under the idealistic assumption that new cycles in the graph lead to independently and uniformly distributed permutations of the fibre $$Q_{1}$$ the expected Betti number needed for completion in Algorithm 3.4 is at most 2.1033. If we assume that augment increases the Betti number by one by adding at most a fixed number of edges then the expected number of tracked paths is linear in d. Remark 4.3 We point out that Babai (1989) proved Dixon’s conjecture stating that the subgroup of $${\mathrm{S}}_{d}$$ generated by two random permutations is $${\mathrm{S}}_{d}$$ or $${\mathrm{A}}_{d}$$ with probability $$1 - d^{-1} + O\left (d^{-2}\right )$$. This shows that other subgroups are rare. However, it is easy to construct families with a transitive monodromy group that is neither full symmetric nor alternating. For example, take $${x_{1}^{2}}-c_{1}={x_{2}^{2}}-c_{2}=0$$ with irreducible solution variety, and four solutions for generic choices of $$c_{1}$$ and $$c_{2}$$. Tracking two solutions with the same $$x_{1}$$ coordinate, as $$c_{1}$$ and $$c_{2}$$ vary, the moving points on the tracked paths will continue to have equal projections to $$x_{1}$$. The monodromy group is $${\mathbb{Z}}_{2}\times {\mathbb{Z}}_{2}$$. We reiterate that generators for the monodromy group are seldomly known a priori. Computing them is likely to be prohibitively expensive, and the probability distribution with which our algorithm picks elements of the monodromy group is unknown, as it is prohibitively hard to analyze. 4.2 Proof of a generalization of Dixon’s theorem Fix any integer j ≥ 2. We wish to prove Theorem 4.1 by estimating the quantity $$t_{d} = \Pr \left[ \langle \sigma_{1}, \sigma_{2}, \ldots, \sigma_{j} \rangle \textrm{ is transitive}\right]\!,$$ where $$\sigma _{1}, \ldots , \sigma _{j}$$ are independent and uniformly distributed on $${\mathrm{S}}_{d}$$. Suppose we partition the set {1, 2, … , d} in such a way that there are $$k_{i}$$ classes of size i for each 1 ≤ i ≤ d. All such partitioning schemes are indexed by the set $$K_{d} = \left\{ \vec{k} \in{\mathbb{N}}^{d} \left|\ \sum \right. i \, k_{i} = d \right\}\!.$$ The number of partitions corresponding to each $$\vec{k} \in K_{d}$$ is $${d!}/\big (\!\prod _{i=1}^{d} (i!\!)^{k_{i}} \cdot k_{i}!\!\big )$$. For each $$\vec{k} \in K_{d},$$ the partition given by the orbits of $$\langle \sigma _{1}, \ldots , \sigma _{j} \rangle$$ is $$\vec{k}$$-indexed precisely when this group acts transitively on all classes of some partition associated to $$\vec{k}.$$ The number of tuples in $${S_{i}^{j}}$$ with coordinates generating a group acting transitively on $$\{1,\dots ,i\}$$ is $$t_{i} \,\left (i!\!\right )^{j}$$. Thus, we may count the set $$\underbrace{ {\mathrm{S}}_{d} \times \cdots \times {\mathrm{S}}_{d}}_{j}$$ as \begin{align*} \left(d!\!\right)^{j} &= \displaystyle\sum_{\vec{k} \in K_{d}} \displaystyle\frac{d!}{\prod_{i=1}^{d} (i!\!)^{k_{i}} \cdot k_{i}!} \cdot \displaystyle\prod_{i=1}^{d} \big( t_{i} \left(i!\!\right)^{j} \big)^{k_{i}} \\[-2pt] &= d! \cdot \displaystyle\sum_{\vec{k} \in K_{d}}\displaystyle\prod_{i=1}^{d} \displaystyle\frac{\left( t_{i} \, (i!\!)^{j-1} \right)^{k_{i}}}{k_{i}!}. \end{align*} Let $$\widehat{F}$$ denote the generating function of the sequence $$F(d) = (d!\!)^{j-1}.$$ Note the formal identity $$\exp \left( \displaystyle\sum_{i=0}^{\infty } y_{i} x^{i} \right) = \displaystyle\sum_{d=0}^{\infty} x^{d} \, \sum_{\vec{k} \in K_{d} } \prod_{i=1}^{d} \frac{y_{i}^{k_{i}}}{k_{i}!},$$ which follows by letting g(x) denote the right-hand side as a formal power series in x, $$f(x) = \sum y_{i} x^{i},$$ and noting the equivalent form f′g = g′ with $$f(0)=y_{0}.$$ We have \begin{align*} \displaystyle\sum_{d=1}^{\infty} d\cdot (d!\!)^{j-1} \, x^{d-1} &= \frac{d}{d x} \, \widehat{F} (x) \\ &= \frac{d}{dx} \, \exp \left( \displaystyle\sum_{i=1}^{\infty } t_{i} \, (i!\!)^{j-1} x^{i} \right)\\ &= \left( \displaystyle\sum_{d=0}^{\infty} (d!\!)^{j-1} \, x^{d} \right) \cdot \displaystyle\sum_{i=1}^{\infty } i \cdot t_{i} \, (i!\!)^{j-1} x^{i-1} \\ &= \displaystyle\sum_{d^{\prime}=1}^{\infty } x^{d^{\prime}-1} \, \left( \displaystyle\sum_{i=1}^{d^{\prime}} i \cdot t_{i} \, \big( i! \cdot \left(d^{\prime}-i\right)! \big)^{j-1} \right)\!, \end{align*} where the first equation follows by formal differentiation of the power series $$\widehat{F}$$, the second from the two identities above with properly substituted values for $$y_{i}$$, the third by applying the chain rule and the definition of $$\widehat{F}$$ and the fourth by rearranging terms by index substitution $$d^{\prime}=i+d$$. Upon equating coefficients of $$x^{d-1}$$ for $$d=1,\dots$$ we obtain $$d = \displaystyle\sum_{i=1}^{d} \binom{d}{i}^{1-j} \, i\, t_{i}.$$ (4.1) Table 1 Numerical approximations of $$t_{d}$$—the probability of the j random permutations acting transitively on a fibre of size d for j = 2, 3, 4. After computing these values for larger j, a numerical approximation of $$\mathrm{E}\left [X_{d}\right ]$$ is extracted d j = 2 j = 3 j = 4 $$\mathrm{E}\left [X_{d}\right ]$$ 1 1 1 1 0 2 0.75 0.875 0.9375 2 3 0.72222222 0.89814815 0.96450617 2.10000000 4 0.73958333 0.93012153 0.98262080 2.10329381 5 0.76833333 0.95334722 0.99115752 2.08926525 10 0.88180398 0.98954768 0.99898972 2.02976996 20 0.94674288 0.99747856 0.99987487 2.00591026 30 0.96536852 0.99888488 0.99996295 2.00245160 d j = 2 j = 3 j = 4 $$\mathrm{E}\left [X_{d}\right ]$$ 1 1 1 1 0 2 0.75 0.875 0.9375 2 3 0.72222222 0.89814815 0.96450617 2.10000000 4 0.73958333 0.93012153 0.98262080 2.10329381 5 0.76833333 0.95334722 0.99115752 2.08926525 10 0.88180398 0.98954768 0.99898972 2.02976996 20 0.94674288 0.99747856 0.99987487 2.00591026 30 0.96536852 0.99888488 0.99996295 2.00245160 Table 1 Numerical approximations of $$t_{d}$$—the probability of the j random permutations acting transitively on a fibre of size d for j = 2, 3, 4. After computing these values for larger j, a numerical approximation of $$\mathrm{E}\left [X_{d}\right ]$$ is extracted d j = 2 j = 3 j = 4 $$\mathrm{E}\left [X_{d}\right ]$$ 1 1 1 1 0 2 0.75 0.875 0.9375 2 3 0.72222222 0.89814815 0.96450617 2.10000000 4 0.73958333 0.93012153 0.98262080 2.10329381 5 0.76833333 0.95334722 0.99115752 2.08926525 10 0.88180398 0.98954768 0.99898972 2.02976996 20 0.94674288 0.99747856 0.99987487 2.00591026 30 0.96536852 0.99888488 0.99996295 2.00245160 d j = 2 j = 3 j = 4 $$\mathrm{E}\left [X_{d}\right ]$$ 1 1 1 1 0 2 0.75 0.875 0.9375 2 3 0.72222222 0.89814815 0.96450617 2.10000000 4 0.73958333 0.93012153 0.98262080 2.10329381 5 0.76833333 0.95334722 0.99115752 2.08926525 10 0.88180398 0.98954768 0.99898972 2.02976996 20 0.94674288 0.99747856 0.99987487 2.00591026 30 0.96536852 0.99888488 0.99996295 2.00245160 Remark 4.4 Equation (1) gives a list of linear equations in the probabilities $$t_{1},t_{2},\dots$$ allowing us to successively determine these values by backward substitution. In Table 1 we list some solutions for j = 2, 3, 4. To complete the proof of Theorem 4.1, we introduce, as in Dixon’s proof, the auxiliary quantities \begin{equation*} r_{d} = d\, \left(1 - t_{d}\right)\ \ \textrm{and}\ \ \ c_{d} = \displaystyle\sum_{i=1}^{d-1} \binom{d}{i}^{1-j} \, i . \end{equation*} Noting that $$\binom{d}{i}^{1-j}i+\binom{d}{d-i}^{1-j}(d-i)={d\over 2}\left (\binom{d}{i}^{1-j}+\binom{d}{d-i}^{1-j}\right )$$, we have \begin{align} {c_{d}\over d} =\frac{1}{2} \cdot \displaystyle\sum_{i=1}^{d-1} \binom{d}{i}^{1-j} &= d^{1-j} + \left[ \binom{d}{2}^{1-j} + \frac{1}{2} \displaystyle\sum_{i=3}^{d-3} \binom{d}{i}^{1-j} \right] \nonumber\\ &\leq d^{1-j} + \left[ \binom{d}{2}^{1-j} + \frac{1}{2} (d-5) \binom{d}{3}^{1-j} \right]\!. \end{align} (4.2) From j ≥ 2 it follows that the bracketed expression in (2) is $$O\left (d^{-j}\right ).$$ Using (1), $$t_{i}=1-{r_{i}\over i}\leq 1-{r_{i}\over d}$$ and the definition of $$c_{d}$$ we may bound $$r_{d}$$: \begin{align*} r_{d} &= d(1-t_{d})=-dt_{d}+d=-dt_{d}+\displaystyle\sum_{i=1}^{d} \binom{d}{i}^{1-j} \, i\cdot t_{i} \\ &=\displaystyle\sum_{i=1}^{d-1} \binom{d}{i}^{1-j} \, i \cdot t_{i} \leq\displaystyle\sum_{i=1}^{d-1} \binom{d}{i}^{1-j} \, i \cdot \left(1-{r_{i}\over d}\right)\nonumber \\ &= \sum_{i=1}^{d-1} \binom{d}{i}^{1-j}i -{1\over d}\displaystyle\sum_{i=1}^{d-1} r_{i} \, \binom{d}{i}^{1-j}= c_{d} -{1\over d}\displaystyle\sum_{i=1}^{d-1} r_{i} \, \binom{d}{i}^{1-j}\leq c_{d}. \end{align*} (3) (4) To bound the error term $$R_{j}(d):=t_{d} - \left (1-d^{1-j}\right )$$, we consider first the case where its sign is positive. Expanding $$t_{d}=1-{r_{d}\over d}$$ using (3) above and $$i\, t_{i}=i-{r_{i}}$$, \begin{align*} t_{d} - \big(1-d^{1-j}\big) &= 1-\sum_{i=1}^{d-1} \binom{d}{i}^{1-j}{it_{i}\over d}-1+d^{1-j}=d^{1-j}-\sum_{i=1}^{d-1} \binom{d}{i}^{1-j}{it_{i}\over d}\\ &=d^{1-j}-\sum_{i=1}^{d-1} \binom{d}{i}^{1-j}{i\over d}+\sum_{i=1}^{d-1} \binom{d}{i}^{1-j}{r_{i}\over d}\\ &=\underbrace{d^{1-j}-d^{1-j}\left({1\over d}+{d-1\over d}\right)}_{=0} \, - \, \sum_{i=2}^{d-2} \binom{d}{i}^{1-j}{i\over d} \, + \, \sum_{i=1}^{d-1} \binom{d}{i}^{1-j}{r_{i}\over d}, \end{align*} and we may focus on the last summation to get, for d ≥ 2, \begin{align*} t_{d} - \big(1-d^{1-j}\big) &\leq d^{1-j}{r_{d-1}\over d-1}+\sum_{i=2}^{d-2} \binom{d}{i}^{1-j}{r_{i}\over d}\qquad\qquad(\text{note }r_{1} = 0)\\ &\leq d^{1-j}{c_{d-1}\over d-1}+\sum_{i=2}^{d-2} \binom{d}{2}^{1-j}{c_{i}\over d}\qquad\qquad\qquad(\text{by }(4))\\ &\leq d^{1-j}{c_{d-1}\over d-1}+\left({d^{2}\over 4}\right)^{1-j}\, \sum_{i=2}^{d-2}{c_{i}\over i}\\ &=O\big(d^{2-2j}\big)+O\big(d^{(2-2j)+(2-j)}\big) \\ &= O\big(d^{-j}\big).\qquad\qquad\qquad\qquad\qquad(\text{since }j\geq 2) \end{align*} The case where $$R_{j}(d)\le 0$$ may be handled similarly, using $$t_{d}=1-{r_{d}\over d}$$, $$r_{d}\leq c_{d}$$ and what we know about the content of the bracket in (2). \begin{align*} -t_{d} + \big(1-d^{1-j}\big)&= -\left(1-{r_{d} \over d}\right)+\big(1-d^{1-j}\big)={r_{d} \over d}-d^{1-j}\leq{c_{d} \over d}-d^{1-j}\\ &\leq \left[ \binom{d}{2}^{1-j} + \frac{1}{2} (d-5) \binom{d}{3}^{1-j} \right] = O\big(d^{-j}\big). \end{align*} 5. Implementation We implement the package MonodromySolver in Macaulay2 (Grayson & Stillman) using the functionality of the package NumericalAlgebraicGeometry (Leykin, 2011). The source code and examples used in the experiments in the next section are available at Duff et al. The main function monodromySolve realizes Algorithms 3.1 and 3.4, see the documentation for details and many options. The tracking of homotopy paths in our experiments is performed with the native routines implemented in the kernel of Macaulay2, however, NumericalAlgebraicGeometry provides an ability to outsource this core task to an alternative tracker (PHCpack or Bertini). Main auxiliary functions—createSeedPair, sparseSystemFamily, sparseMonodromySolve and solveSystemFamily—are there to streamline the user’s experience. The last two are blackbox routines that don’t assume any knowledge of the framework described in this paper. The overhead of managing the data structures is supposed to be negligible compared to the cost of tracking paths. However, since our implementation uses the interpreted language of Macaulay2 for other tasks, this overhead could be sizable (up to 10% for large examples in Section 6). Nevertheless, most of our experiments are focused on measuring the number of tracked paths as a proxy for computational complexity. Remark 5.1 This paper’s discussion focuses on linear parametric systems with a nonempty dominant component. However, the implementation works for other cases where our framework can be applied. For instance, if the system is linear in parameters, but has no dominant component, there may still be a unique ‘component of interest’ with a straightforward way to produce a seed pair. This is so, for instance, in the problem of finding the degree of the variety $${\mathop{SO}}(n)$$, which we use in Table 7. The point x is restricted to $${\mathop{SO}}(n)$$, the special orthogonal group, which is irreducible as a variety. This results in a unique ‘component of interest’ in the solution variety, the one that projects onto $${\mathop{SO}}(n)$$, see the study by Brandt et al. (2017) for details. In a yet more general case of a system that is nonlinear in parameters, it is still possible to use our software. We outline the theoretical issues one would need to consider in 4 of Section 7. 5.1 Randomization Throughout the paper we refer to random choices we make that we assume avoid various nongeneric loci. For implementation purposes we make simple choices. For instance, the vertices of the graph get distributed uniformly in a cube in the base space with the exception of the seeded vertex: createSeedPair picks $$(p_{0},x_{0})\in B\times {\mathbb{C}}^{n}$$ by choosing x uniformly in a cube, then choosing $$p_{0}$$ uniformly in a box in the subspace $$\left \{ p \mid F_{p}(x)=0 \right \}$$. A choice of probability distribution on B translates to some (discrete) distribution on the symmetric group $${\mathrm{S}}_{d}$$. However, it is simply too hard to analyze—there are virtually no studies in this direction. We make the simplest possible assumption of uniform distribution on $${\mathrm{S}}_{d}$$ in order to perform the theoretical analysis in Section 4, and shed some light on why our framework works well. There is an interesting, more involved, alternative to this assumption in the studies by Galligo & Poteaux (2011) and Galligo & Miclo (2012), which relies on the intuition in the case n = 1. 5.2 Solution count The BKK bound, computed via mixed volume, is used as a solution count in the examples of sparse systems in Sections 6.1.1 and 6.1.2. In the latter we compute mixed volume via a closed formula that involves permanents, while the former relies on general algorithms implemented in several software packages. Our current implementation uses PHCpack (Verschelde, 1999), which incorporates the routines of MixedVol further developed in Hom4PS-2 (Lee et al., 2008). Other alternatives are pss5 (Malajovich, 2017) and Gfanlib (Jensen, 2016). While some implementations are randomized the latter uses symbolic perturbations to achieve exactness. The computation of the mixed volume is not a bottleneck in our algorithm. The time spent in that preprocessing stage is negligible compared to the rest of the computation. 5.3 Certification The reader should realize that the numerical homotopy continuation we use is driven partly by heuristics. As a post-processing step we can certify (i.e., formally prove) the completeness and correctness of the solution set to a polynomial system computed with our main method. This is possible in the scenario when the parametric system is square, all solutions are regular (the Jacobian of the system is invertible), and the solution count is known. We can use Smale’s $$\alpha$$-theory (Blum et al., 1998, Section 8) to certify an approximation to a regular solution of a square system. In a Macaulay2 package NumericalCertification, we implement a numerical version of an $$\alpha$$-test after finding an approximate solution to certify that our solution is an approximate zero in a rigorous sense. One of the main functions of NumericalCertification is certifySolutions, which determines whether the given solution is an approximate zero of the given polynomial system. It also produces an upper bound on the distance from the approximation to the exact solution to which it is associated. See paper-examples/example-NashCertify.m2 at (Duff et al.), which is an example of an $$\alpha$$-test application to the solutions of a problem described in Section 6.1.2. In the implementation of certification all arithmetic and linear algebra operations are done over the field of Gaussian rationals, $${\mathbb{Q}}[ {\mathbf{i}}]/\big ( {\mathbf{i}}^{2}+1\big )$$. To use this certification method we first convert the coefficients of the system to Gaussian rationals, then perform certification numerically. See the study by Hauenstein & Sottile (2012) for a stand-alone software package alphaCertified and detailed implementation notes. 6. Experiments In this section we first report on experiments with our implementation and various examples in Sections 6.1 and 6.2. We then investigate the completion rate of Algorithm 3.1 in Section 6.3. Finally, we compare against other software in Section 6.4. Table 2 Cyclic 7 experimental results for the flower strategy (#vertices-1, edge multiplicity) (3,2) (4,2) (5,2) (3,3) (4,3) |E(G)| 6 8 10 9 12 $$\beta _{1}(G)$$ 3 4 5 6 8 |E(G)| ⋅ 924 5544 7392 9240 8316 11088 completion rate 100% 100% 100% 100% 100% Random Edge 5119 6341 7544 6100 7067 potentialLowerBound 5252 6738 8086 6242 7886 potentialE 4551 5626 6355 4698 5674 (#vertices-1, edge multiplicity) (3,2) (4,2) (5,2) (3,3) (4,3) |E(G)| 6 8 10 9 12 $$\beta _{1}(G)$$ 3 4 5 6 8 |E(G)| ⋅ 924 5544 7392 9240 8316 11088 completion rate 100% 100% 100% 100% 100% Random Edge 5119 6341 7544 6100 7067 potentialLowerBound 5252 6738 8086 6242 7886 potentialE 4551 5626 6355 4698 5674 View Large Table 2 Cyclic 7 experimental results for the flower strategy (#vertices-1, edge multiplicity) (3,2) (4,2) (5,2) (3,3) (4,3) |E(G)| 6 8 10 9 12 $$\beta _{1}(G)$$ 3 4 5 6 8 |E(G)| ⋅ 924 5544 7392 9240 8316 11088 completion rate 100% 100% 100% 100% 100% Random Edge 5119 6341 7544 6100 7067 potentialLowerBound 5252 6738 8086 6242 7886 potentialE 4551 5626 6355 4698 5674 (#vertices-1, edge multiplicity) (3,2) (4,2) (5,2) (3,3) (4,3) |E(G)| 6 8 10 9 12 $$\beta _{1}(G)$$ 3 4 5 6 8 |E(G)| ⋅ 924 5544 7392 9240 8316 11088 completion rate 100% 100% 100% 100% 100% Random Edge 5119 6341 7544 6100 7067 potentialLowerBound 5252 6738 8086 6242 7886 potentialE 4551 5626 6355 4698 5674 View Large Table 3 Cyclic 7 experimental results for the completeGraph strategy (#vertices, edge multiplicity) (2,3) (2,4) (2,5) (3,2) (4,1) |E(G)| 3 4 5 6 6 $$\beta _{1}(G)$$ 2 3 4 4 3 |E(G)| ⋅ 924 2772 3698 4620 5544 5544 completion rate 65% 80% 90% 100% 100% Random Edge 2728 3296 3947 4805 5165 potentialLowerBound 2727 3394 3821 4688 5140 potentialE 2692 2964 2957 3886 4380 (#vertices, edge multiplicity) (2,3) (2,4) (2,5) (3,2) (4,1) |E(G)| 3 4 5 6 6 $$\beta _{1}(G)$$ 2 3 4 4 3 |E(G)| ⋅ 924 2772 3698 4620 5544 5544 completion rate 65% 80% 90% 100% 100% Random Edge 2728 3296 3947 4805 5165 potentialLowerBound 2727 3394 3821 4688 5140 potentialE 2692 2964 2957 3886 4380 Table 3 Cyclic 7 experimental results for the completeGraph strategy (#vertices, edge multiplicity) (2,3) (2,4) (2,5) (3,2) (4,1) |E(G)| 3 4 5 6 6 $$\beta _{1}(G)$$ 2 3 4 4 3 |E(G)| ⋅ 924 2772 3698 4620 5544 5544 completion rate 65% 80% 90% 100% 100% Random Edge 2728 3296 3947 4805 5165 potentialLowerBound 2727 3394 3821 4688 5140 potentialE 2692 2964 2957 3886 4380 (#vertices, edge multiplicity) (2,3) (2,4) (2,5) (3,2) (4,1) |E(G)| 3 4 5 6 6 $$\beta _{1}(G)$$ 2 3 4 4 3 |E(G)| ⋅ 924 2772 3698 4620 5544 5544 completion rate 65% 80% 90% 100% 100% Random Edge 2728 3296 3947 4805 5165 potentialLowerBound 2727 3394 3821 4688 5140 potentialE 2692 2964 2957 3886 4380 6.1 Sparse polynomial systems The example families in this subsection have the property that the support of the equations is fixed, while the coefficients can vary freely, as long as they are generic. We run the static graph strategy Algorithm 3.1 on these examples. Our timings do not include the $$\alpha$$-test, which was only applied in Section 6.1.2. 6.1.1 Cyclic roots The cyclic n-roots polynomial system is $$\begin{cases} i = 1, 2, 3, 4, \ldots, n-1: \displaystyle\sum_{j=0}^{n-1} ~ \prod_{k=j}^{j+i-1} x_{k~\textrm{mod}~n}=0 \\ x_{0}x_{1}x_{2} \cdots x_{n-1} - 1 = 0. \end{cases}$$ (6.1) This system is commonly used to benchmark polynomial system solvers. We will study the modified system with randomized coefficients and seek solutions in $$( {\mathbb{C}}\setminus \{0\})^{n}$$. Therefore, the solution count can be computed as the mixed volume of the Newton polytopes of the left-hand sides, providing a natural stopping criterion discussed in Section 3.2.2. This bound is 924 for cyclic 7. Tables 2 and 3 contain averages of experimental data from running 20 trials of Algorithm 3.1 on cyclic 7. The main measurement reported is the average number of paths tracked, as the unit of work for our algorithm is tracking a single homotopy path. The experiments were performed with 10 different graph layouts and three edge-selection strategies. With respect to number of paths tracked, we see that it is an advantage to keep the Betti number high and edge number low. Remark 6.1 Computing the expected success rates (> 99%) using Remark 4.4, we conclude that the resulting permutations do not conform to the model of picking uniformly from $$S_{924}$$. The completion rate depends on the choice of strategy (compare Table 2 to Table 3). Nevertheless, both in theory (assuming uniform distribution as in Section 4) and in practice (with distribution unknown to us), the completion rate does converge to 100% rapidly as the Betti number grows. 6.1.2 Nash equilibria Semimixed multihomogeneous systems arise when one is looking for all totally mixed Nash equilibria (TMNE) in game theory. A specialization of mixed volume using matrix permanents gives a concise formula for a root count for systems arising from TMNE problems (Emiris & Vidunas, 2014). We provide an overview of how such systems are constructed based on the study by Emiris & Vidunas (2014). Suppose there are N players with m options each. For player i ∈ {1, … , N} using option j ∈ {1, … , m} we have the equation $$P_{j}^{(i)} = 0$$, where \begin{align} \begin{aligned} P_{j}^{(i)} = \sum_{\substack{k_{1},\ldots,k_{i-1},\\ k_{i+1},\ldots, k_{N}}} a^{(i)}_{k_{1},\ldots, k_{i-1},\, j, k_{i+1}, \ldots, k_{N}}p_{k_{1}}^{(1)}p_{k_{2}}^{(2)}\cdots p_{k_{i-1}}^{(i-1)}p_{k_{i+1}}^{(i+1)} \cdots p_{k_{N}}^{(N)}. \end{aligned} \end{align} (6.2) The parameters $$a^{(i)}_{k_{1},k_{2},\ldots ,k_{N}}$$ are the pay-off rates for player i when players $$1,\ldots ,i-1,i+1,\ldots , N$$ are using options $$k_{1}, \ldots ,k_{i-1},k_{i+1},\ldots ,k_{N}$$, respectively. Here the unknowns are $$p^{(i)}_{k_{j}}$$, representing the probability that player i will use option $$k_{j} \in \{1,\ldots , m\}$$. There is one constraint on the probabilities for each player i ∈ {1, …, N}, namely the condition that \begin{align} p_{1}^{(i)}+p_{2}^{(i)}+\cdots+p_{m}^{(i)}=1. \end{align} (6.3) The system (4) consists of N ⋅ m equations in N ⋅ m unknowns. Using condition (5) reduces the number of unknowns to N(m − 1). Finally, we eliminate the $$P_{j}^{(i)}$$ by constructing \begin{align} P_{1}^{(i)}=P_{2}^{(i)}, ~ P_{1}^{(i)}=P_{3}^{(i)}, ~ \ldots ~, P_{1}^{(i)}=P_{m}^{(i)},\quad \textrm{ for each } i \in \{1,\ldots, N\}. \end{align} (6.4) The final system is a square system of N(m − 1) equations in N(m − 1) unknowns. For one of our examples (paper-examples/example-Nash.m2 at (Duff et al.)), we chose the generic system of this form for N = 3 players with m = 3 options for each. The result is a system of six equations in six unknowns and 81 parameters with 10 solutions. We also use this example to demonstrate that these solutions can be certified using NumericalCertification (Section 5.3). 6.2 Chemical reaction networks A family of interesting examples arises from chemical reaction network theory. A chemical reaction network considered under the laws of mass action kinetics leads to a dynamical polynomial system, the solutions of which represent all the equilibria for the given reaction network (MacLean et al., 2015; Gross et al., 2016). These polynomial systems are not generically sparse and we cannot easily compute their root count. In our experiments we used the stabilization stopping criterion, terminating the algorithm after a fixed number of iterations that do not deliver new points; the default is 10 fruitless iterations. Figure 4 gives an example of a small chemical reaction network. Fig. 4. View largeDownload slide Chemical reaction network example. Fig. 4. View largeDownload slide Chemical reaction network example. Applying the laws of mass action kinetics to the reaction network above we obtain the polynomial system (7) consisting of the corresponding steady-state and conservation equations. Here the $$k_{i}$$s represent the reaction rates, $$x_{i}$$s represent species concentrations and the $$c_{i}$$s are parameters. \begin{align} \begin{aligned} \dot{x_{A}} & = k_{1}{x_{B}^{2}}-k_{2}x_{A}-k_{3}x_{A}x_{C}+k_{4}x_{D}+k_{5}x_{B}x_{E} \\ \dot{x_{B}} & = 2k_{1}x_{A}-2k_{2}{x_{B}^{2}}+k_{4}x_{D}-k_{5}x_{B}x_{E} \\ \dot{x_{C}} & = -k_{3}x_{A}x_{C}+k_{4}x_{D}+k_{5}x_{B}x_{E} \\ \dot{x_{D}} & = k_{3}x_{A}x_{C}-(k_{4}+k_{6})x_{D} \\ \dot{x_{E}} & = -k_{5}x_{B}x_{E}+k_{6}x_{D} \\ 0 & = 2x_{A}+x_{B}-x_{C}+x_{D}-c_{1} \\ 0 & = -2x_{A}-x_{B}+2x_{C}+x_{E}-c_{2}. \end{aligned} \end{align} (6.5) Typically, systems resulting from chemical reaction networks will be overdetermined. With the current implementation one needs to either square the system or use a homotopy tracker that supports following a homotopy in a space of overdetermined systems. Table 4 Katsura-(n − 1) for the flower strategy (#vertices-1, edge multiplicity) n BKK Bound (3,2) (4,2) (5,2) (3,3) (4,3) $$\beta _{1}=3$$ $$\beta _{1}=4$$ $$\beta _{1}=5$$ $$\beta _{1}=6$$ $$\beta _{1}=8$$ 5 12 96.4% 99.4% 99.6% 100% 99.8% 6 30 98.6% 100% 99.8% 100% 99.6% 7 54 97.6% 98.8% 99.4% 99.4% 98.4% 8 126 99.2% 99.8% 99.6% 99.8% 99.8% 9 240 98.8% 99.6% 98.4% 98.4% 98.6% 10 504 98.6% 98.8% 99.2% 99.4% 98.8% (#vertices-1, edge multiplicity) n BKK Bound (3,2) (4,2) (5,2) (3,3) (4,3) $$\beta _{1}=3$$ $$\beta _{1}=4$$ $$\beta _{1}=5$$ $$\beta _{1}=6$$ $$\beta _{1}=8$$ 5 12 96.4% 99.4% 99.6% 100% 99.8% 6 30 98.6% 100% 99.8% 100% 99.6% 7 54 97.6% 98.8% 99.4% 99.4% 98.4% 8 126 99.2% 99.8% 99.6% 99.8% 99.8% 9 240 98.8% 99.6% 98.4% 98.4% 98.6% 10 504 98.6% 98.8% 99.2% 99.4% 98.8% Table 4 Katsura-(n − 1) for the flower strategy (#vertices-1, edge multiplicity) n BKK Bound (3,2) (4,2) (5,2) (3,3) (4,3) $$\beta _{1}=3$$ $$\beta _{1}=4$$ $$\beta _{1}=5$$ $$\beta _{1}=6$$ $$\beta _{1}=8$$ 5 12 96.4% 99.4% 99.6% 100% 99.8% 6 30 98.6% 100% 99.8% 100% 99.6% 7 54 97.6% 98.8% 99.4% 99.4% 98.4% 8 126 99.2% 99.8% 99.6% 99.8% 99.8% 9 240 98.8% 99.6% 98.4% 98.4% 98.6% 10 504 98.6% 98.8% 99.2% 99.4% 98.8% (#vertices-1, edge multiplicity) n BKK Bound (3,2) (4,2) (5,2) (3,3) (4,3) $$\beta _{1}=3$$ $$\beta _{1}=4$$ $$\beta _{1}=5$$ $$\beta _{1}=6$$ $$\beta _{1}=8$$ 5 12 96.4% 99.4% 99.6% 100% 99.8% 6 30 98.6% 100% 99.8% 100% 99.6% 7 54 97.6% 98.8% 99.4% 99.4% 98.4% 8 126 99.2% 99.8% 99.6% 99.8% 99.8% 9 240 98.8% 99.6% 98.4% 98.4% 98.6% 10 504 98.6% 98.8% 99.2% 99.4% 98.8% Although we may obtain large systems they typically have very low root counts compared to the sparse case. The polynomial system (7) has four solutions. A larger example is the wnt signaling pathway from Systems Biology (Gross et al., 2016) consisting of 19 polynomial equations with nine solutions. All nine solutions are obtained in less than a second with Algorithm 3.1. Table 5 Katsura-(n − 1) for the completeGraph strategy (#vertices, edge multiplicity) n BKK Bound (2,3) (2,4) (2,5) (3,2) (4,1) $$\beta _{1}=2$$ $$\beta _{1}=3$$ $$\beta _{1}=4$$ $$\beta _{1}=4$$ $$\beta _{1}=3$$ 5 12 65.6% 88.2% 95% 99.2% 98% 6 30 77.4% 95.2% 99% 99.8% 99.6% 7 54 74.4% 96.2% 99.2% 99.6% 99.8% 8 126 81.8% 97% 99.2% 100% 99.8% 9 240 85.2% 97.6% 99.4% 99% 98.2% 10 504 89.2% 98.2% 99.2% 99.4% 99% (#vertices, edge multiplicity) n BKK Bound (2,3) (2,4) (2,5) (3,2) (4,1) $$\beta _{1}=2$$ $$\beta _{1}=3$$ $$\beta _{1}=4$$ $$\beta _{1}=4$$ $$\beta _{1}=3$$ 5 12 65.6% 88.2% 95% 99.2% 98% 6 30 77.4% 95.2% 99% 99.8% 99.6% 7 54 74.4% 96.2% 99.2% 99.6% 99.8% 8 126 81.8% 97% 99.2% 100% 99.8% 9 240 85.2% 97.6% 99.4% 99% 98.2% 10 504 89.2% 98.2% 99.2% 99.4% 99% Table 5 Katsura-(n − 1) for the completeGraph strategy (#vertices, edge multiplicity) n BKK Bound (2,3) (2,4) (2,5) (3,2) (4,1) $$\beta _{1}=2$$ $$\beta _{1}=3$$ $$\beta _{1}=4$$ $$\beta _{1}=4$$ $$\beta _{1}=3$$ 5 12 65.6% 88.2% 95% 99.2% 98% 6 30 77.4% 95.2% 99% 99.8% 99.6% 7 54 74.4% 96.2% 99.2% 99.6% 99.8% 8 126 81.8% 97% 99.2% 100% 99.8% 9 240 85.2% 97.6% 99.4% 99% 98.2% 10 504 89.2% 98.2% 99.2% 99.4% 99% (#vertices, edge multiplicity) n BKK Bound (2,3) (2,4) (2,5) (3,2) (4,1) $$\beta _{1}=2$$ $$\beta _{1}=3$$ $$\beta _{1}=4$$ $$\beta _{1}=4$$ $$\beta _{1}=3$$ 5 12 65.6% 88.2% 95% 99.2% 98% 6 30 77.4% 95.2% 99% 99.8% 99.6% 7 54 74.4% 96.2% 99.2% 99.6% 99.8% 8 126 81.8% 97% 99.2% 100% 99.8% 9 240 85.2% 97.6% 99.4% 99% 98.2% 10 504 89.2% 98.2% 99.2% 99.4% 99% Table 6 Rounded expected probability of success assuming uniform distribuition of permutations and full monodromy group d $$\beta _{1}=2$$ $$\beta _{1}=3$$ $$\beta _{1}\geq 4$$ 12 90.5% 99.3% 100.0% 30 96.5% 99.9% 100.0% 54 98.1% 100.0% 100.0% 126 99.2% 100.0% 100.0% 240 99.6% 100.0% 100.0% 504 99.8% 100.0% 100.0% d $$\beta _{1}=2$$ $$\beta _{1}=3$$ $$\beta _{1}\geq 4$$ 12 90.5% 99.3% 100.0% 30 96.5% 99.9% 100.0% 54 98.1% 100.0% 100.0% 126 99.2% 100.0% 100.0% 240 99.6% 100.0% 100.0% 504 99.8% 100.0% 100.0% Table 6 Rounded expected probability of success assuming uniform distribuition of permutations and full monodromy group d $$\beta _{1}=2$$ $$\beta _{1}=3$$ $$\beta _{1}\geq 4$$ 12 90.5% 99.3% 100.0% 30 96.5% 99.9% 100.0% 54 98.1% 100.0% 100.0% 126 99.2% 100.0% 100.0% 240 99.6% 100.0% 100.0% 504 99.8% 100.0% 100.0% d $$\beta _{1}=2$$ $$\beta _{1}=3$$ $$\beta _{1}\geq 4$$ 12 90.5% 99.3% 100.0% 30 96.5% 99.9% 100.0% 54 98.1% 100.0% 100.0% 126 99.2% 100.0% 100.0% 240 99.6% 100.0% 100.0% 504 99.8% 100.0% 100.0% Table 7 Examples with solution count smaller than BKK bound (timings in seconds) problem wnt $${\mathop{SO}}(4)$$ $${\mathop{SO}}(5)$$ $${\mathop{SO}}(6)$$ $${\mathop{SO}}(7)$$ count 9 40 384 4768 111616 MonodromySolver 0.52 4 23 528 42791 Bertini 42 81 10605 out of memory PHCpack 862 103 > one day problem wnt $${\mathop{SO}}(4)$$ $${\mathop{SO}}(5)$$ $${\mathop{SO}}(6)$$ $${\mathop{SO}}(7)$$ count 9 40 384 4768 111616 MonodromySolver 0.52 4 23 528 42791 Bertini 42 81 10605 out of memory PHCpack 862 103 > one day Table 7 Examples with solution count smaller than BKK bound (timings in seconds) problem wnt $${\mathop{SO}}(4)$$ $${\mathop{SO}}(5)$$ $${\mathop{SO}}(6)$$ $${\mathop{SO}}(7)$$ count 9 40 384 4768 111616 MonodromySolver 0.52 4 23 528 42791 Bertini 42 81 10605 out of memory PHCpack 862 103 > one day problem wnt $${\mathop{SO}}(4)$$ $${\mathop{SO}}(5)$$ $${\mathop{SO}}(6)$$ $${\mathop{SO}}(7)$$ count 9 40 384 4768 111616 MonodromySolver 0.52 4 23 528 42791 Bertini 42 81 10605 out of memory PHCpack 862 103 > one day 6.3 Completion rate We investigate the completion rate of Algorithm 3.1 for the Katsura family parametrized by n with fixed support and generically chosen coefficients. Tables 4 and 5 contain the percentage of successes from 500 runs with distinct random seeds. In Table 6 we show the computed expected values using Remark 4.4. For $$\beta _{1}\geq 3$$ the observed success rates approach the expected values of Table 6. We note that the flower strategy is again closest to the estimates. We do not expect the numbers produced in experiments to match the numbers in Table 1, since the assumptions made for that statistical analysisare quite idealistic; however, both the analysis and experiments show that the probability of success approaches 100% rapidly as the number of solutions grows and the first Betti number increases. 6.4 Timings and comparison with other solvers All timings appearing in this section are done on one thread and on the same machine. Remarks 3.2 and 4.2 show that we should expect the number of tracked paths in Algorithms 3.1 and 3.4 to be linear (with a small constant!) in the number of solutions of the system. In this section we highlight the practicality of our approach in two ways. First, the monodromy method dramatically extends our computational ability for systems where the solution count turns out to be significantly smaller than the count corresponding to a more general family, for example, BKK count for sparse systems. This means that the existing blackbox methods, whose complexity relies on a larger count, are likely to spend significantly more time in computation compared to our approach. In Table 7 we collect timings on several challenging examples mentioned in recent literature where smaller solution counts are known, thus providing us with rigorous test cases for our heuristic stopping criterion. The first system in the table is that of the wnt signaling pathway reaction network mentioned in Section 6.2. The others come from the problem of computing the degree of $${\mathop{SO}}(n)$$, the special orthogonal group, as a variety (Brandt et al., 2017). Table 8 Software timings on large examples (in seconds) problem cyclic-10 cyclic-11 noon-10 BKK bound 35940 184756 59029 completeGraph(2,3) 610 7747 failed (107820 paths) (540155 paths) (59001 solutions) completeGraph(2,4) 740 8450 935 (129910 paths) (737432 paths) (236051 paths) PHCpack 538 4256 751 HOM4PS2 62 410 120 problem cyclic-10 cyclic-11 noon-10 BKK bound 35940 184756 59029 completeGraph(2,3) 610 7747 failed (107820 paths) (540155 paths) (59001 solutions) completeGraph(2,4) 740 8450 935 (129910 paths) (737432 paths) (236051 paths) PHCpack 538 4256 751 HOM4PS2 62 410 120 View Large Table 8 Software timings on large examples (in seconds) problem cyclic-10 cyclic-11 noon-10 BKK bound 35940 184756 59029 completeGraph(2,3) 610 7747 failed (107820 paths) (540155 paths) (59001 solutions) completeGraph(2,4) 740 8450 935 (129910 paths) (737432 paths) (236051 paths) PHCpack 538 4256 751 HOM4PS2 62 410 120 problem cyclic-10 cyclic-11 noon-10 BKK bound 35940 184756 59029 completeGraph(2,3) 610 7747 failed (107820 paths) (540155 paths) (59001 solutions) completeGraph(2,4) 740 8450 935 (129910 paths) (737432 paths) (236051 paths) PHCpack 538 4256 751 HOM4PS2 62 410 120 View Large Below is a list of comments on the set-up: For our implementation we chose small graphs with $$\beta _{1}\leq 4$$ and the random edge-selection strategy. The stopping criterion is ‘stabilization’ as discussed in Section 3.2.3. While the blackbox solver of PHCpack ultimately performs polyhedral homotopy continuation, Bertini relies by default on an equation-by-equation technique dubbed regeneration (see Bates et al., 2013). The latter may be faster than the former in certain cases, which this series of examples shows. Secondly, when the solution count is given by the BKK bound our method is a viable alternative to polyhedral homotopy solvers, since the number of paths we track is linear in the number of solutions. The timings on a few large benchmark problems of our current implementation and several other software packages are in Table 8. Our goal in the rest of this section is to show that our running times are in the same ballpark as polyhedral homotopies. Below is a list of comments on the set-up: For our implementation we chose two small graphs and default (random) edge-selection strategy. For PHCpack there is a way to launch a mixed volume computation, with the option of creating a system with the same support and random coefficients together with its solutions. This is the option we are using; the blackbox computation takes a little longer. HOM4PS2 (Lee et al., 2008) is not open source, unlike all other software mentioned here. (We use HOM4PS2 stock examples for all systems and call its blackbox polyhedral homotopies solver.) HOM4PS2 may use just-in-time compilation of straight-line programs used for evaluation, which speeds up computations considerably. (PHCpack does not use this technique; neither does our software, but our preliminary experiments in Macaulay2 show a potential for a 10- to 20-fold speed up over our currently reported timings.) Remark 6.2 For large examples assuming the probabilistic model leading to Theorem 4.1 and Remark 4.2, the probability of success should be extremely close to 100% even for a random graph with $$\beta _{1} = 2$$. The run of noon-10, which is an example of neural network model from the study by Noonburg (1989), demonstrates an unlikely, but possible failure for $$\beta _{1} = 2$$ followed by success at $$\beta _{1} = 3$$. On the examples in Table 8 we also ran the blackbox solvers of Bertini and Numerical Algebraic Geometry (Leykin, 2011), which use the total-degree homotopy. Both were able to finish noon-10 with timings similar to the table, but all other problems took longer than a day. This is expected, as the BKK bound of noon-10 is only slightly sharper than the Bézout bound. Remark 6.3 In comparison with the naive dynamic strategy (Section 3.1) our framework loses slightly only in one aspect: memory consumption. For a problem with d solutions the naive approach stores up to (and typically close to) 2d points. The number of points our approach stores is up to (and typically considerably fewer than) d times the number of vertices. For instance, it is up to 4d points in all runs in Table 8. The number of tracked paths is significantly lower in our framework: for example, the naive strategy tracks about 7500 paths on average for cyclic 7. Even before looking at Table 2 it is clear that running the flower strategy in combination with the incremental dynamic strategy of Section 3.3 guarantees to dominate the naive strategy. 7. Generalizations While we propose a more general algorithmic framework, a concurrent goal of this paper is to demonstrate that significant practical advantages are already apparent when we apply a relatively simple implementation and analysis to simple problems (linearly parametrized families). The following topics thus lie outside the scope of this article, but seem deserving of further study: One advantage of the MS approach is that it can tolerate numerical failures of the underlying homotopy tracker. In fact, we already implemented a simple failure resistant mechanism, and it successfully tolerates a few failures that arise in some runs for large test examples in Section 6.4. A natural extension of this paper’s statistical analysis would be to model the algorithm’s performance in the presence of failures. Ideally, heuristics such as edge potentials should incorporate information such as the failures discussed above. It is also of interest to adapt potentials to the parallel setting discussed below. The parallelization of the MS approach is not as straightforward as that of other homotopy continuation methods. The question of when speed-ups close to linear can be achieved should be addressed. Consider the generalized set-up in which the base space B is an irreducible variety and the family is given by a rational map from P into a space of systems. To apply our general framework, a major requirement is to find an effective way to parametrize a curve between two points of P. This parametrization would conceivably depend on the nature of the problem being considered. Certain other ingredients are also likely to be problem-specific—for instance, even in the case of $$P = {\mathbb{C}}^{m},$$ the construction of the initial seed $$(p_{0},x_{0})$$ is complicated by the possibility that the systems’ coefficients are nonlinear in the parameters. Nonetheless, this is one of the strengths of the MS framework—once all required ‘oracles’ are supplied, the procedures become effective. In the classical language of enumerative geometry, the monodromy groups we consider are isomorphic to Galois groups of incidence varieties (essentially solution varieties in our terminology). For a large class of Schubert problems and other interesting incidence varieties the associated Galois group turns out to be the full symmetric group (Leykin & Sottile, 2009). A suitable modification of our dynamic strategy is one practical approach to verifying this in conjectural cases. Our paper demonstrates the strength of our method relative to other techniques such as polyhedral homotopy and regeneration. Building on our framework one could use polyhedral homotopy as a subroutine to quickly populate a partial solution set (quickly discarding any path that becomes poorly conditioned). Further advantages may be achievable by using different techniques in parallel. These and other hybrid approaches have the potential to produce even faster and more robust blackbox solvers. Funding Research of T.D., C.H., K.L., and A.L. is supported in part by National Science Foundation grants DMS-1151297 and DMS-1719968. References Babai , L. ( 1989 ) The probability of generating the symmetric group . J. Comb. Theory Ser. A , 52 , 148 – 153 . Google Scholar CrossRef Search ADS Bates , D. J. , Hauenstein , J. D. , Sommese , A. J. & Wampler , C. W. ( 2013 ) Numerically Solving Polynomial Systems With Bertini , vol. 25. Philadelphia : SIAM . Bernstein , D. N. ( 1975 ) The number of roots of a system of equations . Funk. Anal. Priložen , 9 , 1 – 4 . Blum , L. , Cucker , F. , Shub , M. & Smale , S. ( 1998 ) Complexity and Real Computation . New York : Springer . Google Scholar CrossRef Search ADS Brandt , M. , Bruce , J. , Brysiewicz , T. , Krone , R. & Robeva , E. ( 2017 ) The degree of SO($$n,\mathbb{C}$$) . Combinatorial Algebraic Geometry (G. Smith & B. Sturmfels eds). Fields Institute Communications, vol. 80 . New York City : Springer , pp. 229 – 246 . Chen , J. & Kileel , J. ( 2016 ) Numerical implicitization for Macaulay2 . arXiv preprint arXiv:1610.03034 . del Campo , A. M. & Rodriguez , J. I. ( 2017 ) Critical points via monodromy and local methods . J. Symbolic Comput. , 79 , 559 – 574 . Google Scholar CrossRef Search ADS Dixon , J. D. ( 1969 ) The probability of generating the symmetric group . Math. Z. , 110 , 199 – 205 . Google Scholar CrossRef Search ADS Duff , T. , Hill , C. , Jensen , A. , Lee , K. , Leykin , A. & Sommars , J. ( 2016 ) MonodromySolver: a Macaulay2 package for solving polynomial systems via homotopy continuation and monodromy . Available at http://people.math.gatech.edu/aleykin3/MonodromySolver. Emiris , I. Z. & Vidunas , R. ( 2014 ) Root counts of semi-mixed systems, and an application to counting nash equilibria . Proceedings of the 39th International Symposium on Symbolic and Algebraic Computation, ISSAC ’14 . New York, NY, USA : Association for Computing Machinery , pp. 154 – 161 . Galligo , A. & Miclo , L. ( 2012 ) On the cut-off phenomenon for the transitivity of randomly generated subgroups . Random Struct. Alg. , 40 , 182 – 219 . Google Scholar CrossRef Search ADS Galligo , A. & Poteaux , A. ( 2011 ) Computing monodromy via continuation methods on random Riemann surfaces . Theor. Comput. Sci. , 412 , 1492 – 1507 . Google Scholar CrossRef Search ADS Grayson , D. R. & Stillman , M. E. Macaulay2, a software system for research in algebraic geometry . Available at http://www.math.uiuc.edu/Macaulay2/. Gross , E. , Harrington , H. A. , Rosen , Z. & Sturmfels , B. ( 2016 ) Algebraic systems biology: a case study for the Wnt pathway . Bull. Math. Biol. , 78 , 21 – 51 . Google Scholar CrossRef Search ADS PubMed Hauenstein , J. D. & Rodriguez , J. I. ( 2015 ) Multiprojective witness sets and a trace test . arXiv preprint arXiv:1507.07069 . Hauenstein , J. D. , Rodriguez , J. I. & Sottile , F. ( 2017 ) Numerical computation of Galois groups . Found. Comput. Math . To appear . Hauenstein , J. D. & Sottile , F. ( 2012 ) Algorithm 921: alphaCertified: certifying solutions to polynomial systems . ACM Trans. Math. Softw. , 38 , 28 . Google Scholar CrossRef Search ADS Huber , B. & Sturmfels , B. ( 1995 ) A polyhedral method for solving sparse polynomial systems . Math. Comp. , 64 , 1541 – 1555 . Google Scholar CrossRef Search ADS Jensen , A. N. ( 2016 ) An implementation of exact mixed volume computation . Mathematical Software - ICMS 2016 - 5th International Conference, Berlin, Germany, July 11–14, 2016, Proceedings . Cham, Switzerland : Springer , pp. 198 – 205 . Lee , T.-L. , Li , T.-Y. & Tsai , C.-H. ( 2008 ) HOM4PS-2.0: a software package for solving polynomial systems by the polyhedral homotopy continuation method . Computing , 83 , 109 – 133 . Google Scholar CrossRef Search ADS Leykin , A. ( 2011 ) Numerical algebraic geometry . J. Softw. Alg. Geom. , 3 , 5 – 10 . Google Scholar CrossRef Search ADS Leykin , A. , Rodriguez , J. I. & Sottile , F. ( 2016 ) Trace test . arXiv preprint arXiv:1608.00540 . Leykin , A. & Sottile , F. ( 2009 ) Galois groups of Schubert problems via homotopy computation . Math. Comp. , 78 , 1749 – 1765 . Google Scholar CrossRef Search ADS Leykin , A. & Verschelde , J. ( 2009 ) Decomposing solution sets of polynomial systems: a new parallel monodromy breakup algorithm . Int. J. Comput. Sci. Eng. , 4 , 94 – 101 . Google Scholar CrossRef Search ADS MacLean , A. L. , Rosen , Z. , Byrne , H. M. & Harrington , H. A. ( 2015 ) Parameter-free methods distinguish Wnt pathway models and guide design of experiment . Proc. Natl. Acad. Sci. , 112 , 2652 – 2657 . Google Scholar CrossRef Search ADS Malajovich , G. ( 2017 ) Computing mixed volume and all mixed cells in quermassintegral time . Found. Comput. Math. , 17 , 1293 – 1334 . Google Scholar CrossRef Search ADS Morgan , A. ( 1987 ) Solving Polynomial Systems Using Continuation for Engineering and Scientific Problems . Englewood Cliffs, NJ : Prentice Hall Inc . Noonburg , V. W. ( 1989 ) A neural network modeled by an adaptive Lotka–Volterra system . SIAM J. Appl. Math. , 49 , 1779 – 1792 . Google Scholar CrossRef Search ADS Sommese , A. J. , Verschelde , J. & Wampler , C. W. ( 2001 ) Using monodromy to decompose solution sets of polynomial systems into irreducible components . Application of Algebraic Geometry to Coding Theory, Physics, and Computation , Dordrecht : Springer Netherlands , pp. 297 – 315 . Google Scholar CrossRef Search ADS Sommese , A. J. , Verschelde , J. & Wampler , C. W. ( 2002 ) Symmetric functions applied to decomposing solution sets of polynomial systems . SIAM J. Numer. Anal. , 40 , 2026 – 2046 . Google Scholar CrossRef Search ADS Sommese , A. J. , Verschelde , J. & Wampler , C. W. ( 2005 ) Introduction to numerical algebraic geometry . Solving Polynomial Equations: Foundations, Algorithms, and Applications (A. Dickenstein & I. Z. Emiris eds). Berlin, Heidelberg : Springer Berlin Heidelberg , pp. 301 – 337 . Sommese , A. J. & Wampler II , C. W. ( 2005 ) The Numerical Solution of Systems of Polynomials . Hackensack, NJ : World Scientific Publishing Co. Pte. Ltd . Google Scholar CrossRef Search ADS Verschelde , J. ( 1999 ) Algorithm 795: PHCpack: a general-purpose solver for polynomial systems by homotopy continuation . ACM Trans. Math. Softw. , 25 , 251 – 276 . Google Scholar CrossRef Search ADS Verschelde , J. , Verlinden , P. & Cools , R. ( 1994 ) Homotopies exploiting Newton polytopes for solving sparse polynomial systems . SIAM J. Numer. Anal. , 31 , 915 – 930 . Google Scholar CrossRef Search ADS © The Author(s) 2018. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices) For permissions, please e-mail: journals. permissions@oup.com

### Journal

IMA Journal of Numerical AnalysisOxford University Press

Published: Apr 13, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations