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

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.

We now can 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$).

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

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

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

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

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

.

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

`HomotopyContinuation.Solving.results`

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

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).
realsolutions = results(R, onlyreal=true, onlynonsingular=true)
```

`HomotopyContinuation.Solving.mapresults`

— Function.`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).
realsolutions = results(solution, R, onlyreal=true)
```

`HomotopyContinuation.Solving.solutions`

— Function.`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]]
```

`HomotopyContinuation.Solving.realsolutions`

— Function.`realsolutions(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 `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]]
```

`HomotopyContinuation.Solving.uniquesolutions`

— Function.`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]]
```

`HomotopyContinuation.Solving.finite`

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

.

`Base.real`

— Method.`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'.

`HomotopyContinuation.Solving.atinfinity`

— Function.`atinfinity(result::AffineResult)`

Get all results where the solutions is at infinity.

`HomotopyContinuation.Solving.singular`

— Function.`singular(result::Results; conditions...)`

Return all `PathResult`

s for which the solution is singular. This is just a shorthand for `results(R; onlysingular=true, conditions...)`

. For the possible `conditions`

see `results`

.

`HomotopyContinuation.Solving.nonsingular`

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

Return all `PathResult`

s for which the solution is non-singular. This is just a shorthand for `results(R; onlynonsingular=true, conditions...)`

. For the possible `conditions`

see `results`

.

`HomotopyContinuation.Solving.failed`

— Function.`failed(result)`

Get all results where the path tracking failed.

`HomotopyContinuation.Solving.multiplicities`

— Function.`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 `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'.

`HomotopyContinuation.Solving.seed`

— Function.`seed(result)`

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.

`HomotopyContinuation.Solving.nresults`

— Function.`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)
```

`HomotopyContinuation.Solving.nfinite`

— Function.`nfinite(affineresult)`

The number of finite solutions.

`HomotopyContinuation.Solving.nreal`

— Function.`nreal(result; tol=1e-6)`

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

.

`HomotopyContinuation.Solving.nsingular`

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

.

`HomotopyContinuation.Solving.nnonsingular`

— Function.`nnonsingular(result; tol=1e-10)`

The number of non-singular solutions.

`HomotopyContinuation.Solving.natinfinity`

— Function.`natinfinity(affineresult)`

The number of solutions at infinity.

`HomotopyContinuation.Solving.nfailed`

— Function.`nafailed(result)`

The number of failed paths.

### PathResult

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

`HomotopyContinuation.Solving.solution`

— Function.`solution(pathresult)`

Get the solution of the path.

`HomotopyContinuation.Solving.residual`

— Function.`residual(pathresult)`

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

`HomotopyContinuation.Solving.startsolution`

— Function.`startsolution(pathresult)`

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

`Base.isreal`

— Method.```
isreal(pathresult; tol=1e-6)
isreal(pathresult, tol)
```

Determine whether infinity norm of the imaginary part of the so

`LinearAlgebra.issuccess`

— Method.`issuccess(pathresult)`

Checks whether the path is successfull.

`HomotopyContinuation.Solving.isfailed`

— Function.`isfailed(pathresult)`

Checks whether the path failed.

`HomotopyContinuation.Solving.isaffine`

— Function.`isaffine(pathresult)`

Checks whether the path result is affine.

`HomotopyContinuation.Solving.isprojective`

— Function.`isprojective(pathresult)`

Checks whether the path result is projective.

`HomotopyContinuation.Solving.isatinfinity`

— Function.`isatinfinity(pathresult)`

Checks whether the path goes to infinity.

`HomotopyContinuation.Solving.issingular`

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

.

`HomotopyContinuation.Solving.isnonsingular`

— Function.`isnonsingular(pathresult; tol=1e10)`

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

## Monodromy Solve

`HomotopyContinuation.Monodromy.monodromy_solve`

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

- f1 is applied on the original solution
`s`

- f2 is applied on
`s`

and the results of 1 - 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)
```

### Strategies

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