# Polynomial Systems, Homotopy Continuation, and Applications

Timothy Duff
Margaret Regan

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 whose solutions are known, called the start system, can be continuously deformed into a system of equations whose solutions we would like to know, called the target system. The following example illustrates some of the key ideas.

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

where are variables representing unknown quantities and 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 , the system of equations has a finite number of solutions which are nonsingular in the sense that , the Jacobian matrix of with respect to , is invertible. The number 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 and find nonsingular solutions to the system . This can be accomplished with a straight-line homotopy as in Example 1. This has the advantage that the solutions of 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 for any sufficiently general choice of . For this task, we utilize the parameter homotopy

for all and some fixed . In particular, has known solutions , computed in the ab initio phase, and one aims to compute the solutions to . For generic values of the constant , the arc for connects to and avoids the complex discriminant locus. Thus, for , defines precisely solution paths connecting the  points in with the solutions to .

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 and are fixed pivots and , , and are the lengths of the moving links.

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 passes through nine generic positions. Alt’s problem amounts to solving a system with complex solutions. These solutions carry the extra structure of a -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 ) 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 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).

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 , and consider two discrete probability distributions,

If we draw samples from each distribution, , , we may record the frequency of all possible outcomes into a matrix of counts . Let be the matrix giving the joint distribution . Our statistical model is the algebraic variety , where is the variety of matrices of rank at most , and is the probability simplex,

Note that random variables and are independent iff . The MLE problem amounts to minimizing the log-likelihood function

Calculating the partial derivatives reveals that they are rational functions in and , and this leads to a polynomial system of equations whose solutions are the critical points of 3 restricted to the model . 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 , , and are tabulated below (table adapted from 12.) Do you see any patterns?

### 2.3. Power flow systems

Let be an integer and consider a finite undirected graph on the vertices . Fix , , and consider the system of equations in unknowns

This is one formulation of the power flow equations used to model a network of agents (also known as buses.) Each bus may represent a power station, customer, or some other entity within an electrical grid. The coefficients for are called susceptances and are assumed to be known. The unknowns are the real and imaginary parts of the voltage at the -th bus. The fixed values 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 “trivial” solutions obtained by fixing all non-reference buses for . The number of “non-trivial” solutions turns out to depend on the topology of the graph . At one extreme, for the complete graph , these equations will have solutions over the complex numbers, as long as the susceptances are sufficiently generic. For the -cycle , the generic complex root count is a more modest . In either case, there is a symmetry on these solutions that sends , allowing us to reduce the cost of homotopy path-tracking by a factor of . Of course, only the solutions where all and are real are of any practical interest. Figure 4 illustrates the distribution of real solutions as a subset of the susceptances vary for the -cycle . The visible regions which are blue, red, green, purple, and yellow correspond to parameter values with and solutions, respectively. We refer to the article 18 for detailed explanations and many other interesting experimental results obtained using homotopy continuation methods.

## 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.

## References

[1]
A. Baskar, S. Bandyopadhyay. An algorithm to compute the finite roots of large systems of polynomial equations arising in kinematic synthesis. Mechanism and Machine Theory 133 (2019), 493–513.
[2]
A. Baskar, C. Hills, M. M. Plecnik, and J. D. Hauenstein. Estimating the complete solution set of the approximate path synthesis problem for four-bar linkages using random monodromy loops. Proceedings of the ASME 2022 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference. (2022)
[3]
Daniel J. Bates, Jonathan D. Hauenstein, and Andrew J. Sommese, Efficient path tracking methods, Numer. Algorithms 58 (2011), no. 4, 451–459, DOI 10.1007/s11075-011-9463-8. MR2854200Show rawAMSref\bib{PathTrackMethods}{article}{ author={Bates, Daniel J.}, author={Hauenstein, Jonathan D.}, author={Sommese, Andrew J.}, title={Efficient path tracking methods}, journal={Numer. Algorithms}, volume={58}, date={2011}, number={4}, pages={451--459}, issn={1017-1398}, review={\MR {2854200}}, doi={10.1007/s11075-011-9463-8}, } Close amsref.
[4]
D. J. Bates, J. D. Hauenstein, A. J. Sommese, and C. W. Wampler. Bertini: Software for Numerical Algebraic Geometry. Available at bertini.nd.edu with permanent doi: dx.doi.org/10.7274/R0H41PB5.
[5]
Daniel J. Bates, Jonathan D. Hauenstein, Andrew J. Sommese, and Charles W. Wampler, Numerically solving polynomial systems with Bertini, Software, Environments, and Tools, vol. 25, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2013. MR3155500Show rawAMSref\bib{BHSW13}{book}{ author={Bates, Daniel J.}, author={Hauenstein, Jonathan D.}, author={Sommese, Andrew J.}, author={Wampler, Charles W.}, title={Numerically solving polynomial systems with Bertini}, series={Software, Environments, and Tools}, volume={25}, publisher={Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA}, date={2013}, pages={xx+352}, isbn={978-1-611972-69-6}, review={\MR {3155500}}, } Close amsref.
[6]
P. Breiding and S. Timme. HomotopyContinuation. jl: A package for homotopy continuation in Julia. International Congress on Mathematical Software. (2018).
[7]
S. Deshpande, A. Purwar. A task-driven approach to optimal synthesis of planar four-bar linkages for extended Burmester problem. Journal of Mechanisms and Robotics 9 (2017), no. 6, 061005.
[8]
A. K. Dhingra, J. C. Cheng, D. Kohli. Synthesis of six-link, slider-crank and four-link mechanisms for function, path and motion generation using homotopy with m-homogenization. Journal of Mechanical Design 116 (1994), no. 4, 1122–1131.
[9]
P. B. Edwards, A. Baskar, C. Hills, M Plecnik, and J. D. Hauenstein. Output mode switching for parallel five-bar manipulators using a graph-based path planner. In preparation. (2022)
[10]
D. Grayson and M. Stillman. Macaulay2, a software system for research in algebraic geometry. Available at http://www.math.uiuc.edu/Macaulay2/.
[11]
Jonathan D. Hauenstein and Margaret H. Regan, Adaptive strategies for solving parameterized systems using homotopy continuation, Appl. Math. Comput. 332 (2018), 19–34, DOI 10.1016/j.amc.2018.03.028. MR3788668Show rawAMSref\bib{HauensteinRegan}{article}{ author={Hauenstein, Jonathan D.}, author={Regan, Margaret H.}, title={Adaptive strategies for solving parameterized systems using homotopy continuation}, journal={Appl. Math. Comput.}, volume={332}, date={2018}, pages={19--34}, issn={0096-3003}, review={\MR {3788668}}, doi={10.1016/j.amc.2018.03.028}, } Close amsref.
[12]
Jonathan Hauenstein, Jose Israel Rodriguez, and Bernd Sturmfels, Maximum likelihood for matrices with rank constraints, J. Algebr. Stat. 5 (2014), no. 1, 18–38, DOI 10.18409/jas.v5i1.23. MR3279952Show rawAMSref\bib{HRS}{article}{ author={Hauenstein, Jonathan}, author={Rodriguez, Jose Israel}, author={Sturmfels, Bernd}, title={Maximum likelihood for matrices with rank constraints}, journal={J. Algebr. Stat.}, volume={5}, date={2014}, number={1}, pages={18--38}, review={\MR {3279952}}, doi={10.18409/jas.v5i1.23}, } Close amsref.
[13]
Jonathan D. Hauenstein and Frank Sottile, Algorithm 921: alphaCertified: certifying solutions to polynomial systems, ACM Trans. Math. Software 38 (2012), no. 4, Art. 28, 20, DOI 10.1145/2331130.2331136. MR2972672Show rawAMSref\bib{alphaCertified}{article}{ author={Hauenstein, Jonathan D.}, author={Sottile, Frank}, title={Algorithm 921: alphaCertified: certifying solutions to polynomial systems}, journal={ACM Trans. Math. Software}, volume={38}, date={2012}, number={4}, pages={Art. 28, 20}, issn={0098-3500}, review={\MR {2972672}}, doi={10.1145/2331130.2331136}, } Close amsref.
[14]
Jonathan D. Hauenstein and Charles W. Wampler, Isosingular sets and deflation, Found. Comput. Math. 13 (2013), no. 3, 371–403, DOI 10.1007/s10208-013-9147-y. MR3047005Show rawAMSref\bib{Isosingular}{article}{ author={Hauenstein, Jonathan D.}, author={Wampler, Charles W.}, title={Isosingular sets and deflation}, journal={Found. Comput. Math.}, volume={13}, date={2013}, number={3}, pages={371--403}, issn={1615-3375}, review={\MR {3047005}}, doi={10.1007/s10208-013-9147-y}, } Close amsref.
[15]
P. Hruby, T. Duff, A. Leykin, and T. Pajdla, T. Learning to Solve Hard Minimal Problems. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (pp. 5532-5542). (2022)
[16]
Anton Leykin, Jan Verschelde, and Ailing Zhao, Newton’s method with deflation for isolated singularities of polynomial systems, Theoret. Comput. Sci. 359 (2006), no. 1-3, 111–122, DOI 10.1016/j.tcs.2006.02.018. MR2251604Show rawAMSref\bib{Deflation}{article}{ author={Leykin, Anton}, author={Verschelde, Jan}, author={Zhao, Ailing}, title={Newton's method with deflation for isolated singularities of polynomial systems}, journal={Theoret. Comput. Sci.}, volume={359}, date={2006}, number={1-3}, pages={111--122}, issn={0304-3975}, review={\MR {2251604}}, doi={10.1016/j.tcs.2006.02.018}, } Close amsref.
[17]
Anton Leykin, Numerical algebraic geometry, J. Softw. Algebra Geom. 3 (2011), 5–10, DOI 10.2140/jsag.2011.3.5. MR2881262Show rawAMSref\bib{NAG4M2}{article}{ author={Leykin, Anton}, title={Numerical algebraic geometry}, journal={J. Softw. Algebra Geom.}, volume={3}, date={2011}, pages={5--10}, issn={1948-7916}, review={\MR {2881262}}, doi={10.2140/jsag.2011.3.5}, } Close amsref.
[18]
J. Lindberg, A. Zachariah, N. Boston, and B. Lesieutre. The Distribution of the Number of Real Solutions to the Power Flow Equations. IEEE Transactions on Power Systems. (2022).
[19]
M. M. Plecnik, R. S. Fearing. Finding only finite roots to large kinematic synthesis systems. Journal of Mechanisms and Robotics 9 (2017), no. 2, 021005.
[20]
M. M. Plecnik, J. M. McCarthy. Computational design of Stephenson II six-bar function generators for 11 accuracy points. Journal of Mechanisms and Robotics 8 (2016), no. 1, 011017.
[21]
M. M. Plecnik, J. M. McCarthy. Kinematic synthesis of Stephenson III six-bar function generators. Mechanism and Machine Theory 97 (2016), 112–126.
[22]
B. Roth, F. Freudenstein. Synthesis of path-generating mechanisms by numerical methods. Journal of Engineering for Industry 85 (1963), no. 3, 298–304.
[23]
Andrew J. Sommese, Jan Verschelde, and Charles W. Wampler, Introduction to numerical algebraic geometry, Solving polynomial equations, Algorithms Comput. Math., vol. 14, Springer, Berlin, 2005, pp. 301–335, DOI 10.1007/3-540-27357-3_8. MR2161992Show rawAMSref\bib{SVW}{article}{ author={Sommese, Andrew J.}, author={Verschelde, Jan}, author={Wampler, Charles W.}, title={Introduction to numerical algebraic geometry}, conference={ title={Solving polynomial equations}, }, book={ series={Algorithms Comput. Math.}, volume={14}, publisher={Springer, Berlin}, }, date={2005}, pages={301--335}, review={\MR {2161992}}, doi={10.1007/3-540-27357-3\_8}, } Close amsref.
[24]
Andrew J. Sommese and Charles W. Wampler II, The numerical solution of systems of polynomials, World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ, 2005. Arising in engineering and science, DOI 10.1142/9789812567727. MR2160078Show rawAMSref\bib{SW05}{book}{ author={Sommese, Andrew J.}, author={Wampler, Charles W., II}, title={The numerical solution of systems of polynomials}, note={Arising in engineering and science}, publisher={World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ}, date={2005}, pages={xxii+401}, isbn={981-256-184-6}, review={\MR {2160078}}, doi={10.1142/9789812567727}, } Close amsref.
[25]
Seth Sullivant, Algebraic statistics, Graduate Studies in Mathematics, vol. 194, American Mathematical Society, Providence, RI, 2018, DOI 10.1090/gsm/194. MR3838364Show rawAMSref\bib{Sul}{book}{ author={Sullivant, Seth}, title={Algebraic statistics}, series={Graduate Studies in Mathematics}, volume={194}, publisher={American Mathematical Society, Providence, RI}, date={2018}, pages={xiii+490}, isbn={978-1-4704-3517-2}, review={\MR {3838364}}, doi={10.1090/gsm/194}, } Close amsref.
[26]
I. Ullah and S. Kota. Optimal synthesis of mechanisms for path generation using Fourier descriptors and global search methods. Journal of Mechanical Design 119 (1997), no. 4, 504–510.
[27]
C. W. Wampler, A. Morgan, A. J. Sommese. Complete solution of the nine-point path synthesis problem for four-bar linkages. Journal of Mechanical Design 114 (1992), no. 1, 153–159.

## Credits

Figure 1 is courtesy of Silviana Amethyst.

Figures 2 and 3 are courtesy of Jonathan D. Hauenstein.

Figure 4 is courtesy of Julia Lindberg.