JuliaHomotopyContinuation / HomotopyContinuation.jl

A Julia package for solving systems of polynomials via homotopy continuation.
https://www.JuliaHomotopyContinuation.org
MIT License
178 stars 30 forks source link

Paths lost during tracking, `terminated_step_size_too_small` #492

Closed jkosata closed 2 years ago

jkosata commented 2 years ago

Hello, here is an example of the polynomials I am trying to solve. I know this should have (at least) 9 solutions. These are indeed found by the initial solve call, but the subsequent tracking gradually loses them.

Initialization vars = @var x1, x2, x3, x4 pars = @var a

eq1= 4.5*x1 + (-9/2)*x1*a^2 + 1.5*x1*x3^2 + 1.5*x1*x4^2 + 0.00015*x2*a + 19.5*x2^2*x1 + 0.00375*x2^3*a + 0.15*x4^2*x3 + 0.00375*x2*x1^2*a + 19.5*x1^3 - 0.05*x3^3 eq2 =4.5*x2 - 0.00015*x1*a - 0.00375*x1^3*a + (-9/2)*x2*a^2 + 19.5*x2*x1^2 + 1.5*x2*x3^2 + 1.5*x2*x4^2 - 0.15*x4*x3^2 - 0.00375*x2^2*x1*a + 19.5*x2^3 + 0.05*x4^3 eq3 = -1.5e-05 + (1/2)*x3 - 0.15*x1*x3^2 + 0.15*x1*x4^2 + 1.5*x1^2*x3 + 1.5*x2^2*x3 + (-1/2)*x3*a^2 + 1e-05*x4*a + (3/8)*x4^2*x3 + 2.5e-05*x4^3*a - 0.3*x2*x4*x3 + 2.5e-05*x4*x3^2*a + (3/8)*x3^3 eq4 = (1/2)*x4 + 1.5*x1^2*x4 - 0.15*x2*x3^2 + 0.15*x2*x4^2 + 1.5*x2^2*x4 - 1e-05*x3*a - 2.5e-05*x3^3*a + (-1/2)*x4*a^2 + (3/8)*x4*x3^2 + 0.3*x1*x4*x3 - 2.5e-05*x4^2*x3*a + (3/8)*x4^3

eqs = [eq1, eq2, eq3, eq4] s = System(equations, variables=[vars...], parameters = [a]) range = [[x] for x in LinRange(1.0,1.02, 20)]

Solving The starting parameter is not random but anyway recovers 9 solutions: warmup_a = 1.0 warmup = solve(s, target_parameters=range[warmup_a]) ---> 9 non-singular solutions (1 real)

This should now track the 9 solutions: all = solve(e, solutions(warmup), start_parameters=[warmup_a], target_parameters=range)

However, the solutions are gradually 'lost' : length.(solutions.(getindex.(all,1))) ---> 9, 5, 5, 3, 3, 3, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0

Inspecting the tracking: path_results(all[2][1]) shows that the lost paths have return_code → :terminated_step_size_too_small

Setting custom TrackerOptions with much smaller min_step_size and higher max_steps does not help.

Tracking each point separately Works fine but takes longer solve(s, start_system=:total_degree, target_parameters=range) ---> 9 non-singular solutions for each parameter

I will be glad for any hints as to what might be wrong.

PBrdng commented 2 years ago

Hi, the problem in your example is that the complement of the discriminant around a=1 is too small to meaningfully track solutions.

The discriminant consists of polynomial systems which have a singular solutions. This is the numerically infeasible locus, meaning that we can't track solutions through the discriminant. Over the real numbers the discriminant is of codimension 1, over the complex numbers it is of codimension 2. This is why we never track through the real numbers (and the reason why you should always choose complex start parameters).

I made a change of variables in a to blow up the region around a in which tracking works. When I do this, I can solve all your parameters:

using HomotopyContinuation

#Initialization
vars = @var x1, x2, x3, x4
pars = @var a

eq1= 4.5*x1 + (-9/2)*x1*a^2 + 1.5*x1*x3^2 + 1.5*x1*x4^2 + 0.00015*x2*a + 19.5*x2^2*x1 + 0.00375*x2^3*a + 0.15*x4^2*x3 + 0.00375*x2*x1^2*a + 19.5*x1^3 - 0.05*x3^3
eq2 =4.5*x2 - 0.00015*x1*a - 0.00375*x1^3*a + (-9/2)*x2*a^2 + 19.5*x2*x1^2 + 1.5*x2*x3^2 + 1.5*x2*x4^2 - 0.15*x4*x3^2 - 0.00375*x2^2*x1*a + 19.5*x2^3 + 0.05*x4^3
eq3 = -1.5e-05 + (1/2)*x3 - 0.15*x1*x3^2 + 0.15*x1*x4^2 + 1.5*x1^2*x3 + 1.5*x2^2*x3 + (-1/2)*x3*a^2 + 1e-05*x4*a + (3/8)*x4^2*x3 + 2.5e-05*x4^3*a - 0.3*x2*x4*x3 + 2.5e-05*x4*x3^2*a + (3/8)*x3^3
eq4 = (1/2)*x4 + 1.5*x1^2*x4 - 0.15*x2*x3^2 + 0.15*x2*x4^2 + 1.5*x2^2*x4 - 1e-05*x3*a - 2.5e-05*x3^3*a + (-1/2)*x4*a^2 + (3/8)*x4*x3^2 + 0.3*x1*x4*x3 - 2.5e-05*x4^2*x3*a + (3/8)*x4^3

eqs = [eq1, eq2, eq3, eq4]
eqs = [subs(e, a => 1e-12 * a + 1) for e in eqs] # blow up region around a
s = System(eqs, variables=[vars...], parameters = [a])
range = [[1e12*(x-1)] for x in LinRange(1.0,1.02, 20)] # change range accordingly

#Solving
warmup_a = [randn(ComplexF64)]
warmup = solve(s, target_parameters=warmup_a)

#This should now track the 9 solutions:
all_sols = solve(s,
            solutions(warmup),
            start_parameters=warmup_a, target_parameters=range)

Does this help?

jkosata commented 2 years ago

Hi, thanks this did the trick! I wrongly assumed that having found the correct warmup solutions was sufficient, but now I see why tracking through the real numbers is undesirable.