Solving Polynomial Systems

At the heart of the package is the solve function. It takes a bunch of different input combinations and returns an AffineResult or ProjectiveResult depending on the input.

The solve function works in 3 stages.

It takes the input and constructs a homotopy $H(x,t)$ such that $H(x,1)=G(x)$ and $H(x,0)=F(x)$ as well as start solutions $\mathcal{X}$ where for all $x_1 ∈ \mathcal{X}$ we have $H(x_1, 1) ≈ 0$. This step highly depends on the concrete input you provide.
We now can Then all start solutions $x(1) = x_1 ∈ \mathcal{X}$ will be tracked to solutions $x({t_e})$ such that $H(x({t_e}), t_e) ≈ 0$ using a predictor-corrector scheme where $t_e$ is a value between $0$ and $1$ (by default $0.1$).
From these intermediate solutions $x({t_e})$ we start the so called endgame. This is an algorithm to predict the value $x(0)$ for each intermediate solution.

The reason for step 3 is that the final solutions can be singular which provides significant challenges for our gradient based predictor-corrector methods. For more background also check the FAQ.

The solve function

HomotopyContinuation.solve — Function.

solve(F; options...)

Solve the system F using a total degree homotopy. F can be

Vector{<:MultivariatePolynomials.AbstractPolynomial} (e.g. constructed by @polyvar)
Systems.AbstractSystem (the system has to represent a homogenous polynomial system.)

Example

Assume we want to solve the system $F(x,y) = (x^2+y^2+1, 2x+3y-1)$.

@polyvar x y
solve([x^2+y^2+1, 2x+3y-1])

If you polynomial system is already homogenous, but you would like to consider it as an affine system you can do

@polyvar x y z
solve([x^2+y^2+z^2, 2x+3y-z], homvar=z)

This would result in the same result as solve([x^2+y^2+1, 2x+3y-1]).

To solve $F$ by a custom Systems.AbstractSystem you can do

@polyvar x y z
# The system `F` has to be homgoenous system
F = Systems.SPSystem([x^2+y^2+z^2, 2x+3y-z]) # Systems.SPSystem <: Systems.AbstractSystem
# To solve the original affine system we have to tell that the homogenization variable has index 3
solve(F, homvar=3)

or equivalently (in this case) by

solve([x^2+y^2+z^2, 2x+3y-z], system=Systems.SPSystem)

Start Target Homotopy

solve(G, F, start_solutions; options...)

Solve the system F by tracking the each provided solution of G (as provided by start_solutions).

Example

@polyvar x y
G = [x^2-1,y-1]
F = [x^2+y^2+z^2, 2x+3y-z]
solve(G, F, [[1, 1], [-1, 1]])

Parameter Homotopy

solve(F::Vector{<:MultivariatePolynomials.AbstractPolynomial}, parametervariables, startparameters, targetparameters, startsolutions)

Construct a parameter homotopy $H(x,t) := F(x; tp₁+(1-t)p₀)$ where p₁ is startparameters, p₀ is targetparameters and the parametervariables are the variables of F which should be considerd parameters.

Example

We want to solve a parameter homotopy $H(x,t) := F(x; t[1, 0]+(1-t)[2, 4])$ where

\[F(x; a) := (x₁^2-a₁, x₁x₂-a₁+a₂)\]

and let's say we are only intersted in tracking of $[1,1]$. This can be accomplished as follows

@polyvar x[1:2] a[1:2]
F = [x[1]^2-a[1], x[1]*x[2]-a[1]+a[2]]
p₁ = [1, 0]
p₀ = [2, 4]
startsolutions = [[1, 1]]
solve(F, a, p₁, p₀, startsolutions)

Abstract Homotopy

solve(H::Homotopies.AbstractHomotopy, start_solutions; options...)

Solve the homotopy H by tracking the each solution of $H(⋅, t)$ (as provided by start_solutions) from $t=1$ to $t=0$. Note that H has to be a homotopy between homogenous polynomial systems. If it should be considered as an affine system indicate which is the index of the homogenization variable, e.g. solve(H, startsolutions, homvar=3) if the third variable is the homogenization variable.

Options

General options:

system::Systems.AbstractSystem: A constructor to assemble a Systems.AbstractSystem. The default is Systems.FPSystem. This constructor is only applied to the input of solve. The constructor is called with system(polynomials, variables) where polynomials is a vector of MultivariatePolynomials.AbstractPolynomials and variables determines the variable ordering.
homotopy::Systems.AbstractHomotopy: A constructor to construct a Homotopies.AbstractHomotopy. The default is StraightLineHomotopy. The constructor is called with homotopy(start, target) where start and target are homogenous Systems.AbstractSystems.
seed::Int: The random seed used during the computations.
homvar::Union{Int,MultivariatePolynomials.AbstractVariable}: This considers the homogenous system F as an affine system which was homogenized by homvar. If F is an AbstractSystem homvar is the index (i.e. Int) of the homogenization variable. If F is an AbstractVariables (e.g. created by @polyvar x) homvar is the actual variable used in the system F.
endgame_start=0.1: The value of t for which the endgame is started.
report_progress=true: Whether a progress bar should be printed to STDOUT.
threading=true: Enable or disable multi-threading.

Pathtracking specific:

corrector::Correctors.AbstractCorrector: The corrector used during in the predictor-corrector scheme. The default is Correctors.Newton.
corrector_maxiters=2: The maximal number of correction steps in a single step.
predictor::Predictors.AbstractPredictor: The predictor used during in the predictor-corrector scheme. The default is Predictors.RK4.
refinement_maxiters=corrector_maxiters: The maximal number of correction steps used to refine the final value.
refinement_tol=1e-11: The precision used to refine the final value.
tol=1e-7: The precision used to track a value.
initial_steplength=0.1: The initial step size for the predictor.
steplength_increase_factor=2.0: The factor with which the step size is increased after steplength_consecutive_successes_necessary consecutive successes.
steplength_decrease_factor=inv(increase_factor): The factor with which the step size is decreased after a step failed.
steplength_consecutive_successes_necessary=5: The numer of consecutive successes necessary until the step size is increased by steplength_increase_factor.
maximal_steplength=max(0.1, initial_steplength): The maximal step length.
minimal_steplength=1e-14: The minimal step size. If the size of step is below this the path is considered failed.

Endgame specific options

cauchy_loop_closed_tolerance=1e-3: The tolerance for which is used to determine whether a loop is closed. The distance between endpoints is normalized by the maximal difference between any point in the loop and the starting point.
cauchy_samples_per_loop=6: The number of samples used to predict an endpoint. A higher number of samples should result in a better approximation. Note that the error should be roughly $t^n$ where $t$ is the current time of the loop and $n$ is cauchy_samples_per_loop.
egtol=1e-10: This is the tolerance necessary to declare the endgame converged.
maxnorm=1e5: If our original problem is affine we declare a path at infinity if the infinity norm with respect to the standard patch is larger than maxnorm.
maxwindingnumber=15: The maximal windingnumber we try to find using Cauchys integral formula.
max_extrapolation_samples=4: During the endgame a Richardson extrapolation is used to improve the accuracy of certain approximations. This is the maximal number of samples used for this.
minradius=1e-15: A path is declared false if the endgame didn't finished until then.
sampling_factor=0.5: During the endgame we approach $0$ by the geometric series $h^kR₀$ where $h$ is sampling_factor and R₀ the endgame start provided in runendgame.

Solving Polynomial Systems

The solve function

The result of solve

PathResult