Paths lost during tracking, `terminated_step_size_too_small`
jkosata opened this issue · 2 comments
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.
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?
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.