Analysis of Normal-Form Algorithms for Solving Systems of Polynomial Equations

TL;DR

The paper investigates normal-form, eigenvalue-based rootfinding methods for multivariate polynomial systems by analyzing two Macaulay-based reduction strategies (direct Macaulay and null-space Macaulay) and practical speedups (degree-by-degree construction, random combinations). It provides a comprehensive temporal-complexity framework, comparing QR and SVD variants, and reveals a fundamental stability hurdle when many roots are close, exemplified by the Noferini–Townsend devastator example where eigenvalue conditioning scales as $\kappa(\lambda,M_g)=\Omega(\varepsilon^{-n})$ despite the root conditioning being $\kappa(z,f)=\varepsilon^{-1}$. Despite these instabilities, the methods perform well for well-separated, low-dimensional root sets, with the degree-by-degree and SVD variants offering favorable speed-accuracy tradeoffs in practical regimes. The work highlights the need for preconditioning or basis-rescaling strategies to mitigate conditioning problems and outlines future work integrating these methods with Chebyshev-based proxy approaches for real-root computation in higher dimensions. Overall, it provides a rigorous performance map for these rootfinding techniques and practical guidance on when to apply which variant.

Abstract

We examine several of the normal-form multivariate polynomial rootfinding methods of Telen, Mourrain, and Van Barel and some variants of those methods. We analyze the performance of these variants in terms of their asymptotic temporal complexity as well as speed and accuracy on a wide range of numerical experiments. All variants of the algorithm are problematic for systems in which many roots are very close together. We analyze performance on one such system in detail, namely the 'devastating example' that Noferini and Townsend used to demonstrate instability of resultant-based methods.

Figures (8)

• Figure 1: Comparison of the FLOPs used with a simple (no speedups) null space reduction versus the degree-by-degree null space reduction. The top and bottom rows show results in dimensions three and four, respectively. The panels on the left show the total number of FLOPs required for the two constructions. The panels on the right show the ratio of operations between the simple construction and the degree-by-degree construction. Notice that the peak, where the savings of degree-by-degree is most significant, moves to the right as dimension increases.
• Figure 2: Log-scale plots comparing the Direct Macaulay QRP and SVD methods in terms of average residuals (left panel) and Macaulay condition numbers (right panel) for random, dense polynomials (in the power basis) of total degree $2$ over dimensions $3$ through $7$. Both the average residuals and the average eigenvalue condition number for the SVD method are smaller than for the QRP method. Similar improvements of SVD over QRP can be observed for a fixed dimension and varying degree.
• Figure 3: Average solve time for random, dense, polynomial systems using the SVD variants of direct Macaulay, Macaulay null space reduction, and the degree-by-degree method. On the left, we have quadratic systems of varying dimension. On the right, we have Dimension 3 polynomial systems of varying degree. Although the QRP method generally performs faster, the SVD method gives more accurate results (see Figure \ref{['fig:QRPvsSVD_res/eigs']}). The results for systems of dimensions and degrees not visible on these plots are solved quickly enough to make their inclusion unhelpful for comparison using a linear scale.
• Figure 4: Average residuals for random, dense, polynomial systems using the SVD variants of direct Macaulay, Macaulay null space reduction, and the degree-by-degree methods. On the left, we have quadratic systems of varying dimension. On the right, we have Dimension 3 polynomial systems of varying degree. Comparing the residuals of QRP yields a similar trend, but the SVD consistently outperforms QRP in terms of accuracy.
• Figure 5: Numerically calculated conditioning ratios for devastating and random quadratic systems solved using the direct Macaulay SVD method. The conditioning ratios of random systems show very slow exponential growth in dimension compared to the devastating example. Orange: Line of best fit for conditioning ratios of $n$ dimensional devastating systems with a randomly chosen $Q$. All of the computed conditioning ratios were within $0.015\%$ of the theoretical value, and the computed growth rate was $9.001$. Blue: Random systems of quadratic polynomials with coefficients drawn from the standard normal distribution. The violin and box plots show the distributions of the conditioning ratios of theses systems. The dotted black lines represent the tail ends of these distributions out to the most extreme observed conditioning ratios. The line of best fit to the base-10 logarithm of the conditioning ratios is also shown, with a growth rate of $g \approx 0.102$.

Theorems & Definitions (37)

• Proposition 4.1
• proof
• Definition 5.1
• Lemma 5.2
• proof
• Lemma 5.3
• proof
• Lemma 6.1
• proof
• Lemma 6.2