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

Missed solution #502

Open laurentbartholdi opened 2 years ago

laurentbartholdi commented 2 years ago

I have a system with 84 solutions, but HomotopyContinuation returns only 83 when I use InterpretedSystem:

julia> function hc_bicrit(m,n)
    @assert 2≤m≤n-1
    @var a,b, z[2:n-1]
    cycle = [0; 1; z[1:m-2]...; :∞; z[m:n-2]...]

    f(v,w) = (v==:∞ ? a-w*b : w==:∞ ? 1+b*v^2 : (1+a*v^2)-w*(1+b*v^2))

    eqs = [f(cycle[i],cycle[(i%n)+1]) for i=1:n]
    solve(InterpretedSystem(System(eqs)))
end
julia> count(a->!(a[1]≈a[2]),solutions(hc_bicrit(7,9)))
83
julia> count(a->!(a[1]≈a[2]),solutions(hc_bicrit(2,9))) # solutions should be in bijection with previous ones
84

while replacing InterpretedSystem(System(eqs)) with System(eqs) gives the correct count 84.

I'm using the latest release, 2.6.4, in Julia master (Version 1.9.0-DEV.97, Commit 822aea53c3*).

saschatimme commented 2 years ago

I looked into this a little bit. The problem is not the difference between InterpretedSystem(System(eqs)) and System(eqs) (these should hit ultimately the same core code) but there is a randomness in the solving of the system that can yield to different results (which is ultimately a failure of the used heuristics).

Solving the system 100 times and testing for the 84 solutions show that we only end up with all 84 solutions in 56 of the instances

julia> results = [count(a->!(a[1]≈a[2]),solutions(hc_bicrit(7,9))) for i in 1:100];
julia> count(n -> n == 84, results)
56

A possible workaround (until we improve the underlying heuristics) is to solve the system multiple times, merge all solutions and then pick all unique solutions

function hc_bicrit_sols(m,n)
          s1 = solutions(hc_bicrit(m,n))
          s2 = solutions(hc_bicrit(m,n))
           s3 = solutions(hc_bicrit(m,n))
           # merge solutions and filter for unique ones
          unique_points([s1;s2;s3])
   end

With this I get

julia> results = [count(a->!(a[1]≈a[2]),(hc_bicrit_sols(7,9))) for i in 1:100]
julia> count(n -> n == 84, results)
100