6.2 Numerical methods for polynomial system solving

Numerical methods for polynomial system solving are powerful tools for finding solutions to complex equations. These techniques, like homotopy continuation and Newton's method, use clever mathematical tricks to transform difficult problems into more manageable ones.

By tracking solution paths or iteratively refining guesses, these methods can find all solutions, even complex or multiple ones. Understanding their strengths and limitations is crucial for tackling real-world problems in engineering, physics, and other fields that involve polynomial systems.

Fundamentals of homotopy continuation

Basic concepts and formulation

• Homotopy continuation is a numerical method for solving systems of polynomial equations by continuously deforming a simple start system into the target system
• The basic idea is to construct a homotopy $H(x,t) = (1-t)F(x) + tG(x)$, where $F(x)$ is the target system, $G(x)$ is a start system with known solutions, and $t$ is a parameter varying from 0 to 1
  • As $t$ varies from 0 to 1, the solutions of $H(x,t) = 0$ trace out solution paths that connect the solutions of $G(x) = 0$ to those of $F(x) = 0$
  • The solution paths are tracked numerically using predictor-corrector methods, such as Euler's method for prediction and Newton's method for correction
• Homotopy continuation can be used to find all isolated solutions of a polynomial system, including complex solutions and multiple solutions

Start systems and solution paths

• The start system $G(x)$ is chosen to be easily solvable and to have at least as many solutions as the target system $F(x)$
  • Common choices include the total degree homotopy and the polyhedral homotopy
    • Total degree homotopy: $G(x)$ is a system of polynomials with the same degrees as $F(x)$ but with random coefficients
    • Polyhedral homotopy: $G(x)$ is a system of monomials determined by the Newton polytopes of $F(x)$
• The solution paths defined by the homotopy $H(x,t)$ have some important properties
  • The paths are smooth curves parameterized by $t \in [0,1]$
  • The paths do not intersect except possibly at the end points ($t=0$ or $t=1$)
  • The number of paths is equal to the number of solutions of the start system $G(x)$
• The gamma trick is a technique for handling singular solutions by slightly perturbing the homotopy to avoid singularities
  • The homotopy is modified as $H(x,t) = (1-t)F(x) + \gamma t(1-t)G(x)$, where $\gamma$ is a small random complex number
  • The gamma trick ensures that the solution paths avoid singular points and reach the target solutions

Newton's method for polynomial systems

Newton's iteration and its variants

• Newton's method is an iterative algorithm for finding roots of a system of equations, including polynomial systems
• The basic idea is to start with an initial guess and iteratively refine it by solving a linear approximation of the system at the current point
  • For a polynomial system $F(x) = 0$, the Newton iteration is $x_{k+1} = x_k - J(x_k)^{-1} F(x_k)$, where $J(x_k)$ is the Jacobian matrix of $F$ at $x_k$
  • The convergence of Newton's method is quadratic for simple roots but may be slower for multiple roots or near singular solutions
• Variants of Newton's method can improve convergence for multiple roots and singular solutions
  • The modified Newton's method uses a higher-order approximation of the system to achieve better convergence for multiple roots
    • The iteration is modified as $x_{k+1} = x_k - (I - P(x_k))J(x_k)^{-1} F(x_k)$, where $P(x_k)$ is a projection matrix that depends on the multiplicity structure of the root
  • The deflation method successively eliminates found solutions to improve the conditioning of the system and the convergence of Newton's method
    • After finding a solution $x^$, the system is deflated by replacing $F(x)$ with $(x - x^)^{-1} F(x)$ and solving the deflated system for the remaining solutions

Globalization strategies

• Globalization strategies can be used to ensure convergence of Newton's method from poor initial guesses
• Line search methods choose the step size in the Newton direction to ensure sufficient decrease in the norm of the system
  • The step size is chosen as $\alpha_k = \arg\min_{\alpha > 0} |F(x_k - \alpha J(x_k)^{-1} F(x_k))|$
  • Common line search conditions include the Armijo condition and the Wolfe conditions
• Trust region methods restrict the Newton step to a trusted region around the current iterate where the linear approximation is accurate
  • The trust region subproblem is $\min_{d} |J(x_k) d + F(x_k)|$ subject to $|d| \leq \Delta_k$, where $\Delta_k$ is the trust region radius
  • The trust region radius is updated based on the ratio of actual to predicted reduction in the norm of the system

Convergence and stability of algorithms

Conditioning and error analysis

• The convergence of numerical algorithms for polynomial system solving depends on factors such as the conditioning of the system, the multiplicity of the roots, and the choice of initial guesses
• The condition number of a polynomial system measures its sensitivity to perturbations in the coefficients or the solutions
  • Systems with high condition numbers are ill-conditioned and may lead to numerical instability
  • The condition number can be estimated using the Jacobian matrix or the higher-order derivatives of the system
• Backward error analysis measures the stability of a numerical solution by computing the smallest perturbation of the input data that makes the solution exact
  • The backward error of a computed solution $\tilde{x}$ is defined as $\eta(\tilde{x}) = \min{\varepsilon : F(x) + \Delta F(x) = 0, |\Delta F| \leq \varepsilon|F|, x = \tilde{x}}$
  • A small backward error indicates that the computed solution is the exact solution of a nearby problem
• Forward error analysis measures the difference between the computed solution and the true solution
  • The forward error of a computed solution $\tilde{x}$ is defined as $|x - \tilde{x}|$, where $x$ is the true solution
  • The forward error can be bounded in terms of the backward error and the condition number of the system

Convergence analysis and bifurcation

• The multiplicity of a root affects the convergence rate of iterative methods like Newton's method
  • For a root of multiplicity $m$, the convergence rate of Newton's method is linear with a factor of $1 - 1/m$
  • Higher-order methods or deflation techniques may be necessary for multiple roots to achieve faster convergence
• The basin of attraction of a root is the set of initial guesses that lead to convergence to that root under an iterative method
  • The size and shape of the basins can provide insight into the stability and efficiency of the method
  • Fractal boundaries between the basins can indicate sensitivity to initial guesses and potential difficulties in convergence
• Bifurcation analysis can be used to study how the solution structure of a polynomial system changes with respect to parameters
  • A bifurcation occurs when the qualitative structure of the solutions changes as a parameter varies
  • Examples of bifurcations in polynomial systems include fold bifurcations (creation or annihilation of solutions) and Hopf bifurcations (birth of periodic solutions)
  • Phenomena such as solution path jumping and divergence can occur near bifurcation points

Path-tracking strategies in homotopy methods

Predictor-corrector methods

• Path tracking is the core numerical procedure in homotopy continuation methods, responsible for following the solution paths from the start system to the target system
• Predictor-corrector methods are commonly used for path tracking, where a predictor step estimates the next point on the path and a corrector step refines it to lie on the path
  • Euler's method and higher-order Runge-Kutta methods can be used for prediction, while Newton's method is typically used for correction
    • Euler's method: $x_{k+1}^{(0)} = x_k + h \frac{dx}{dt}|_{t=t_k}$, where $h$ is the step size and $\frac{dx}{dt}$ is the tangent vector to the path
    • Runge-Kutta methods: Higher-order approximations of the path using multiple evaluations of the derivative
  • The step size $h$ and the tolerance for the corrector iteration are key parameters that control the accuracy and efficiency of path tracking
• Adaptive step size control can be used to automatically adjust the step size based on the local path geometry and the performance of the corrector iteration
  • Error-based step size control: The step size is adjusted based on the estimated local truncation error of the predictor
  • Path curvature-based step size control: The step size is adjusted based on the estimated curvature of the path to maintain a desired resolution

End games and singular solutions

• End-game strategies are used to accurately compute the solutions when approaching the end of the paths, especially near singular or multiple solutions
• Fractional power series approximation is a technique for computing singular solutions with known multiplicity structure
  • The solution path is approximated by a fractional power series $x(t) = x_0 + x_1 t^{1/p} + x_2 t^{2/p} + \ldots$, where $p$ is the multiplicity of the solution
  • The coefficients of the series are computed by solving a sequence of linear systems obtained by substituting the series into the homotopy and equating coefficients
• The Cauchy integral method is a technique for computing singular solutions without prior knowledge of the multiplicity structure
  • The solution path is represented as a Cauchy integral $x(t) = \frac{1}{2\pi i} \oint_{\Gamma} \frac{x(s)}{s-t} ds$, where $\Gamma$ is a small loop around $t$ in the complex plane
  • The integral is approximated numerically using quadrature rules and the values of $x(s)$ on the loop are obtained by solving a system of linear equations
• Path jumping can occur when the Jacobian matrix becomes ill-conditioned near singular solutions, causing the path to jump to another solution
  • Detection techniques, such as monitoring the condition number of the Jacobian or the size of the Newton step, can be used to identify potential path jumping
  • Correction techniques, such as the gamma trick (adding a small random perturbation to the homotopy) and singular endgames (using specialized methods for singular solutions), can be used to handle path jumping

Parallel and adaptive strategies

• Parallel path tracking can be used to exploit the inherent parallelism in homotopy continuation methods, where multiple paths can be tracked independently on different processors or cores
  • The paths are distributed among the available processors and tracked concurrently, with occasional communication to share information and detect completed paths
  • Load balancing strategies, such as dynamic work stealing or path redistribution, can be used to ensure efficient utilization of the parallel resources
• Adaptive strategies can be used to dynamically adjust the path-tracking parameters or switch between different path-tracking methods based on the characteristics of the problem and the performance of the current method
  • Adaptive precision control: The precision of the arithmetic (e.g., double, multiple, or arbitrary precision) is adjusted based on the numerical stability and accuracy requirements of the problem
  • Adaptive method switching: Different path-tracking methods (e.g., predictor-corrector, Runge-Kutta, or endgame methods) are used for different parts of the paths based on their suitability and efficiency
• Different path-tracking strategies can be compared in terms of their accuracy, efficiency, robustness, and scalability
  • Accuracy: The ability to compute the solutions to a desired level of precision and avoid path jumping or other numerical artifacts
  • Efficiency: The computational cost (time and memory) required to track the paths and find the solutions
  • Robustness: The ability to handle a wide range of problems, including ill-conditioned systems, singular solutions, and high-degree polynomials
  • Scalability: The ability to scale to larger problems (in terms of the number of variables and equations) and to exploit parallel computing resources