The solve function

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

The solve function takes many different arguments and options depending on your specific situation, but in the end it always returns a Result containing the result of the computations. Depending on the prodived arguments different kind of homotopies are constructed. In particular it is possible to construct the following homotopies:

• Total degree homotopy
• Polyhedral homotopy
• Parameter homotopy
• Multi-homogeneous homotopy
• Start target homotopy

If the input is a homogeneous polynomial system, solutions in projective space are computed. Otherwise affine solutions are computed.

Options

The solve routines takes many different options. In particular all options to CoreTracker and PathTracker are allowed. Additionally the following options are allowed:

• affine_tracking::Bool=true: Indicate whether path tracking should happen in affine space.
• early_stop_cb: Here it is possible to provide a function (or any callable struct) which accepts a PathResult r as input and returns a Bool. If early_stop_cb(r) is true then no further paths are tracked and the computation is finished. This is only called for successfull paths unless save_all_paths is true. This is for example useful if you only want to compute one solution of a polynomial system. For this early_stop_cb = _ -> true would be sufficient.
• 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.
• 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.
• path_jumping_check::Bool=true: Enable a check whether one of the paths jumped to another one.
• path_result_details=:default: The amount of information computed in each path result. Possible values are :minimal (minimal details), :default (default) and :extensive.
• projective_tracking::Bool=false: Indicate whether path tracking should happen in projective space. The flag affine_tracking is dominant.
• seed: The random seed used during the computations. The seed is also reported in the result. For a given random seed the result is always identical.
• show_progress (default true): Whether a progress bar should be printed to report the progress of the current computation.
• 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.
• system_scaling (default :equations) Whether to apply an automatic scaling of the equations (:equations), of the equations and variables (:equations_and_variables) or no scaling at all (nothing).
• threading (default true): Enable or disable multi-threading. The number of threads used is controlled by the environment variable JULIA_NUM_THREADS.
• variable_ordering: Provide a custom ordering of the variables.

Examples

Total Degree Homotopy

solve(F; options...)

Solve the system F using a start system computed from the degrees of the entries of F. The number of paths to track is equal to the total degree d₁⋯dⱼ, where dᵢ is the degree of the ith entry of F. F can be

• Vector{<:MultivariatePolynomials.AbstractPolynomial} (e.g. constructed by using the exported @polyvar)
• A composition of polynomial systems constructed by compose
• Any AbstractSystem

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 solutions
==================================
• 2 non-singular solutions (0 real)
• 0 singular solutions (0 real)
• 2 paths tracked
• random seed: 661766

Polyhedral Homotopy

solve(F; start_system = :polyhedral, only_torus=false, options...)

Solve the system F using a start system computed from the Newton Polytopes of the ntries F. The number of paths to track is equal to the mixed volume of the Newton Polytopes of the entries of F. The mixed volume is at most the total degree of F. F can be

• Vector{<:MultivariatePolynomials.AbstractPolynomial} (e.g. constructed by @polyvar)
• A composition of polynomial systems constructed by compose. Note that the composition will not preserved.

If only_torus == true then only solutions in the algebraic torus $(ℂ\setminus \{0\})^n$ will be computed.

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]; start_system = :polyhedral)
Result with 2 solutions
==================================
• 2 non-singular solutions (0 real)
• 0 singular solutions (0 real)
• 2 paths tracked
• random seed: 222880

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 γ₀=γ₁=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.

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₀)

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.

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

Homogeneous Systems

If F has is homogeneous, we return results in projective space

julia> @polyvar x y z;
julia> solve([x^2+y^2+z^2, 2x+3y-z])
Result{PVector{Complex{Float64},1}} with 2 solutions
====================================================
• 2 non-singular solutions (0 real)
• 0 singular solutions (0 real)
• 2 paths tracked
• random seed: 490575

It your polynomial system is not homogeneous, you can homogenize it as follows

@polyvar x y
g = [x^2+y^2+1, 2x+3y-1]
f = homogenize(g)

It is also possible to specify the homogenizing variable.

@polyvar x y z
g = [x^2+y^2+1, 2x+3y-1]
f = homogenize(g, z)

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

Multi-homogeneous Systems

By exploiting the multi-homogeneous 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.

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

Many parameters

solve(
    F, start_solutions;
    start_parameters = ...,
    target_parameters = ..., # multiple parameters
    options...,
    many_parameters_options...
)

Solve the system paramterized system F for many different target parameters. If no further options are specified (resp. only those which are explained above) then the result of this funciton is similar to

map(params) do p
    res = solve(F, S; start_parameters = start_params, target_parameters = p, options...)
    (res, p)
end

Note if threading = true then the parallelization is done on a parameter by parameter basis (and not as above on a path by path basis). Note that order of the results is not guaranteed. The many_parameters_options allow to apply transformations on the results

Many parameter options

• flatten::Bool: Flatten the output of transform_result. This is useful for example if transform_result returns a vector of solutions, and you only want a single vector of solutions as the result (instead of a vector of vector of solutions).
• transform_parameters::Function: Transform a parameters values p before passing it to target_parameters = .... By default this identity.
• transform_result::Function: A function taking two arguments, the result and the parameters p. By default this returns the tuple (result, p).

Many parameter examples

We start with some setup. We intersect a circle with a many different lines.

@polyvar x y
f = x^2 + y^2 - 1

@polyvar a b c
l = a * x + b * y + c
F = [f, l]

# Compute start solutions S₀ for given start parameters p₀
p₀ = randn(ComplexF64, 3)
S₀ = solutions(solve(subs(F, [a, b, c] => p₀)))
# The parameters we are intersted in
params = [rand(3) for i = 1:100]

Here are the some examples of the different options.

julia> result1 = solve(
           F,
           S₀;
           parameters = [a, b, c],
           start_parameters = p₀,
           target_parameters = params,
       );

julia> typeof(result1)
Array{Tuple{Result{Array{Complex{Float64},1}},Array{Float64,1}},1}

julia> result1[1]
(Result{Array{Complex{Float64},1}} with 2 solutions
==================================================
• 2 non-singular solutions (2 real)
• 0 singular solutions (0 real)
• 2 paths tracked
• random seed: 814544
, [0.6323567622396351, 0.15381245892007867, 0.44070720179305134])

julia> # Only keep real solutions
       result2 = solve(
           F,
           S₀;
           parameters = [a, b, c],
           start_parameters = p₀,
           target_parameters = params,
           transform_result = (r, p) -> real_solutions(r),
       );

julia> typeof(result2)
Array{Array{Array{Float64,1},1},1}

julia> result2[1:3]
3-element Array{Array{Array{Float64,1},1},1}:
 [[-0.4840916699908168, -0.8750172884266356], [-0.831904725657279, 0.554918487193468]]
 [[-0.9979849536743425, 0.06345102236859917], [0.5252178600449688, -0.8509678016762935]]
 [[-0.8602211106475295, 0.5099212103809087], [0.7844127443972155, -0.6202391848530926]]

julia> # Now instead of an Array{Array{Array{Float64,1},1},1} we want to have an
       # Array{Array{Float64,1},1}
       result3 = solve(
           F,
           S₀;
           parameters = [a, b, c],
           start_parameters = p₀,
           target_parameters = params,
           transform_result = (r, p) -> real_solutions(r),
           flatten = true
       );

julia> typeof(result3)
Array{Array{Float64,1},1}

julia> result3[1:3]
3-element Array{Array{Float64,1},1}:
 [-0.4840916699908168, -0.8750172884266356]
 [-0.831904725657279, 0.554918487193468]
 [-0.9979849536743425, 0.06345102236859917]

julia> # The passed `params` do not directly need to be the target parameters.
       # Instead they can be some more concrete informations (e.g. an index)
       # and we can them by using the `transform_parameters` method
       result4 = solve(
           F,
           S₀;
           parameters = [a, b, c],
           start_parameters = p₀,
           target_parameters = 1:100,
           transform_result = (r, p) -> (real_solutions(r), p),
           transform_parameters = _ -> rand(3)
       );

julia> typeof(result4)
Array{Tuple{Array{Array{Float64,1},1},Int64},1}

julia> result4[1:3]
3-element Array{Tuple{Array{Array{Float64,1},1},Int64},1}:
 ([[-0.9459225556889465, -0.3243925378902725], [0.39166252055983675, -0.9201089446303168]], 1)
 ([[-0.838341089344409, 0.5451460519134572], [0.4299213699880814, -0.9028663332008627]], 68)
 ([[-0.8261686468010706, -0.5634229024834612], [-0.19003485030903675, -0.9817773452611452]], 69)

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(prob::AbstractProblem)

Get the random seed used for the problem prob.

seed(result)

The random seed used in the computation.

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

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

singular(R::Results; tol=1e10, multiple_results=false, kwargs...)

Return all PathResults for which the solution is singular. A solution is labeled singular if the condition number is greater than singular_tol, or if the winding number is > 1. If multiple_results=false only one point from each cluster of multiple solutions is returned. If If multiple_results=true all singular solutions in R are returned. For the possible kwargs see results.

results(result; only_real=false, real_tol=1e-6, only_nonsingular=false,
            onlysigular=false, singular_tol=1e10, onlyfinite=true, multiple_results=false)

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).
real_solutions = results(R, only_real=true, only_nonsingular=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).
real_solutions = mapresults(solution, R, only_real=true)

mapresults(f, result::MonodromyResult; only_real=false, real_tol=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).
real_solutions = mapresults(solution, R, only_real=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(result::MonodromyResult; only_real=false, real_tol=1e-6)

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

Example

real_solutions = solutions(R, only_real=true)

solutions(MS::MonodromySolver)

Get the solutions of the loop.

real_solutions(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 only_real is always true and real_tol is now tol.

Example

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

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

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

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 is_real for details regarding the determination of 'realness'.

at_infinity(result::AffineResult)

Get all results where the solutions is at infinity.

failed(result)

Get all results where the path tracking failed.

multiplicities!(result::Result; tol=1e-6)

Compute the multiplicities of the solutions in result with respect to the given tolerance.

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;
    only_real = false,
    real_tol = 1e-6,
    only_nonsingular = false,
    singular_tol = 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, only_real=true, real_tol=1e-8, only_nonsingular=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;
    singular_tol = 1e10,
    multiplicitytol = 1e-5,
    counting_multiplicities = false,
    kwargs...,
)

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. If counting_multiplicities=true the number of singular solutions times their multiplicities is returned.

nnonsingular(result; tol=1e-10)

The number of non-singular solutions.

nat_infinity(result)

The number of solutions at infinity.

nafailed(result)

The number of failed paths.

bezout_number(F::MPPolys; variable_groups=[variables(F)], homvars=nothing, parameters=nothing)
bezout_number(multidegrees, groups::VariableGroups)

Compute the multi-homogeneous bezout number associated to the given system and variable groups.

mixed_volume(F::Vector{<:MP.AbstractPolynomialLike}; show_progress=true, algorithm=:regeneration)
mixed_volume(𝑨::Vector{<:Matrix}; show_progress=true, algorithm=:regeneration)

Compute the mixed volume of the given polynomial system F resp. represented by the support 𝑨. There are two possible values for algorithm:

• :total_degree: Use the total degree homotopy algorithm described in Section 7.1
• :regeneration: Use the tropical regeneration algorithm described in Section 7.2