Homotopy Techniques for Analytic Combinatorics in Several Variables

From the beginning of its modern period in work of Pemantle and Wilson [21], the goal of ACSV has always been to develop methods explicit enough to be implemented in a computer algebra system. The ‘surgery’ approach of [21], which applies to generating functions with smooth singular sets that form manifolds, essentially computes a residue in one variable to obtain a $(d-1)$-dimensional integral that is approximated using the saddle-point method. Although this surgery method does not require much theory beyond univariate analytic combinatorics, it requires strong conditions on the locations of minimal points that can be computationally expensive to verify. Later techniques, using cones of hyperbolicity [3] and multivariate residue and homology computations [2], rely on more advanced theory but simplify the assumptions that need to be verified for the results to hold. In the simplest cases, which hold for the majority of examples encountered in combinatorial applications, it suffices to determine which of the critical points are minimal and then add explicit asymptotic contributions corresponding to the (finite number of) minimal critical points. The most expensive step in such an analysis is almost always checking minimality.

The first systematic algorithmic study of ACSV methods was conducted by Melczer and Salvy [18], who encoded critical points using a symbolic-numeric data structure known as a Kronecker or rational univariate representation and then reduced checking minimality to rigorously approximating the roots of certain univariate polynomials to sufficiently high accuracy. Those authors created a preliminary implementation of their work, which does not certify numeric computations to provide rigorous proofs and requires combinatorial rational functions, in the Maple computer algebra system. A rational function $F(\mathbf{z})$ is combinatorial if all of its power series coefficients are non-negative: although this condition is satisfied for any multivariate generating function, in many combinatorial examples only one diagonal of $F$ enumerates a combinatorial class and the non-diagonal entries have negative coefficients. It is an open problem, even in the univariate case, whether it is decidable to detect when a rational function is combinatorial (see [20] for some open problems in this area). Although Melczer and Salvy [18] detail a method that, in principle, yields an algorithm for asymptotics that does not require combinatorality, in practice an implementation in Maple would not halt in reasonable time beyond low degree examples in two or three variables.

Instead of continuing with the Kronecker representation approach of Melczer and Salvy, in this paper we exploit homotopy continuation methods to certify minimality of critical points, and ultimately determine asymptotics of $\mathbf{r}$-diagonals of rational functions. Using the HomotopyContinuation.jl Julia package [7] for polynomial system solving, we provide the first implementation of ACSV methods under assumptions that often hold in practice. Our implementation is efficient enough to work even without the assumption of combinatorality, although when the user knows a priori that their input rational function is combinatorial then the computation is greatly reduced. In addition, we describe two heuristic methods to classify minimal critical points using numerical approximations that are extremely efficient, and are the only implemented algorithms we currently know of that can aid in the search for minimal points in more than three variables.

The rest of this paper proceeds as follows. Section 1 gives a quick recap of the methods of ACSV and the high-level problems that need to be decided to find asymptotics, with a description of numerical algebraic geometry methods for polynomial system solving given in Section 2. Section 3 uses this background material to detail our ACSVHomotopy.jl Julia package, while Section 4 illustrates the package on a wide variety of combinatorial examples, including benchmarks between different algorithms. Although our algorithms always terminate, due to the nature of homotopy continuation methods they may not always provide a rigorous proof of asymptotics -- Section 5 discusses this issue and describes situations in which the algorithms do give rigorous proofs. Finally, Section 6 concludes with some extensions that we believe should be addressed next.

The hardest work in applying Theorem 1.4 is computing the critical points, defined implicitly by (1), and determining which, if any, are minimal.

(Combinatorial Case) Recall that a function is called combinatorial if its power series expansion contains only a finite number of negative coefficients. When $F$ is combinatorial there is a simple test for minimal critical points.

In the combinatorial case we firstly use Lemma 1.6 to characterize the minimal critical points with positive coordinates by studying the solutions to (2). From them, find the solutions to (1) with the same coordinate-wise modulus. The following algorithm summarizes this approach.

(General Case) If $F$ is not combinatorial, or if we don’t know a priori that $F$ is combinatorial, then it is no longer sufficient to consider only the critical points with positive real coordinates to check minimality. In order to express the moduli of coordinates as algebraic equations, we write $H(\mathbf{x}+i\mathbf{y})=H^{\mathfrak{R}}(\mathbf{x},\mathbf{y})+iH^{\mathfrak{I}}(\mathbf{x},\mathbf{y})$ for real variables $\mathbf{x},\mathbf{y}\in\mathbb{R}^{d}$ and polynomials $H^{\mathfrak{R}},H^{\mathfrak{I}}\in\mathbb{R}[\mathbf{x},\mathbf{y}]$. Translating the smooth critical point equations (1) into these new coordinates gives that $\mathbf{z}=\mathbf{a}+i\mathbf{b}$ with $\mathbf{a},\mathbf{b}\in\mathbb{R}^{d}$ is critical if and only if

for some $\lambda_{R},\lambda_{I}\in R$, where $1\leq j\leq d$ in each equation. To test minimality of these critical points we add the equations

for $1\leq j\leq d$, and verify there is no real solution to (3)-(7) with $0<t<1$. Generically (3)-(7) have a finite set of real solutions, corresponding to the generically finite number of critical points of $F$, but because this system contains $3d+4$ equations in $4d+3$ variables it will never have a (non-zero) finite number of solutions over the complex numbers. By considering critical values of the projection map onto the $t$ coordinate, Melczer and Salvy [18] proved that minimality can be tested by adding the additional equations

for $1\leq j\leq d$. When $\nu_{1}\neq 0$ then we can scale by $\nu_{1}$ and introduce the equations

to (3)-(7), resulting in a square system with $4d+4$ variables and equations. The case when $\nu_{1}=0$ is dealt with separately by adding the equations

for $1\leq j\leq d$. We determine the minimal critical points by finding $\mathbf{p}+i\mathbf{q}$ such that equations (3)-(5) have a real solution with $(\mathbf{a},\mathbf{b})=(\mathbf{p},\mathbf{q})$ but neither (3)-(8) nor (3)-(8’) have a real solution with $(\mathbf{a},\mathbf{b})=(\mathbf{p},\mathbf{q})$ and $0<t<1$. This process provides the following algorithm.

Unfortunately, to moving $4d+4$ variables makes verifying minimality much less practical than the combinatorial case. In essence, Lemma 1.6 states that to prove minimality in the combinatorial case it is sufficient to consider specific line segments in $\mathbb{R}^{d}$, while to prove minimality in the general case one must consider a much larger set of points in $\mathbb{C}^{d}$ whose coordinate-wise moduli lie on specific line segments in $\mathbb{R}^{d}$.

Homotopy continuation is one method to find numerical approximations of solutions to an $n\times n$ square system $\mathcal{F}=(f_{1},\dots,f_{n})$ of polynomial equations with $n$ variables. From the system $\mathcal{F}$ we construct an $n\times n$ polynomial system $\mathcal{G}$ whose solutions are known a priori. The system $\mathcal{G}$ is called a start system and the system $\mathcal{F}$ is called the target system: connecting $\mathcal{F}$ and $\mathcal{G}$ using a homotopy $\mathcal{H}(x,t)$ such that $\mathcal{H}(x,0)=\mathcal{G}$ and $\mathcal{H}(x,1)=\mathcal{F}$, we obtain solutions of $\mathcal{F}$ by tracking homotopy paths from $t=0$ to $t=1$. To track the homotopy paths, a numerical ordinary differential equation solving technique called the Davidenko equation and Newton iteration are used. These tracking techniques are typically referred to as predictor-corrector methods. For details, see [24, Chapter 2]. Homotopy continuation is implemented in Bertini [4], HomotopyContinuation.jl [7], and NAG4M2 [16].

The complexity of solving a polynomial system using homotopy continuation is determined by the number of homotopy paths to track. It is thus important to track a number of paths that are at least as large as the number of solutions of the system (so that all solutions can be found) but is not too much larger (to save computation). For polytopes $Q_{1},\dots,Q_{n}$ the Euclidean volume $\text{Vol}(a_{1}Q_{1}+\cdots+a_{n}Q_{n})$ of the Minkowski sum $a_{1}Q_{1}+\cdots+a_{n}Q_{n}$ is a homogeneous polynomial in the $n$ variables $a_{1},\dots,a_{n}$, whose coefficient of $a_{1}a_{2}\cdots a_{n}$ is the mixed volume $\text{MVol}(Q_{1},\dots,Q_{n})$ of $Q_{1},\dots,Q_{n}$.

The polyhedral homotopy continuation method established by Huber and Sturmfels [13] is one common way to construct a start system whose solutions form a set with the size of the mixed volume of a system. Consider a polynomial

where $A$ is a collection of integer lattice points. Multiplying each monomial $\mathbf{x}^{\mathbf{a}}$ of $f$ by some term $t^{w(\mathbf{a})}$ for a lifting function $w:A\rightarrow\mathbb{Z}$, we obtain the lifted polynomial

Suppose that a target system $\mathcal{F}$ consists of polynomials $f_{1},\dots,f_{n}$ supported on $A_{f_{1}},\dots,A_{f_{n}}$, respectively. Lifting all polynomials $f_{1},\dots,f_{n}$ in $\mathcal{F}$ gives a lifted system $\overline{\mathcal{F}}(\mathbf{x},t)$ satisfying $\overline{\mathcal{F}}(\mathbf{x},1)=\mathcal{F}$. The solutions of $\overline{\mathcal{F}}$ can be expressed by Puiseux series $\mathbf{x}(t)=(x_{1}(t),\dots,x_{n}(t))$ where

for some $\alpha_{i}\in\mathbb{Q}$ and nonzero constant $y_{i}$, and substituting $\mathbf{x}(t)$ back into our polynomials gives

For a suitable choice of $w$, the constants $\mathbf{y}$ and exponents $\mathbf{a}$ can be computed at each branch of $\overline{\mathcal{F}}$, ultimately describing a start system $\mathcal{G}=\overline{\mathcal{F}}(\mathbf{x},0)$ with the right number of solutions. The polyhedral homotopy continuation is implemented in HOM4PS2 [15], HomotopyContinuation.jl [7], and PHCpack [26].

As seen in Section 1, we typically have some solutions of a polynomial system representing critical points and want to determine additional solutions to rule out those that are non-minimal. This ‘bootstrapping’ can be accomplished by monodromy.

For $m,n\in\mathbb{N}$, consider the complex linear space of $n\times n$ square systems $\mathcal{F}_{p}=(f_{p}^{1},\dots,f_{p}^{n})$ depending on some coefficient parameters $p\in\mathbb{C}^{m}$, where the monomial support for each polynomial $f_{p}^{i}$ is fixed. If we consider an affine linear map $\varphi:p\mapsto F_{p}$ for $p\in\mathbb{C}^{m}$ then we can write $\varphi(\mathbb{C}^{m})=B$, where $B$ is a parametrized linear variety of systems, and we define the solution variety $V=\{(F_{p},x)\in B\times\mathbb{C}^{n}\mid F_{p}(x)=0\}$ and projection map $\pi:V\rightarrow B$.

Assume that the fiber $\pi^{-1}(F_{p})$ only has finitely many points for a generic choice of $p$. The set $D$ of systems in $B$ with non-generic fiber is called the branch locus of $\pi$. Each element in the fundamental group $\pi_{1}(B\setminus D)$ of loops in $B\setminus D$ modulo homotopy equivalence induces a permutation on the fiber $\pi^{-1}(F_{p})$, which is called a monodromy action. To find all solutions of a system $F_{p}\in\pi(V)$ with generic $p$, one can first find a seed solution $(p_{0},x_{0})\in V$ and numerically compute the monodromy action to find all solutions of $F_{p}$. When the solution variety $V$ is irreducible then the monodromy action is transitive. This method for finding solutions of polynomial systems is studied and implemented in [7, 9].

By construction, numerical methods return approximations, so some kind of certification is necessary for rigorous results. Specifically, a user needs a certificate that an approximation obtained by the homotopy method is properly approximating a solution of a system. A numerical approximation is called certified if it can be refined to an actual solution of the system to an arbitrary precision by applying iterative operators (such as Newton iteration). Software providing such certification includes alphaCertified [12], the function certify implemented in HomotopyContinuation.jl [6] and NumericalCertification [14]. In our implementation, we use the function certify in HomotopyContinuation.jl exploiting the Krawczyk’s method via interval arithmetic [19, Chapter 8].

For the combinatorial case we first compute the distinct solutions to (1) using a polyhedral homotopy with certification by Krawczyk’s method. We then solve and certify (2) with the added equation $(1-t)\mu-1=0$ to eliminate all solutions with $t=1$ (there are never any solutions with $t=0$ as this would imply $H(\mathbf{0})=0$, contradicting $F$ having a power series expansion). Since we no longer have solutions where $t=1$, by refining the solutions to sufficient precision we can determine the solutions with positive real coordinates where $0<t<1$, match the projection onto the $\mathbf{z}$ variables of each to a distinct solution of (1), and thus rule out all non-minimal critical points with positive coordinates. We then find all critical points with the same coordinate-wise moduli and return that set.

As in the work of Melczer and Salvy, the most expensive operation occurs when trying to group roots with the same coordinate-wise modulus as a known minimal critical point (step (4) in Algorithm 1). When using a symbolic-numeric method it is possible to compute minimal polynomials for the values of the coordinates and use this to identify points with the same coordinate-modulus by computing numerically to $O(h\delta^{3d})$ bits of precision, where $h$ is a bound on the bitsize of the coefficients of the denominator $H$ and $\delta$ is the degree of $H$ (see [18, Corollary 54]). Combining past bounds in the literature [23, 8] allows us to identify an explicit precision such that if two solutions have coordinate-wise moduli agreeing to this precision then their coordinate-wise moduli are exactly equal. Our bound is also of $O(h\delta^{3d})$ bits, however the constant in front is worse than that when using minimal polynomials. After a minimal critical point is identified, our code continually refines precision using Newton iteration until any points with the same coordinate-wise moduli are found. In practice, any points with different coordinate-wise moduli are identified at precision much lower than the worst case bound, but we must compute up to the bound when there are points with the same coordinate-wise moduli. Unfortunately, the extreme precision involved with a large number of variables means that we cannot always rigorously check to the required accuracy, and the code returns a warning to the user with its output when this occurs.

In general we must consider the extended systems (3)-(8) and (3)-(8’), which essentially doubles the number of variables under consideration. Mirroring the combinatorial case, we can solve (3)-(8) and (3)-(8’) using a polyhedral homotopy, certify the results using Krawczyk’s method, and refine to a sufficient precision to determine when $0<t<1$ to rule out non-minimal points.

As seen in the examples of Section 4, the high number of variables in the extended system even in low dimensional examples means it does not terminate within reasonable time for polynomials with four or more variables. In order to speed up our solvers, we can numerically approximate the distinct solutions to the small system (3)-(5) and then substitute each of these solutions as parameters into the extended equations (6)-(8) and (6)-(8’). In the implementation, it is done by running the function find_min_crits with the flag approx_crit = true.

It is also possible to use the approximations to the critical points as a start system to solve (6)-(8) using the monodromy method. More precisely, for any $(\mathbf{a},\mathbf{b})$ solving (3)-(5) we set $(\mathbf{x},\mathbf{y})=(\mathbf{a},\mathbf{b})$ and $t=1$ in (6)-(8) and then compute a corresponding start value of $\nu$ by computing the left kernel of the Jacobian matrix of (6)-(7) with respect to variables $\mathbf{x},\mathbf{y}$ and $t$. From a given parameter value $(\mathbf{a},\mathbf{b})$ and the initial solution $(\mathbf{x},\mathbf{y},\nu,t)$, we collect real solutions from the monodromy method and check if $t\in(0,1)$: if it is then we remove the parameter value $(\mathbf{a},\mathbf{b})$ as it is non-minimal. Interestingly, it appears that monodromy cannot detect solutions where $\nu=0$ when starting with a non-zero value of $\nu$, and vice-versa, suggesting that solution variety is the union of components corresponding to these cases. We thus repeat this process separately for the cases where $\nu=0$ and when (8’) replaces (8). Finally, we return the values of $(\mathbf{a},\mathbf{b})$ that are not disregarded.

We believe that further study of the geometric properties of the extended system (6)-(8) could help make this monodromy approach a powerful tool for ACSV analysis.