Solving general systems

# Solving general polynomial systems

The `solve` function is the most convenient way to solve general polynomial systems. For the mathematical background take a look at our introduction guide.

``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 and at the end we list all possible options.

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 `d1⋯dn`, where `di` is the degree of the `i`th entry of `F`. `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.)

Example

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 entries `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`.
• `AbstractSystem` (the system has to represent a homogeneous polynomial system.)

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

Example

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

Homogeneous Systems

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

Examples

``````julia> @polyvar x y z;
julia> solve([x^2+y^2+z^2, 2x+3y-z])
Result with 2 solutions
==================================
• 2 non-singular solutions (0 real)
• 0 singular solutions (0 real)
• 2 paths tracked
• random seed: 291729``````

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

Multihomogeneous Systems

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

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

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

\[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]]
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.AbstractPolynomial`s 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 `AbstractSystem`s.
• `affine_tracking::Bool=true`: Indicate whether path tracking should happen in affine space.
• `projective_tracking::Bool=false`: Indicate whether path tracking should happen in projective space. The flag `affine_tracking` is dominant.
• `path_jumping_check::Bool=true`: Enable a check whether one of the paths jumped to another one.

Path tracking specific options:

• `accuracy=1e-7`: The accuracy required during the path tracking.
• `corrector::AbstractCorrector`: The corrector used during in the predictor-corrector scheme. The default is `NewtonCorrector`.
• `initial_step_size=0.1`: The size of the first step.
• `max_corrector_iters=2`: The maximal number of correction steps in a single step.
• `max_lost_digits::Real`: The tracking is terminated if we estimate that we loose more than `max_lost_digits` in the linear algebra steps. This threshold depends on the `precision` argument.
• `max_refinement_iters=5`: The maximal number of correction steps used to refine the final value.
• `max_steps=1_000`: The maximal number of iterations the path tracker has available. Note that this changes to `10_000` for parameter homotopies.
• `max_step_size=Inf`: The maximal step size.
• `min_step_size=1e-14`: The minimal step size.
• `precision::PrecisionOption=PRECISION_FIXED_64`: The precision used for evaluating the residual in Newton's method.
• `predictor::AbstractPredictor`: The predictor used during in the predictor-corrector scheme. The default is `Heun`()`.
• `refinement_accuracy=1e-8`: The precision required for the final value.
• `simple_step_size_alg=false`: Use a more simple step size algorithm.
• `steps_jacobian_info_update::Int=1`: Every n-th step a linear system will be solved using a QR factorization to obtain an estimate for the condition number of the Jacobian.
• `terminate_ill_conditioned::Bool=true`: Indicates whether the path tracking should be terminated for ill-conditioned paths. A path is considerd ill-conditioned if the condition number of the Jacobian is larger than ≈1e14 or if it is larger than 1e`max_lost_digits`.

Endgame specific options:

• `accuracy_eg::Float64=min(accuracy, 1e-5))`: It is possible to change the accuracy during the path tracking. Usually you want lower the accuracy.
• `cond_eg_start::Float64=1e4`: The endgame is only started if the condition of the Jacobian is larger than this threshold.
• `max_winding_number::Int=12`: This limits the maximal number of loops taken in applying Cauchy's formula.
• `min_cond_at_infinity::Float64=1e7`: A path is declared as going to infinity only if it's Jacobian is also larger than this threshold.
• `samples_per_loop::Int=12`: To compute singular solutions Cauchy's integral formula is used. The accuracy of the solutions increases with the number of samples per loop.
• `t_eg_start::Float64=0.1`: The endgame starts only if `t` is smaller than this threshold.
• `tol_val_inf_accurate::Float64=1e-4`: A valuation which would result in a path declared as going to infinity is only accepted if the estimated accuracy of the valuation is less than this threshold.
• `tol_val_finite_accurate::Float64=1e-3`: A valuation which would result in a proper solution is only accepted if the estimated accuracy of the valuation is less than this threshold. This is only affects solutions where the path has at some point near 0 a condition number larger than `cond_eg_start`.

It is recommended to also take a look at the `PathTracker` documentation for some context.

Overdetermined system specific options:

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

## Result

A call to `solve` returns a `Result`:

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

The nonsingular solutions are obtained as follows.

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

Return all `PathResult`s 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`.

The singular solutions are returned by using the following.

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

Return all `PathResult`s 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`.

In order to analyse a `Result` we provide the following additional helper functions

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

Return all `PathResult`s 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 `PathResult`s 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 `Vector`s) 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 `SVector`s) 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 `Vector`s 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 `PathResult`s 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 `PathResult`s whose solutions appear with multiplicities greater 1 in 'V'. Two solutions are regarded as equal, when their pairwise distance is less than 'tol'.

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

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

Also make sure to check the documentation for `PathResult`.

## Estimate the complexity

We provide methods to compute the maximal number of solutions of polynomial systems.

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