Solving Polynomial Systems

solve(args...; options...)::Result

The solve function takes many different arguments and options depending on your specific situation, but in the it always returns a Result containing the result of the computations. In the following we show the different inputs solve takes.

Total Degree Homotopy

solve(F; options...)

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

• Vector{<:MultivariatePolynomials.AbstractPolynomial} (e.g. constructed by @polyvar)
• A composition of polynomial systems constructed by compose.
• AbstractSystem (the system has to represent a homogeneous polynomial system.)

Additionally if F has a multi-homogenous structure you can provide variable groups to use a multi-homogenous totaldegree homotopy.

Examples

We can solve the system $F(x,y) = (x^2+y^2+1, 2x+3y-1)$ in the following way:

julia> @polyvar x y;
julia> solve([x^2+y^2+1, 2x+3y-1])
Result with 2 tracked paths
==================================
• 2 non-singular finite solutions (0 real)
• 0 singular finite solutions (0 real)
• 0 solutions at infinity
• 0 failed paths
• random seed: 448703

If your polynomial system is already homogeneous, 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 yields the same result as solve([x^2+y^2+1, 2x+3y-1]).

By exploiting the multi-homogenous structure of a polynomial system it is possible to decrease the number of paths necessary to track.

@polyvar x y
# Use variable groups to only track 2 paths instead of 4
solve([x*y - 6, x^2 - 5], variable_groups=[(x,), (y,)])

To check whether a certain variable grouping is beneficial you can use the bezout_number function.

Start Target Homotopy

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

This constructs the homotopy $H(x,t) = tG(x)+(1-t)F(x)$ to compute solutions of the system F. start_solutions is a list of solutions of G which are tracked to solutions of F.

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, startsolutions; parameters, start_parameters, target_parameters, start_gamma=nothing, target_gamma=nothing)

Solve the parameter homotopy

where $p₁$ (=start_parameters) and $p₀$ (=target_parameters) are vectors of parameter values for $F$ and $γ₁$ (=start_gamma) and $γ₀$ (=target_gamma) are complex numbers. If start_parameters or target_parameters is nothing, it is assumed that γ₁ and γ₀ are $1$. The input parameters specifies the variables of F which should be considered as parameters. Necessarily we have length(parameters) == length(p₁) == length(p₀).

solve(F, startsolutions; parameters, p₁, p₀, γ₁=nothing, γ₀=nothing)

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

Example

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

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]]
p₁ = [1, 0]
p₀ = [3im, 0.5+2im]
solve(F, startsolutions; parameters=a, start_parameters=p₁, target_parameters=p₀)
# If you like unicode this is also possible
solve(F, startsolutions; parameters=a, p₁=p₁, p₀=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 homogeneous 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:

• seed::Int: The random seed used during the computations.
• show_progress=true: Whether a progress bar should be printed to standard out.
• threading=true: Enable or disable multi-threading.
• path_result_details=:default: The amount of information computed in each path result. Possible values are :minimal (minimal details), :default (default) and :extensive (all information possible).
• homvar::Union{Int,MultivariatePolynomials.AbstractVariable}: This considers the homogeneous 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.
• system::AbstractSystem: A constructor to assemble a AbstractSystem. The default is SPSystem. 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. If you experience significant compilation times, consider to change system to FPSystem.
• homotopy::AbstractHomotopy: A constructor to construct a AbstractHomotopy for the totaldegree and start target homotopy. The default is StraightLineHomotopy. The constructor is called with homotopy(start, target) where start and target are homogeneous AbstractSystems.
• affine_tracking::Bool=false: Indicate whether path tracking should happen in affine space rather than projective space. Currently this is only supported for parameter homotopies.
• path_jumping_check::Bool=true: Enable a check whether one of the paths jumped to another one.

Path tracking specific options:

• corrector::AbstractCorrector: The corrector used during in the predictor-corrector scheme. The default is NewtonCorrector.
• max_corrector_iters=3: The maximal number of correction steps in a single step.
• initial_step_size=0.1: The step size of the first step.
• max_steps=1_000: The maximal number of iterations the path tracker has available.
• min_step_size=1e-14: The minimal step size.
• max_step_size=Inf: The maximal step size.
• maximal_lost_digits::Real=-(log₁₀(eps) + 3): The tracking is terminated if we estimate that we loose more than maximal_lost_digits in the linear algebra steps.
• predictor::AbstractPredictor: The predictor used during in the predictor-corrector scheme. The default is Heun()`.
• max_refinement_iters=10: The maximal number of correction steps used to refine the final value.
• refinement_accuracy=1e-8: The precision used to refine the final value.
• accuracy=1e-7: The precision used to track a value.
• auto_scaling=true: This only applies if we track in affine space. Automatically regauges the variables to effectively compute with a relative accuracy instead of an absolute one.
• overdetermined_min_accuracy=1e-5: The minimal accuracy a non-singular solution needs to have to be considered a solution of the original system.
• overdetermined_min_residual=1e-3: The minimal residual a singular solution needs to have to be considered a solution of the original system.

Endgame specific options:

• at_infinity_check::Bool=true: Whether the path tracker should stop paths going to infinity early.
• min_step_size_endgame_start=1e-10: The endgame only starts if the step size becomes smaller that the provided value.
• samples_per_loop::Int=5: To compute singular solutions Cauchy's integral formula is used. The accuracy of the solutions increases with the number of samples per loop.
• max_winding_number::Int=12: The maximal number of loops used in Cauchy's integral formula.
• max_affine_norm::Float64=1e6: A fallback heuristic to decide whether a path is going to infinity.
• min_val_accuracy::Float64=0.001: A tolerance used to decide whether we are in the endgame zone.

Result{V<:AbstractVector}

The result of solve. This is a wrapper around the results of each single path (PathResult) and it contains some additional informations like a random seed to replicate the result.

seed(result)

The random seed used in the computation.

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

Return all PathResults for which the given conditions apply.

Example

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.

Example

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

mapresults(f, result::MonodromyResult; onlyreal=false, realtol=1e-6)

Apply the function f to all entries of MonodromyResult for which the given conditions apply.

Example

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

solutions(result; conditions...)

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

Example

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]]

solutions(loop::Loop)

Get the solutions of the loop.

solutions(result::MonodromyResult; onlyreal=false, realtol=1e-6)

Return all solutions (as SVectors) for which the given conditions apply.

Example

realsolutions = solutions(R, onlyreal=true)

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.

Example

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 less than 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.

Example

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

atinfinity(result::AffineResult)

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.

failed(result)

Get all results where the path tracking failed.

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

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

The number of solutions which satisfy the corresponding predicates.

Example

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)

nfinite(result)

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 less than 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.

natinfinity(result)

The number of solutions at infinity.

nafailed(result)

The number of failed paths.