Polynomial Systems, Homotopy Continuation, and Applications

Systems of multivariate polynomial equations are ubiquitous throughout mathematics. They also appear prominently in scientific applications such as kinematics 2022, computer vision 1115, power flow engineering 18, and statistics 12. Numerical homotopy continuation methods are a fundamental tool for both solving these systems and determining more refined information about their structure.

In this article, we offer a brief glimpse of polynomial homotopy continuation methods: the general theory, a few applications, and some software packages that implement these methods. Our aim is to spark the reader’s interest in this exciting and broad area of research. We invite those looking to learn more to join us at the AMS Short Course: Polynomial systems, homotopy continuation, and applications, to be held January 2–3 at the 2023 Joint Mathematics Meetings in Boston.

1. Homotopy Continuation

Many types of homotopy continuation methods exist, but they all are based on the same strategy. A system of equations $g(z)=0$ whose solutions are known, called the start system, can be continuously deformed into a system of equations $f(z) = 0$ whose solutions we would like to know, called the target system. The following example illustrates some of the key ideas.

Figure 1.

Graphical summary of homotopy continuation.

Example 1.

Consider the target polynomial system

$$\begin{equation*} f(z) = \left[\begin{array}{c} 4 (z_1^2 z_2 + z_1 z_2) - 3 z_2\\z_2^2 - 2 z_1 + 1 \end{array}\right]. \end{equation*}$$

The corresponding total degree start system

$$\begin{equation*} g(z) = \left[\begin{array}{c} z_1^3 -1 \\z_2^2 -1 \end{array}\right] \end{equation*}$$

has 6 solutions of the form $(z_1, z_2) = (\omega , \pm 1)$, where $\omega$ is a third root of unity. These will serve as the smooth start points depicted in Figure 1.

Using the straight-line homotopy

$$\begin{equation*} H(z; t) = (1-t)f(z) + t g(z) = 0 \end{equation*}$$

to deform $g$ into $f$, each solution for $t\approx 1$ may be estimated from the start points by numerical predictor/corrector methods, also known as homotopy path-tracking 24, Ch. 2. Iterating this procedure, we find for this example there are six smooth homotopy paths of the form $z(t) : (0, 1] \to \mathbb{C}^2$, where $z(1)$ is one of the smooth start points. As $t\to 0^+$, one solution path will diverge towards infinity, and three will come together at the singular endpoint $(z_1, z_2) = (1/2, 0)$. The remaining two endpoints $(z_1, z_2) = (-1/2, \pm 2 i)$ are finite and nonsingular. We see that $f$ has three solutions, five counting multiplicity.

Path-tracking methods are well-studied in numerical analysis, and are especially potent when applied to a parametrized polynomial system

$$\begin{equation} \begin{split} f(z;p) &= f(z_1,\dots ,z_N;p_1,\dots ,p_k)\\ &= \left[\begin{array}{c} f_1(z_1,\dots ,z_N;p_1,\dots ,p_k) \\\vdots \\f_n(z_1,\dots ,z_N;p_1,\dots ,p_k) \end{array}\right] \end{split} {\tag{1}} \end{equation}$$

where $z_1, \ldots , z_N$ are variables representing unknown quantities and $p_1, \ldots , p_k$ are parameters representing physical measurements. Section 2 gives few examples of such systems appearing in applications.

A general parameter continuation theorem 24, Theorem 7.1.1 is based on the fact that for almost all parameter values $p^* \in \mathbb{C}^k$, the system of equations $f(z; p^*)=0$ has a finite number $d$ of solutions $z^* \in \mathbb{C}^N$ which are nonsingular in the sense that $J_z f(z^*;p^*)$, the $N\times N$ Jacobian matrix of $f$ with respect to $z$, is invertible. The number $d$ is sometimes called the generic root count of the system 1.

The essential observation of the parameter continuation theorem is that all isolated solutions can be computed via tailor-made homotopies which operate in a problem’s natural parameter space. These parameter homotopies involve two phases, summarized below. See 5, Chap. 6 for more details.

Ab initio phase

The first step for a parameter homotopy is to fix parameter values $p^*\in \mathbb{C}^k$ and find $d$ nonsingular solutions to the system $f(z;p^*) = 0$. This can be accomplished with a straight-line homotopy as in Example 1. This has the advantage that the solutions of $H(z;1) = g(z)=0$ are trivial to compute. More sophisticated methods allow us to track potentially fewer paths in this phase. These include multihomogeneous homotopies 24, Sec. 8.4.2, polyhedral homotopy 24, Sec. 8.5, and methods based on monodromy 24, Sec. 15.4.

Parameter homotopy phase

With the ab initio phase complete, the “online” parameter homotopy phase aims to solve $f(z;p) = 0$ for any sufficiently general choice of $p\in \mathbb{C}^k$. For this task, we utilize the parameter homotopy

$$\begin{multline} H(z;t) = f(z;\tau (t)\cdot p^* + (1-\tau (t))\cdot p) = 0\\ \text{where}\quad \tau (t) = \dfrac{\gamma t}{1+(\gamma -1)t} {\tag{2}} \end{multline}$$

for all $t\in [0,1]$ and some fixed $\gamma \in \mathbb{C}$. In particular, $H(z;1) = f(z;p^*)=0$ has known solutions $S$, computed in the ab initio phase, and one aims to compute the solutions to $H(z;0) = f(z;p)=0$. For generic values of the constant $\gamma \in \mathbb{C}$, the arc $\tau (t)\cdot p^* + (1-\tau (t))\cdot p$ for $t\in [0,1]$ connects $p^*$ to $p$ and avoids the complex discriminant locus. Thus, for $t\in [0,1]$, $H(z;t) =0$ defines precisely $d$ solution paths connecting the $d$ points in $S$ with the $d$ solutions to $f(z;p)=0$.

Our discussion of homotopy continuation methods in this section is necessarily incomplete. Here we list a few additional topics falling under the rubric of general methods. One important topic is numerical algebraic geometry 23, which allows us to study positive-dimensional algebraic varieties. In the opposite case of an overdetermined system, several techniques allow us to reduce to the case of a well-constrained parametrized system of the form 1; see 11 and the references therein. Lastly, we mention numerical certification methods which can prove that approximated solutions will converge to exact solutions (see, e.g., 13), and deflation methods for regularizing systems with singular solutions 1416.

2. Applications

Polynomial homotopy continuation has been a key to advances in various applications. We summarize three that will be featured in our short course.

2.1. Kinematics

Mechanical linkage systems of interest have constrained motions that are naturally modeled with systems of polynomial equations. Such polynomial formulations cover a wide breadth of mechanisms including planar, spherical, and spatial types.

For example, consider the 4-bar mechanism in Figure 2 where $A$ and $B$ are fixed pivots and $\ell _1$, $\ell _2$, and $\ell _3$ are the lengths of the moving links.

Figure 2.

Schematic of a four-bar linkage.

The path synthesis problem associated with this mechanism seeks to find all possible linkages that meet certain design requirements. Wampler et al. 27 solved the exact path synthesis problem for 4-bars, also known as Alt’s problem, which imposes that the coupler trace point $X$ passes through nine generic positions. Alt’s problem amounts to solving a system with $8652$ complex solutions. These solutions carry the extra structure of a $6$-fold symmetry group, including the 3-fold Roberts’ cognate triplets.

Homotopy continuation has been used for these exact synthesis problems by finding the roots of the corresponding system of polynomial equations 1819202122. These are large-scale (up to $10^6$) root-finding problems, where homotopy continuation is the only method capable of computing complete solution sets at such scales.

Other methods focus on approximate kinematic synthesis, relying on optimization techniques to accommodate any number of design specifications. For example, in 2, the approximate path synthesis problem using optimization yields about 303,249 $\pm$ 713 Roberts’ cognate triples as critical points with 95$\%$ confidence for the 4-bar linkage. This offers an advantage over exact methods, with the downside being large computational effort. Numerical homotopy continuation methods were central to the use of monodromy loops that made this computation possible.

The use of homotopy continuation within optimization problems in kinematic design has also enabled the study of the configurations of the parallel 5-bar mechanism, which displays more nonlinearity that the serial 5-bar mechanism. Figure 3 shows this complicated configuration space. In 9, homotopy continuation is used to quantify transmission quality using the curves of input and output singularities. This enables developing a path that switches between non-neighboring output modes (i.e., solution sheets).

Figure 3.

Workspace curve of input singularities and velocity ellipse for a five-bar linkage.

In general, homotopy continuation methods have led to the analysis and solving of much more complicated problems in kinematics.

2.2. Algebraic statistics

Maximum likelihood estimation (MLE) is a fundamental technique of statistical inference, in which the likelihood function associated to a data set is maximized over a space of all possible parameters that specifies a statistical model. A major theme in the field of algebraic statistics 25 is that the space of model parameters will often be an algebraic variety. In this case, homotopy continuation methods can be used for global optimization of the likelihood function. This complements the widely used EM algorithm for MLE, which has the advantage of being easy to implement, but is generally susceptible to local minima.

To make these ideas expressed above concrete, we consider a discrete statistical model from 12. Fix positive integers $m,n\ge 1$, and consider two discrete probability distributions,

$$\begin{align*} p_X &: \{ 1, \ldots , m \} \to [0,1],\\ p_Y &: \{ 1 , \ldots , n \} \to [0,1]. \end{align*}$$

If we draw $k$ samples from each distribution, $X_1, \ldots , X_k \sim p_X$, $Y_1, \ldots , Y_k \sim p_Y$, we may record the frequency of all possible $mn$ outcomes into a $m\times n$ matrix of counts $U$. Let $P$ be the $m\times n$ matrix giving the joint distribution $p_{i j} = \Pr \left(X =i, \, Y=j \right)$. Our statistical model is the algebraic variety $\mathcal{V}_r \cap \Delta _{mn-1}$, where $\mathcal{V}_r$ is the variety of matrices of rank at most $r$, and $\Delta _{mn-1}$ is the probability simplex,

$$\begin{equation*} \Delta _{mn-1} = \{ P \in \mathbb{R}^{m\times n} \mid p_{i j} \ge 0, \, \displaystyle \sum p_{i j} = 1 \}. \end{equation*}$$

Note that random variables $X$ and $Y$ are independent iff $P\in \mathcal{V}_r \cap \Delta _{mn-1}$. The MLE problem amounts to minimizing the log-likelihood function

$$\begin{equation} \ell (P ; U) = \log \left( \displaystyle \prod _{\substack{1\le i \le m\\1\le j \le n}} p_{i j}^{u_{i j}} \right) - k\, \log \left( \left(\displaystyle \sum _{\substack{1\le i \le m\\1\le j \le n}}p_{i j}\right) \right). {\tag{3}} \end{equation}$$

Calculating the partial derivatives $\frac{\partial \ell }{\partial p_{i j}}$ reveals that they are rational functions in $u_{i j}$ and $p_{i j}$, and this leads to a polynomial system of equations whose solutions are the critical points of 3 restricted to the model $\mathcal{V}_r \cap \Delta _{mn-1}$. Among these critical points is the maximum-likelihood estimate. The total number of critical points is known as the ML degree of the model. Using parameter homotopies, we can track exactly this number of homotopy paths to find all critical points. The ML degrees for small $m$, $n$, and $r$ are tabulated below (table adapted from 12.) Do you see any patterns?

$$\begin{equation*} \begin{matrix} & (m,n) = & (3,3) & (3,4) & (3,5) & (4,4) & (4,5) & (4,6)\\ r=1 & & 1 & 1 & 1 & 1 & 1 & 1 \\ r=2 & & 10 & 26 & 58 & 191 & 843 & 3119 \\ r=3 & & 1 &1 &1 & 191 & 843 & 3119 \\ r=4 & & & & &1 &1 & 1 \end{matrix} \end{equation*}$$

2.3. Power flow systems

Let $n\ge 2$ be an integer and consider a finite undirected graph $G$ on the vertices $\{ 1, \ldots , n \}$. Fix $x_1 = 1$, $y_1=0$, and consider the system of $2(n-1)$ equations in $2(n-1)$ unknowns

$$\begin{equation} \begin{split} x_i^2 + y_i^2 - 1 &= 0 \quad i = 2, \, \ldots , \, n,\\ \displaystyle \sum _{(i,j)\in E(G)} b_{i j} (x_i y_j - x_j y_i) &= 0 \quad i = 2, \, \ldots , \, n. \displaystyle \end{split} {\tag{4}} \end{equation}$$

This is one formulation of the power flow equations used to model a network of $n$ agents (also known as buses.) Each bus may represent a power station, customer, or some other entity within an electrical grid. The coefficients $b_{i j} \in \mathbb{R}$ for $(i,j) \in E(G)$ are called susceptances and are assumed to be known. The unknowns $(x_i, y_i)$ are the real and imaginary parts of the voltage at the $i$-th bus. The fixed values $(x_1, y_1)=(1,0)$ determine the reference bus.

Solving the power flow equations plays an important role in operating and controlling electrical networks. It is common for engineers to approach this problem with local, iterative algorithms such as Newton’s method, which will return a single real solution.

But what can we say about all solutions to 4? Notice that there are $2^{n-1}$ “trivial” solutions obtained by fixing all non-reference buses $(x_i, y_i) = (\pm 1, 0)$ for $i=2, \ldots , n$. The number of “non-trivial” solutions turns out to depend on the topology of the graph $G$. At one extreme, for the complete graph $G=K_n$, these equations will have $\binom{2n-2}{n-1}$ solutions over the complex numbers, as long as the susceptances are sufficiently generic. For the $n$-cycle $G=C_n$, the generic complex root count is a more modest $n \cdot \binom{n-1}{\lfloor (n-1)/2 \rfloor }$. In either case, there is a symmetry on these solutions that sends $(x_i, y_i) \mapsto (x_i, - y_i)$, allowing us to reduce the cost of homotopy path-tracking by a factor of $2$. Of course, only the solutions where all $x_i$ and $y_i$ are real are of any practical interest. Figure 4 illustrates the distribution of real solutions as a subset of the susceptances vary for the $5$-cycle $C_5$. The visible regions which are blue, red, green, purple, and yellow correspond to parameter values with $0, 2, 4, 6,$ and $8$ solutions, respectively. We refer to the article 18 for detailed explanations and many other interesting experimental results obtained using homotopy continuation methods.

Figure 4.

Distribution of real solutions to equations 4 for the cycle graph $G=C_5$, with three susceptances $b_{i j}$ drawn uniformly from the unit sphere in $\mathbb{R}^{2}$ and the other two fixed.

3. Software

A wide variety of software packages implementing polynomial homotopy continuation methods exists. Here we highlight three that will be used during the upcoming short course:

1.

Bertini 4 is a standalone software package, whose functionality includes many of the standard homotopy methods for isolated solutions, as well as numerical irreducible decomposition for positive-dimensional solutions.

2.

HomotopyContinuation.jl 6 is a software package written for the Julia language, a programming language designed for high-performance numerical computing.

3.

Macaulay2 10 is a computer algebra system focused on computational commutative algebra and algebraic geometry. Primarily a tool for symbolic computation, it also has a growing number of numerical algebraic geometry packages 17.