Introduction. The search for a robust, global, and numerically stable algorithm for solving multivariate polynomial systems has persisted for decades [2, 3, 5, 6, 7, 8, 9]. The goal is to compute all the solutions to zero-dimensional polynomial systems of the form:

\begin{equation} \begin{pmatrix} p_1(x_1, \ldots , x_d) \\ \vdots \\ p_d(x_1, \ldots , x_d) \end {pmatrix} = \mathbf {0}, \label {eq:polysystem} \end{equation}

where \(d \geq 2\) and \(p_1, \ldots , p_d\) are polynomials in \(x_1,\ldots , x_d\) with complex coefficients. All the promising solvers are based on an elegant approach of converting the multivariate rootfinding problem in (1) into one or more eigenproblems At first this approach appears to be a practitioner’s dream as a difficult rootfinding problem can be solved by the robust QR or QZ algorithm. For this reason, these methods have received considerable research attention from the scientific computing community. However, we are currently stuck waiting for new ideas to emerge from algebraic geometry or hoping for novel structured eigensolvers from numerical linear algebra.

A popular class of techniques known as hidden variable resultant methods are notoriously difficult—and maybe impossible—to make numerically robust [6]. Naive implementations are plagued with unwanted spurious solutions, inaccurate roots, and miss zeros. In this talk, we will discuss the ongoing quest for a numerically stable rootfinder using Sylvester resultants, Gröbner bases, Möller–Stetter matrices, and Macaulay resultants. Our focus will be on understanding the sources of the instability in these approaches, in the hope that they can be circumvented.

Motivation for Eigenvalue-Based Approaches. Given a well-conditioned rootfinding problem, we would like to derive a stable algorithm to solve it. Roughly speaking, an algorithm is stable if it computes an accurate solution to well-conditioned problems. The search for a stable algorithm for multivariate polynomial rootfinding is motivated by the existence of stable algorithms for many related problems. All the univariate problems, such as eigenproblems, univariate polynomial rootfinding, and matrix polynomial eigenproblems, have stable algorithms. Likewise, there are stable algorithms to solve linear systems of the form \(A \mathbf {x} = \mathbf {b}\), which are multivariate.

For univariate rootfinding, instead of solving a rootfinding problem directly, one can first construct an eigenproblem whose eigenvalues match the desired roots. The companion matrix of a polynomial \(p(x)\) is an example of this, as its characteristic polynomial is \(p\), so its eigenvalues are the roots of \(p\). One can solve the companion eigenproblem using an eigensolver, which is one of the most reliable algorithms in numerical linear algebra. For roots in \([-1,1]\), a provably stable algorithm for univariate polynomial rootfinding is based on the colleague matrix [4]. For multivariable rootfinding, we attempt the same conversion, i.e., we try to convert (1) into one or more generalized eigenvalue problems (GEPs). For polynomial systems in (1) in \(d\) variables, one usually constructs \(d\) GEPs, the eigenvalues of which give the coordinates of each root. The Macaulay resultant method is an exception as it constructs a single GEP and extracts the roots from the eigenvectors, not eigenvalues. Analogously to the univariate case, one hopes that the eigenproblems are as well-conditioned as the original rootfinding problem. Unfortunately, this is not always the case.

Gröbner bases, Möller–Stetter matrices, rational univariate reductions, multiparameter eigenvalue problems, and Macaulay resultants all convert (1) into one or more GEPs, either directly or by way of a univariate rootfinding problem. For each method, we show that either a constructed eigenproblem or an intermediate univariate rootfinding problem can be more ill-conditioned than the original rootfinding problem by a factor that is exponentially large in \(d\).

A devastating example. The analysis of the instability of Sylvester and Cayley resultant method appears in [5, 6] and studied the following “devastating” polynomial rootfinding problem.

By a conditioning analysis, one should expect to find the root at \((0,\ldots ,0)\) to within \(\approx \sigma u\), where \(u\) is the unit roundoff on a computer. However, it has been shown that Sylvester resultants can only achieve \(\approx \sigma ^{-2} u\) when \(d=2\) and Cayley can only achieve \(\approx \sigma ^{-d} u\) [6]. The eigenproblems constructed by these methods can be far more sensitive to perturbations than the original rootfinding problem, which is a hallmark of an unstable algorithm. Similar examples show that Gröbner bases, Möller–Stetter matrices, rational univariate reductions, multiparameter eigenvalue problems, and Macaulay resultants can also generate highly ill-conditioned eigenproblems.

Theoretical results supported by numerical experiments. One might be hopeful that these worst-case examples are rarely and never realized in practice. Unfortunately, this is not the case in practice. We regularly observe the ill-conditioning in numerical experiments, causing solvers to generate spurious solutions, miss roots, or compute them inaccurately. There is a small glimmer of hope in some new approaches that construct structured eigenproblems and solve them with eigensolvers that respective that structure [1].

We hope that our work inspires new approaches that circumvent the limitations of current techniques and provide robust solutions to the fundamental multivariate polynomial roofinding problem.

References

• [1]  H. He, D. Kressner, and B. Plestenjak. Randomized methods for computing joint eigenvalues, with applications to multiparameter eigenvalue problems and root finding. arXiv preprint arXiv:2409.00500, 2024.

• [2]  G. F. Jónsson and S. A. Vavasis. Accurate solution of polynomial equations using Macaulay resultant matrices. Math. Comp., 74(249):221–262, 2005.

• [3]  B. Mourrain, S. Telen, and M. Van Barel. Truncated normal forms for solving polynomial systems: generalized and efficient algorithms. J. Symbolic Comput., 102:63–85, 2021.

• [4]  Y. Nakatsukasa and V. Noferini. On the stability of computing polynomial roots via confederate linearizations. Math. Comp., 85(301):2391–2425, 2016.

• [5]  Y. Nakatsukasa, V. Noferini, and A. Townsend. Computing the common zeros of two bivariate functions via Bézout resultants. Numer. Math., 129(1):181–209, 2015.

• [6]  V. Noferini and A. Townsend. Numerical instability of resultant methods for multidimensional rootfinding. SIAM J. Numer. Anal., 54(2):719–743, 2016.

• [7]  B. Plestenjak and M. E. Hochstenbach. Roots of bivariate polynomial systems via determinantal representations. SIAM J. Sci. Comput., 38(2):A765–A788, 2016.

• [8]  F. Rouillier. Solving zero-dimensional systems through the rational univariate representation. Appl. Algebra Engrg. Comm. Comput., 9(5):433–461, 1999.

• [9]  S. Telen, B. Mourrain, and M. Van Barel. Solving polynomial systems via truncated normal forms. SIAM J. Matrix Anal. Appl., 39(3):1421–1447, 2018.