Logo
ScienceStack
Table of Contents
Fetching ...
Sign In

Numerical Instability of Algebraic Rootfinding Methods

Emil Graf, Alex Townsend

Abstract

We demonstrate that the most popular variants of all common algebraic multidimensional rootfinding algorithms are unstable by analyzing the conditioning of subproblems that are constructed at intermediate steps. In particular, we give multidimensional polynomial systems for which the conditioning of a subproblem can be worse than the conditioning of the original problem by a factor that grows exponentially with the number of variables.

Numerical Instability of Algebraic Rootfinding Methods

Abstract

We demonstrate that the most popular variants of all common algebraic multidimensional rootfinding algorithms are unstable by analyzing the conditioning of subproblems that are constructed at intermediate steps. In particular, we give multidimensional polynomial systems for which the conditioning of a subproblem can be worse than the conditioning of the original problem by a factor that grows exponentially with the number of variables.
Paper Structure (20 sections, 13 theorems, 57 equations, 5 figures, 2 tables)

This paper contains 20 sections, 13 theorems, 57 equations, 5 figures, 2 tables.

Table of Contents

  1. Introduction
  2. Motivation for Eigenvalue-Based Approaches
  3. Background
  4. Generalized Eigenproblems and Matrix Polynomials
  5. Condition Numbers
  6. Zero-Dimensional Ideals and Varieties
  7. Gröbner Basis Elimination
  8. Rational Univariate Representation
  9. A Devastating Example for the Rational Univariate Representation
  10. Multiparameter Eigenproblem Solved by Operator Determinants
  11. Conditioning of the Operator Determinants Method
  12. Normal Form Methods
  13. Conditioning of the Möller--Stetter Eigenproblems
  14. Macaulay Resultant
  15. Conditioning of the Macaulay Resultant Eigenproblem
  16. ...and 5 more sections

Key Result

Theorem 4.2

Consider solving the system in ex:Hypercube via the rational univariate representation using a linear projection $t(x_1,\ldots,x_d) = u_1x_1 + \cdots + u_dx_d$, with $\sum_{i=1}^d \lvert u_i \rvert^2 \leq 1$, and let $f(x) = \prod_{i=1}^r (x-t(\mathbf{x}_i))$. Then, the root $x^* = t(\mathbf{x}^*)$

Figures (5)

  • Figure 1: Performance on a bivariate version of the system in \ref{['ex:dev']} for all methods except for Gröbner basis elimination (\ref{['ex:GBdev']}) and the rational univariate representation (\ref{['ex:Hypercube']}), with the root shifted to $(\frac{1}{3},\frac{1}{3})$. The deviation of the practical performance from our theory is explained by the extreme proximity of the roots when $\sigma$ is very small, which indicates that we should not expect conditioning analysis to be a good predictor. The methods are extremely inaccurate when $\sigma$ is small.
  • Figure 1: Performance of Gröbner basis elimination on \ref{['ex:GBdev']} for $d \geq 2$ and $\sigma = \frac{1}{2}$. We plot the practical performance against the theoretical performance of a stable algorithm and the prediction given by the analysis of \ref{['ex:GBdev']}. The deviation of the practical performance from our theory is explained by the extreme proximity of the roots for large values of $d$.
  • Figure 1: Performance of the multiparameter eigenproblem method on \ref{['ex:MEPdev']} for $d \geq 2$ and $\sigma = \frac{1}{100}$. We plot the practical performance against the theoretical performance of a stable algorithm and the prediction given by \ref{['thm:multipareig']}.
  • Figure 1: Performance of a normal form method and the Macaulay resultant method on \ref{['ex:notdev']}. We plot the practical performance against the theoretical performance of a stable algorithm. The deviation of observations from our prediction is because of the proximity of roots for small values of $\sigma$.
  • Figure 2: Performance of a normal form method on \ref{['ex:3Ddev']}. We plot the practical performance against the theoretical performance of a stable algorithm, the predicted performance given by the line $\sigma^{-2}$, and the prediction given by the Jacobian, which demonstrates that this example can be solved more accurately than would be predicted by a direct analogy of \ref{['thm:LowDimMS']}.

Theorems & Definitions (34)

  • Example 1.1
  • Example 3.1
  • Example 4.1
  • Theorem 4.2
  • Proof 1
  • Proposition 5.1
  • Proof 2
  • Theorem 5.2
  • Proof 3
  • Example 5.3
  • ...and 24 more