Solving Polynomial Systems

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.

  1. 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.

  2. 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$).

  3. 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

solve(F; options...)

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

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


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 AbstractSystem you can do

@polyvar x y z
# The system `F` has to be homgoenous system
F = SPSystem([x^2+y^2+z^2, 2x+3y-z]) # SPSystem <: 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=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).


@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

    startsolutions; parameters::Vector{<:MP.AbstractVariable}, p₁, p₀, γ₁=nothing, γ₀=nothing)

Solve the parameter homotopy

\[H(x, t) = F(x, (tγ₁p₁+(1-t)γ₀p₀) / (tγ₁+(1-t)γ₀))\]

, where p₁ and p₀ are a vector of parameter values for $F$ and γ₁ and γ₀ are complex numbers. If γ₁ or γ₀ is nothing, it is assumed that γ₁ and γ₀ are $1$. The input parameters specifies the parameter variables of F which should be considered as parameters. Neccessarily, $length(parameters) == length(p₁) == length(p₀)$.

        startsolutions; parameters::Vector{<:MP.AbstractVariable},
        startparameters, targetparameters,
        startgamma=randn(ComplexF64), targetgamma=randn(ComplexF64))

This is a non-unicode variant where γ₁=start_parameters, γ₀=target_parameters, γ₁=start_gamma, γ₀=target_gamma.


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]]
startsolutions = [[1, 1]]
solve(F, startsolutions, parameters=a, p₁=p₁, p₀=p₀)
# If you don't like unicode this is also possible
solve(F, startsolutions, parameters=a, startparameters=p₁, targetparameters=p₀)

Abstract Homotopy

solve(H::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.


General options:

  • system::AbstractSystem: A constructor to assemble a AbstractSystem. The default is 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::AbstractHomotopy: A constructor to construct a AbstractHomotopy. The default is StraightLineHomotopy. The constructor is called with homotopy(start, target) where start and target are homogenous 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::AbstractCorrector: The corrector used during in the predictor-corrector scheme. The default is Newton.
  • corrector_maxiters=2: The maximal number of correction steps in a single step.
  • predictor::AbstractPredictor: The predictor used during in the predictor-corrector scheme. The default is Heun.
  • refinement_maxiters=corrector_maxiters: The maximal number of correction steps used to refine the final value.
  • refinement_tol=1e-8: 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.
  • minimal_steplength=1e-14: The minimal step size. If the size of step is below this the path is considered failed.
  • maxiters=1000: The maximal number of steps per path.

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.
  • maxiters_per_step=100: The maximal number of steps bewtween two samples.

The result of solve

Depending on the input solve returns one of the following types

AffineResult <: Result

The result of an (non-homogenous) system of polynomials.

ProjectiveResult <: Result

The result of a homogenous system of polynomials.

A Result is a wrapper around the results of each single path (PathResult) and it contains some additional informations like the used random seed for the computation.

In order to analyze a Result we provide the following helper functions

results(result; onlyreal=false, realtol=1e-6, onlynonsingular=false, onlysigular=false, singulartol=1e14, onlyfinite=true)

Return all PathResults for which the given conditions apply.


R = solve(F)

# This gives us all PathResults considered non-singular and real (but still as a complex vector).
realsolutions = results(R, onlyreal=true, onlynonsingular=true)
mapresults(f::Function, result; conditions...)

Apply the function f to all PathResults for which the given conditions apply. For the possible conditions see results.


# This gives us all solutions considered real (but still as a complex vector).
realsolutions = results(solution, R, onlyreal=true)
solutions(result; conditions...)

Return all solution (as Vectors) for which the given conditions apply. For the possible conditions see results.


julia> @polyvar x y
julia> result = solve([(x-2)y, y+x+3]);
julia> solutions(result)
[[2.0+0.0im, -5.0+0.0im], [-3.0+0.0im, 0.0+0.0im]]

Get the solutions of the loop.


Returns the solutions of the result.

realsolutions(result; tol=1e-6, conditions...)

Return all real solution (as Vectors of reals) for which the given conditions apply. For the possible conditions see results. Note that onlyreal is always true and realtol is now tol.


julia> @polyvar x y
julia> result = solve([(x-2)y, y+x+3]);
julia> realsolutions(result)
[[2.0, -5.0], [-3.0, 0.0]]
realsolutions(res::MonodromyResult; tol = 1e-6)

Returns the solutions of res whose imaginary part has norm < 1e-6.

uniquesolutions(R::Result; tol=1e-6, multiplicities=false, conditions...)

Return all unique solutions. If multiplicities is true, then all unique solutions with their correspnding multiplicities as pairs (s, m) where s is the solution and m the multiplicity are returned. For the possible conditions see results.


julia> @polyvar x;
julia> uniquesolutions([(x-3)^3*(x+2)], multiplicities=true)
[([3.0+0.0im], 3), ([-2.0+0.0im], 1)]
julia> uniquesolutions([(x-3)^3*(x+2)])
[[3.0+0.0im], [-2.0+0.0im]]
finite(result::AffineResults; conditions...)

Return all PathResults for which the solution is finite. This is just a shorthand for results(R; onlyfinite=true, conditions...). For the possible conditions see results.

real(result, tol=1e-6)

Get all results where the solutions are real with the given tolerance tol. See isreal for details regarding the determination of 'realness'.


Get all results where the solutions is at infinity.

singular(result::Results; conditions...)

Return all PathResults for which the solution is singular. This is just a shorthand for results(R; onlysingular=true, conditions...). For the possible conditions see results.

nonsingular(result::Results; conditions...)

Return all PathResults for which the solution is non-singular. This is just a shorthand for results(R; onlynonsingular=true, conditions...). For the possible conditions see results.


Get all results where the path tracking failed.

multiplicities(vectors, tol, distance)

Returns an array of arrays of integers. Each vector v in vectors contains all indices i,j such that V[i] and V[j] have distance at most tol.

multiplicities(V::Results; tol=1e-6)

Returns a Vector of Vector{PathResult}s grouping the PathResults whose solutions appear with multiplicities greater 1 in 'V'. Two solutions are regarded as equal, when their pairwise distance is less than 'tol'.


The random seed used in the computation.

If you are interested in the number of solutions of a certain kind we also provide the following helper functions.

nresults(result; onlyreal=false, realtol=1e-6, onlynonsingular=false, singulartol=1e14, onlyfinite=true)

The number of solutions which satisfy the corresponding predicates.


result = solve(F)
# Get all non-singular results where all imaginary parts are smaller than 1e-8
nresults(result, onlyreal=true, realtol=1e-8, onlynonsingular=true)

The number of finite solutions.

nreal(result; tol=1e-6)

The number of real solutions where all imaginary parts of each solution are smaller than tol.

nreal(res::MonodromyResult; tol = 1e-6)

Counts how many solutions of res have imaginary part norm < 1e-6.

nsingular(result; tol=1e10)

The number of singular solutions. A solution is considered singular if its windingnumber is larger than 1 or the condition number is larger than tol.

nnonsingular(result; tol=1e-10)

The number of non-singular solutions.


The number of solutions at infinity.


The number of failed paths.


For each path we return a PathResult containing the detailed information about the single path.

PathResult(startvalue, pathtracker_result, endgamer_result, solver)

A PathResult is the result of the tracking of a path (inclusive endgame). Its fields are

  • returncode: One of :success, :at_infinity or any error code from the EndgamerResult
  • solution::Vector{T}: The solution vector. If the algorithm computed in projective space and the solution is at infinity then the projective solution is given. Otherwise an affine solution is given if the startvalue was affine and a projective solution is given if the startvalue was projective.
  • residual::Float64: The value of the infinity norm of H(solution, 0).
  • newton_residual: The value of the 2-norm of $J_H(\text{solution})^{-1}H(\text{solution}, 0)$
  • condition_number: This is the condition number of the row-equilibrated Jacobian at the solution. A high condition number indicates a singularity.
  • windingnumber: The estimated winding number
  • angle_to_infinity: The angle to infinity is the angle of the solution to the hyperplane where the homogenizing coordinate is $0$.
  • real_solution: Indicates whether the solution is real given the defined tolerance at_infinity_tol (from the solver options).
  • startvalue: The startvalue of the path
  • iterations: The number of iterations the pathtracker needed.
  • endgame_iterations: The number of steps in the geometric series the endgamer did.
  • npredictions: The number of predictions the endgamer did.
  • predictions: The predictions of the endgamer.

The following helper functions are provided


Get the solution of the path.


Get the residual of the solution $x$ of the path, i.e., $||H(x, t)||_{\infty}$.


Get the start solution of the solution $x$ of the path.

isreal(pathresult; tol=1e-6)
isreal(pathresult, tol)

Determine whether infinity norm of the imaginary part of the so


Checks whether the path is successfull.


Checks whether the path failed.


Checks whether the path result is affine.


Checks whether the path result is projective.


Checks whether the path goes to infinity.

issingular(pathresult, tol)

Checks whether the path result is singular. This is true if the winding number > 1 or if the condition number of the Jacobian is larger than tol.


Checks whether the path result is non-singular. This is true if it is not singular.

Monodromy Solve

monodromy_solve(F, sols, p; parameters=..., options..., pathtrackerkwargs...)

Solve a polynomial system F(x;p) with specified parameters and initial solutions sols by monodromy techniques. This makes loops in the parameter space of F to find new solutions.


  • target_solutions_count=nothing: The computations are stopped if this number of solutions is reached.
  • done_callback=always_false: A callback to end the computation early. This function takes 2 arguments. The first one is the new solution x and the second one are all current solutions (including x). Return true if the compuation is done.
  • maximal_number_of_iterations_without_progress::Int=10: The maximal number of iterations (i.e. loops generated) without any progress.
  • group_action=nothing: A function taking one solution and returning other solutions if there is a constructive way to obtain them, e.g. by symmetry.
  • strategy: The strategy used to create loops. By default this will be Triangle with weights if F is a real system.
  • showprogress=true: Enable a progress meter.
  • tol::Float64=1e-7: The tolerance with which paths are tracked and with which it is decided whether two solutions are identical.
  • group_actions=GroupActions(group_action): If there is more than one group action you can use this to chain the application of them.
  • group_action_on_all_nodes=false: By default the group_action(s) are only applied on the solutions with the main parameter p. If this is enabled then it is applied for every parameter q.
  • parameter_sampler=independent_normal: A function taking the parameter p and returning a new random parameter q. By default each entry of the parameter vector is drawn independently from the unviraite normal distribution.
  • timeout=float(typemax(Int)): The maximal number of seconds the computation is allowed to run.
  • minimal_number_of_solutions: The minimal number of solutions before a stopping heuristic is applied. By default this is half of target_solutions_count if applicable otherwise 2.


If there is a group acting on the solution set of the polynomial system this can provided with the group_action keyword for single group actions or with the group_actions keyword for compositions of group actions.


Store a bunch of group actions (f1, f2, f3, ...). Each action has to return a tuple. The actions are applied in the following sense

  1. f1 is applied on the original solution s
  2. f2 is applied on s and the results of 1
  3. f3 is applied on s and the results of 1) and 2)

and so on


julia> f1(s) = (s * s,);

julia> f2(s) = (2s, -s, 5s);

julia> f3(s) = (s + 1,);

julia> GroupActions(f1)(3)

julia> GroupActions(f1,f2)(3)
(9, 18, -9, 45)

julia> GroupActions(f1,f2, f3)(3)
(9, 18, -9, 45, 10, 19, -8, 46)

A group action which returns the elementwise complex conjugated solutions.



A triangle is a loop consisting of the main node and two addtional nodes. If weights is true the edges are equipped with additional random weights. Note that this is usually only necessary for real parameters.


A petal is a loop consisting of the main node and one other node connected by two edges with different random weights.