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

Starting solution is reported as invalid but should be valid? #475

Open iliailmer opened 3 years ago

iliailmer commented 3 years ago

I have a system with large integers as a solution, some solutions are of 1e+150 in magnitude. When using monodromy_solve, I keep getting Warning: None of the provided solutions is a valid start solution.

Is it possible that such large integers are not properly processed by HomotopyContinuation.jl?

saschatimme commented 3 years ago

Double precision floating point range is only up to around 10^308. So even just squaring a coordinate in your case would result in an overflow (inf). Can you rescale your system such that the solutions are of reasonable magnitude?

iliailmer commented 3 years ago

Thank you for your response @saschatimme. So I changed all large coefficients to 1's in the system below:

 S=[1 - x1_0 - u₁*x2_0
 x1_1 - a*u₂*x1_0 + b*x1_0*x2_0
 c*x2_0 + u₃*x2_1 - d*x1_0*x2_0
 -u₄ - x1_1 - x2_1
 x1_2 + x1_1*(-a + b*x2_0) + b*u₅*x1_0*x2_1
 u₆*x2_2 + x2_1*(c - d*x1_0) - d*x1_1*x2_0
 -u₇ - x1_2 - x2_2
 x1_3 - a*x1_2 + b*u₈*(x1_0*x2_2 + 2*x1_1*x2_1 + x1_2*x2_0)
 c*x2_2 + d*(-x1_0*x2_2 - 2*x1_1*x2_1 - x1_2*x2_0) + u₉*x2_3
 1 - x1_3 - x2_3*u₁₀
 x1_4 - a*x1_3 + (x1_0*x2_3 + 3*x1_1*x2_2 + 3*x1_2*x2_1 + x1_3*x2_0)*b*u₁₁
 x2_4 + (-x1_0*x2_3 - 3*x1_1*x2_2 - 3*x1_2*x2_1 - x1_3*x2_0)*d + c*x2_3*u₁₂
 1 - x1_4 - x2_4*u₁₃
 x1_5 - a*x1_4 + b*(4*x1_1*x2_3 + 6*x1_2*x2_2 + 4*x1_3*x2_1 + x1_4*x2_0 + x1_0*x2_4*u₁₄)
 x2_5 + c*x2_4 + d*(-4*x1_1*x2_3 - 6*x1_2*x2_2 - 4*x1_3*x2_1 - x1_4*x2_0 - x1_0*x2_4*u₁₅)
 1 - x1_5 - x2_5*u₁₆
 x1_6 - a*x1_5 + (x1_0*x2_5 + 5*x1_1*x2_4 + 10*x1_2*x2_3 + 10*x1_3*x2_2 + 5*x1_4*x2_1 + x1_5*x2_0)*b*u₁₇
 x2_6 + c*x2_5 + d*(-5*x1_1*x2_4 - 10*x1_2*x2_3 - 10*x1_3*x2_2 - 5*x1_4*x2_1 - x1_5*x2_0 - x1_0*x2_5*u₁₈)
 -1 - x1_6 - x2_6*u₁₉
 -1 + z_aux + w_aux*u₂₀]

The system is created as:

F = System(S; parameters=u)

where u is declared as @var u[1:20]. I chose starting solutions for parameters u to be u_ = [0, 2, 0, -2, -1, 1, -2, 0, 3, 0, 0, 7, 0, -15, -13, 0, 0, -29, -2, 0] and everything else is [1 for i in 1:20]. Sanity check:

sum(S(F.variables=>[1 for i in 1:20])(u=>u_))

returns 0 (that is, the solutions are valid). But running monodromy_solve(F, [1 for i in 1:20], u_) gives

┌ Warning: None of the provided solutions is a valid start solution.
└ @ HomotopyContinuation ~/.julia/packages/HomotopyContinuation/XnHM3/src/monodromy.jl:782
MonodromyResult
===============
• return_code → :invalid_startvalue
• 0 solutions
• 0 tracked loops
• random_seed → 0xc087eb4f

Am i doing something wrong? I was following this example.

PBrdng commented 3 years ago

Hi, can you send the complete code (including variable declaration)?

saschatimme commented 3 years ago

The problem is that the Jacobian of your system at the given solution and u_ does not have full rank

julia> svdvals(jacobian(F, pt, u_))
20-element Vector{Float64}:
 51.9256162568611
 25.6244168604164
 14.839201646797381
  3.54015254081851
  2.853223531890554
  2.393072648775708
  2.2992341558184135
  2.126029876348194
  1.9211463650450664
  1.5850687817577318
  1.264164003299951
  1.151668682037086
  1.0
  0.8426994519663353
  0.7053305670837358
  0.45779806231391557
  0.3361014969100419
  0.22698447068524652
  0.03863718539889924
  3.5354215785959044e-16

where pt = [1 for i in 1:20]. A start solution is only valid if the Jacobian has full column-rank and the given start point is an approximate zero.

iliailmer commented 3 years ago

@saschatimme thank you for pointing this out, I will check the same on another system. Is there any particular rule for adding parameters u?