Fix SemialgebraicSets support for system with redundancy

[blegat]

I've tried using HomotopyContinuation as solver for extracting the solutions from the moment matrix for finding the solution of polynomial optimisation problem here: https://jump.dev/SumOfSquares.jl/v0.4.6/generated/Polynomial%20Optimization/ It didn't work because the system generated has 3 equations for 2 variables although the third equation is redundant. Just removing the check as in this PR seems to resolve the issue. It finds a solution at infinity and an excess solution but they are correctly filtered out.

[saschatimme]

The check is there to avoid the case that you have more variables than equations (in which case there are infinitely many [complex] solutions). Your test case also passes without the changes.

I think there could be a problem if the equations are homogenous, i.e., 2 homogenous equations in 3 variables is fine. Was this maybe the particular case you had?

[blegat]

Indeed, that's weird, at the moment, I have trouble trying to reproduce what I obtained because of https://github.com/JuliaHomotopyContinuation/ProjectiveVectors.jl/issues/6

[saschatimme]

Sorry I missed the notifications for the projective vectors. I just tagged a new release there.

[blegat]

Ok, I think this PR is indeed not what's needed. Here is the issue:

julia> p = -x*y + 0.4403063500907135*y^2 - 0.4410018802732259*y + 0.0002957369925164957
-xy + 0.4403063500907135y² - 0.4410018802732259y + 0.0002957369925164957

julia> V = SemialgebraicSets.algebraicset([-x - y + 1, p, -x^2 + 1.0000005149411058*y^2 - 2y + 1], solver)
Algebraic Set defined by 3 equalities
 -x - y + 1.0 = 0
 -x*y + 0.4403063500907135*y^2 - 0.4410018802732259*y + 0.0002957369925164957 = 0
 -x^2 + 1.0000005149411058*y^2 - 2.0*y + 1.0 = 0

julia> S = collect(V)
2-element Vector{Vector{Float64}}:
 [0.9997947277610838, 0.00020527223893177982]
 [-0.00027790021927181114, 1.000277144120604]

julia> V = SemialgebraicSets.algebraicset([-x - y + 1, p, -x^2 + 1.0000005149411058*y^2 - 2.003y + 1], solver)
Algebraic Set defined by 3 equalities
 -x - y + 1.0 = 0
 -x*y + 0.4403063500907135*y^2 - 0.4410018802732259*y + 0.0002957369925164957 = 0
 -x^2 + 1.0000005149411058*y^2 - 2.003*y + 1.0 = 0

julia> S = collect(V)
1-element Vector{Vector{Float64}}:
 [0.9997948787834195, 0.00020504151977317816]

julia> V = SemialgebraicSets.algebraicset([-x - y + 1, p, -x^2 + 1.0000005149411058*y^2 - 2.003y + 1.001], solver)
Algebraic Set defined by 3 equalities
 -x - y + 1.0 = 0
 -x*y + 0.4403063500907135*y^2 - 0.4410018802732259*y + 0.0002957369925164957 = 0
 -x^2 + 1.0000005149411058*y^2 - 2.003*y + 1.001 = 0

julia> S = collect(V)
Vector{Float64}[]

So the issue is that the coefficients of the polynomials of the system obtained from the moment matrix have some noise. What I do with the Groeber basis approach is that I use a loose tolerance when doing comparisons. Is there a way to do make HomotopyContinuation obtain the solutions from a noisy system ?

Looking at

julia> F = System(V)
System of length 3
 2 variables: x, y

 1.0 - 1.0*x - 1.0*y
 0.000295736992516496 - 0.441001880273226*y - 1.0*x*y + 0.440306350090714*y^2
 1.001 - 2.003*y - 1.0*x^2 + 1.00000051494111*y^2

julia> solve(F; compile=false)
Result with 0 solutions
=======================
• 4 paths tracked
• 0 non-singular solutions (0 real)
• 4 excess solutions
• random_seed: 0xa14c4063
• start_system: :polyhedral

julia> S = solve(F; compile=false)
Result with 0 solutions
=======================
• 4 paths tracked
• 0 non-singular solutions (0 real)
• 4 excess solutions
• random_seed: 0x12172890
• start_system: :polyhedral

julia> S.path_results[1].solution
2-element Vector{ComplexF64}:
 -3.6530388456922798 - 4.908176105437279im
 -1.4895170878420623 - 4.746529831928619im

julia> S.path_results[2].solution
2-element Vector{ComplexF64}:
    1.0002917134459328 + 0.00029786216451824004im
 2.3659080570625917e-5 - 0.00031371889963652147im

julia> S.path_results[3].solution
2-element Vector{ComplexF64}:
 -0.2941563998978482 - 0.3870297932317766im
  1.8724168759244917 - 0.22221001278245503im

julia> S.path_results[4].solution
2-element Vector{ComplexF64}:
 0.00033549285931581925 + 0.0016667974297138283im
      0.997986737126781 - 0.0014261704608962012im

it seems the second and fourth solution are what we want, we should just use looser tolerance for filtering the excess solution.

[PBrdng]

Hi, in general it is not possible to solve overdetermined systems (= more equations than variables) with noise. The reason is in the space of overdetermined systems only a lower dimensional subset has zeros. This is because for an overdetermined system to have zeros, the hypersurfaces defined by the single polynomials must intersect non-transversally (this is a polynomial condition on the coefficients). One approach could be to take a square subsystem, but then the question is which one. I don't know to be honest.

Maybe knowing a bit more about the structure of your system can help (for instance, finding the "closest" system having solutions to the noisy system, but I would not know how to described the locus of overdetermined systems having solutions).

[blegat]

I don't know of any structure of the systems. The solutions are then checked again: we try to solve a linear system to make the moments match with the moment of the atomic measures centered at the solutions so if there is excess solutions returned it's no big deal. I added an option to allow reconsidering excess solution in https://github.com/JuliaHomotopyContinuation/HomotopyContinuation.jl/pull/488, this should cover the use case of SumOfSquares.