Solving Polynomial Systems

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},
    startsolutions; parameters::Vector{<:MP.AbstractVariable}, p₁, p₀, γ₁=nothing, γ₀=nothing)

Solve the parameter homotopy

, 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₀)$.

solve(F::Vector{<:MultivariatePolynomials.AbstractPolynomial},
        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.

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]]
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::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-8: The precision used to refine the final value.
• tol=1e-6: 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.

AffineResult <: Result

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

ProjectiveResult <: Result

The result of a homogenous system of polynomials.

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

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

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

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

seed(result)

The random seed used in the computation.

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(affineresult)

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.

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(affineresult)

The number of solutions at infinity.

nafailed(result)

The number of failed paths.

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

solution(pathresult)

Get the solution of the path.

residual(pathresult)

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

startsolution(pathresult)

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

issuccess(pathresult)

Checks whether the path is successfull.

isfailed(pathresult)

Checks whether the path failed.

isaffine(pathresult)

Checks whether the path result is affine.

isprojective(pathresult)

Checks whether the path result is projective.

isatinfinity(pathresult)

Checks whether the path goes to infinity.

issingular(pathresult; tol=1e10)
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.

isnonsingular(pathresult; tol=1e10)

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

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.

Options

• 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-5: The tolerance with which it is decided whether two solutions are identical. The paths are tracked with tolerance tol*1e-2.
• 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.

GroupActions(actions::Function...)

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

Example

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

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

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

julia> GroupActions(f1)(3)
(9,)

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

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

Triangle(;weights=true)

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.

Petal()

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